{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Efficiency Warning: Interpolation will be slow, use setup.py!\n",
      "\n",
      "            python setup.py build_ext --inplace\n",
      "    \n",
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: pylab import has clobbered these variables: ['linalg']\n",
      "`%matplotlib` prevents importing * from pylab and numpy\n"
     ]
    }
   ],
   "source": [
    "from SimPEG import *\n",
    "import simpegDCIP as DC\n",
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cs = 25.\n",
    "hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]\n",
    "hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]\n",
    "hz = [(cs,7, -1.3),(cs,20)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mesh = Mesh.TensorMesh([hx, hy, hz], 'CCN')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "blk1 = Utils.ModelBuilder.getIndicesBlock(np.r_[-50, 75, -50], np.r_[75, -50, -150], mesh.gridCC)\n",
    "sighalf = 1e-3\n",
    "sigma = np.ones(mesh.nC)*sighalf\n",
    "sigma[blk1] = 1e-1\n",
    "sigmahomo = np.ones(mesh.nC)*sighalf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<matplotlib.collections.QuadMesh at 0x157f41d0>,\n",
       " <matplotlib.lines.Line2D at 0x15631320>)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEZCAYAAAC5AHPcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFbVJREFUeJzt3X+w5XV93/HnSyhCO9uIg5Ufu4V1eu10K20RZXWMJSaB\n2aRVcBoDTksUmbaTrYGZmmiQSQeT6UaTxlRtoT+QIgklYTQgNkBYE0vMNGSF0ki6QcBhKXsjmhAM\nNo1xF97943xv9nA959yzH8653/vj+Zi5s9/z/ny+n+/nzt1z3/f7fX/O95uqQpKko/WivicgSVqf\nTCCSpCYmEElSExOIJKmJCUSS1MQEIklqYgKRZiDJO5N8fuj1N5Kc0d+MpPkzgUhTSvKdSf5Hkq8n\neSrJbyV5zai+VbWlqg7M+PjvTnJfkm8m+S/L2v5xl7SWvv40yXNJzprlHKRhJhBpCkn+KvDfgI8A\nJwKnAR8A/nwVp7EI/BRw/fKGqrqpS1pbqmoLsBv4clU9sIrz0yZjApGm80qgquqXa+CbVbW3qh4c\n1bn76/8V3fYJSX4uyYHu7OXzSY7v2l7XndU8neR/JTl33ASq6taq+jTw1BTzfSdw41F/l9JRMIFI\n0/kS8GySG5LsSnLiUez7b4CzgNcDLwV+DHguyWkMzmp+sqpOBH4U+FSSk1YYLxMbk9OBN2IC0ZyZ\nQKQpVNU3gO8ECvjPwNeSfDrJX5u0X5IXAZcCV1TVV6rquaq6t6q+BfwT4I6quqs7xmeB+4DvX2k6\nK7T/EPCbVfX4it+Y9AKYQKQpVdVDVXVpVW0DXgWcCvzbFXY7CTge+PKIttOBt3WXr55O8jTwBuDk\nFcaceAbCIIF8YoU+0gtmApEaVNWXGPySftUKXf8I+CbwN0a0/R/gF6rqxKGvLVX1MysdflxDkjcA\npwCfXGEM6QUzgUhTSPI3k/zLrm5Bkm3A24HfnrRfVT3HYNXUh5OckuSYJK9Pchzwi8Cbk5zfxY9P\n8l1Lxxgxh2O64vuxwDFJXpzkmGXd3gF8sqr+9IV9x9LKTCDSdL4B7AR+J8n/ZZA4vgi8p2svnn9m\nMLz9o8CDwBcYrKD6aeBFVXUQuAB4P/A1Bmck72H8+/IngP8HvI9B/eTPgKuWGrvk8ja8fKVVko3y\nQKkkuxhcjz4GuK6qPtTzlCRpQ9sQCaQ7jf8S8L0MPmz1BeDtVfX7vU5MkjawjXIJ6xzg0ao6UFWH\ngF9icGlAkjQnGyWBnAY8MfT6YBeTJM3JRkkg6/86nCStM8f2PYEZWQS2Db3exuAs5C8kMclIUoOq\nGvnh1Y2SQO4DFrrnL/wBcBGDNfrLXL2KU1ruc8Cb1sj4K/VtnetqzWFU23qPLW2vhbnMKzZN2zTt\nrX2PZt/VmsO0VuMY41w9tmVDJJCqOpzk3cCvMVjG+3FXYEnSfG2IBAJQVXcCd/Y9D0naLDZKEV1H\n5Yy+J6Bvc0bfE9DznNH3BNaFDfFBwmkMiuhjn9UjSRrpng1fRJ9SX0UoWFtF9HmNZRHd2KTYNG3T\ntLf2neW+8xin72OMc8/YFi9hSZKamEAkSU1MIJKkJhbRJUkTWETvWESf71gW0QHO5eq5H/Vqzp3x\niBbR21lElyTpqJhAJElNTCCSpCYW0aUZO3fCNeNZucf/y1o1FtE7FtHnO5ZF9IF7VuEIzHhEi+jt\nLKJLknRUTCCSpCYmEElSE4vo0oxZRNfGYhG9YxF9vmNZRB+wiG4RfdYsokuSNhATiCSpiQlEktTE\nBCJJauIqLGnGXIWljcVVWB1XYc13LFdhDbgKy1VYs+YqLEnSBmICkSQ1MYFIkppYRJckTWARvWMR\nfb5jWUQ3Nik2Tds07a19Z7nvPMbp+xjjWESXJM2YCUSS1MQEIklqYhFdkjTBOiqiJ/lZ4B8C3wK+\nDFxaVX/StV0JvAt4Fri8qu7u4mcDNwDHA3dU1RWjR7eIPt+xLKIbmxSbpm2a9ta+s9x3HuP0fYxx\n1lcR/W7gb1fV3wUeBq4ESLIDuAjYAewCrkmylBWvBS6rqgVgIcmu1Z+2JG0uay6BVNXeqnque/k7\nwNZu+wLg5qo6VFUHgEeBnUlOAbZU1b6u343Ahas5Z0najNZcAlnmXcAd3fapwMGhtoPAaSPii11c\nkjRHvRTRk+wFTh7R9P6q+kzX5yrg1VX1j7rXHwPuraqbutfXAXcCB4APVtV5XfyNwHur6s3Ljllw\n+lDkJd2XJOmIr3dfSx5fW0X0pV/24yR5J/D9wPcMhReBbUOvtzI481jkyGWupfji6JEvPeq5zo5F\n9KPraxF948WmaZumvbXvLPedxzh9H2Ocq8e2rLlLWF0B/MeAC6rqm0NNtwMXJzkuyXZgAdhXVU8C\nzyTZ2RXVLwFuW/WJS9Ims+aW8QIfA44D9naLrH67qnZX1f4ktwD7gcPA7jpy/W03g2W8JzBYxnvX\n6k9bkjaXNZdAuqW449r2AHtGxO8HzpznvCRJz+cn0SVJE6yjT6LPl59En+9YFtGNTYpN0zZNe2vf\nWe47j3H6PsY46+uT6JKkdcAEIklqYgKRJDUxgUiSmrgKS5I0gauwOq7Cmu9YrsIyNik2Tds07a19\nZ7nvPMbp+xjjuApLkjRjJhBJUhMTiCSpiUV0SdIEFtE7FtHnO5ZFdGOTYtO0TdPe2neW+85jnL6P\nMY5FdEnSjJlAJElNTCCSpCYW0SVJE1hE71hEn+9YFtGNTYpN0zZNe2vfWe47j3H6PsY4FtElSTNm\nApEkNTGBSJKaWESXJE1gEb1jEX2+Y1lENzYpNk3bNO2tfWe57zzG6fsY41hElyTNmAlEktTEBCJJ\namICkSQ1cRWWJGkCV2F1XIU137FchWVsUmyatmnaW/vOct95jNP3McZxFZYkacZMIJKkJiYQSVIT\ni+iSpAnWYRE9yXuAnwVOqqo/7mJXAu8CngUur6q7u/jZwA3A8cAdVXXF6FEtos93LIvoxibFpmmb\npr217yz3ncc4fR9jnHVWRE+yDTgPeHwotgO4CNgB7AKuSbKUFa8FLquqBWAhya5VnrIkbTprMoEA\nHwbeuyx2AXBzVR2qqgPAo8DOJKcAW6pqX9fvRuDCVZupJG1Say6BJLkAOFhVX1zWdCpwcOj1QeC0\nEfHFLi5JmqNeiuhJ9gInj2i6Cng/cH5VPZPkMeA1VfVUko8B91bVTd0Y1wF3AgeAD1bVeV38jcB7\nq+rNy45ZcPpQ5CXdlyTpiK93X0seX1tF9KVf9ssleRWwHfjdrryxFbg/yU4GZxbbhrpvZXDmsdht\nD8cXRx/50hc48xfCIvrR9bWIvvFi07RN097ad5b7zmOcvo8xztVjW9bUJayq+r2qenlVba+q7QwS\nxKur6qvA7cDFSY5Lsh1YAPZV1ZPAM0l2dkX1S4DbevsmJGmTWLPLeDt/cX2tqvYnuQXYDxwGdteR\n62+7GSzjPYHBMt67VnuikrTZrOkEUlWvWPZ6D7BnRL/7gTNXa16SJD+JLkmaaB1+En0+/CT6fMey\niG5sUmyatmnaW/vOct95jNP3McZZZ59ElyStfSYQSVITE4gkqYlFdEnSBBbROxbR5zuWRXRjk2LT\ntE3T3tp3lvvOY5y+jzGORXRJ0oyZQCRJTUwgkqQmJhBJUhNXYUmSJnAVVsdVWPMdy1VYxibFpmmb\npr217yz3ncc4fR9jHFdhSZJmzAQiSWpiApEkNbGILkmawCJ6xyL6fMeyiG5sUmyatmnaW/vOct95\njNP3McaxiC5JmjETiCSpiQlEktTEBCJJarJiAknyG0n+wbLYf5rflCRJ68GKy3iTPAY8Afx6VX2g\niz1QVWetwvxmxmW8ktTihS3j/Trw3cBHk3wGuGSWU1tdLuOd71gu4zU2KTZN2zTtrX1nue88xun7\nGOO8wGW8VXW4qnYDnwI+D7xsNhOTJK1X05yB/Ieljaq6IcmDwL+Y35QkSevBigmkqv7jstf3A++a\n24wkSeuCy3glSU1MIJKkJiYQSVITb+cuSZpgnd3OPcmPALuBZ4Ffrar3dfErGRTwnwUur6q7u/jZ\nwA3A8cAdVXXF6JH9HMh8x/JzIMYmxaZpm6a9te8s953HOH0fY5zxnwNZcwkkyZuAtwB/p6oOJXlZ\nF98BXATsAE4DPptkoQanUNcCl1XVviR3JNlVVXf19T1I0mawFmsgPwz8dFUdAqiqP+ziFwA3V9Wh\nqjoAPArsTHIKsKWq9nX9bgQuXOU5S9KmsxYTyALw95Pcm+S/J3lNFz8VODjU7yCDM5Hl8cUuLkma\no14uYSXZC5w8oukqBnM6sapel+S1wC3AK1ZzfpKklfWSQKrqvHFtSX4Y+JWu3xeSPJfkJAZnFtuG\num5lcOax2G0PxxdHj/65oe0zgO1HP3lJ2tAeAw5M1XPNFdGB2xjc/feeJK8EjquqP0pyO/Bfk3yY\nwSWqBWBfVVWSZ5LsBPYxuFvwR0cP3ecqLElaD7bz/D+u19EqLOB64Prupo3fAn4IoKr2J7kF2A8c\nBnbXkQ+x7GawjPcEBst4XYElSXO25hJIt/pq5DNHqmoPsGdE/H7gzDlPTZI0ZC2uwpIkrQMmEElS\nE++FJUmaYJ3dC2t+vBfWfMfyXljGJsWmaZumvbXvLPedxzh9H2OcF/hMdEmSljOBSJKamEAkSU1M\nIJKkJiYQSVITE4gkqYkJRJLUxAQiSWpiApEkNTGBSJKamEAkSU1MIJKkJiYQSVITb+cuSZrA27l3\nvJ37fMfydu7GJsWmaZumvbXvLPedxzh9H2Mcb+cuSZoxE4gkqYkJRJLUxAQiSWpiApEkNTGBSJKa\nmEAkSU1MIJKkJiYQSVITE4gkqYkJRJLUxAQiSWpiApEkNTGBSJKarLkEkuScJPuSPJDkC0leO9R2\nZZJHkjyU5Pyh+NlJHuzaPtLPzCVpc1lzCQT4GeAnquos4F91r0myA7gI2AHsAq5JsvSQk2uBy6pq\nAVhIsmv1py1Jm8taTCBfAb6j234JsNhtXwDcXFWHquoA8CiwM8kpwJaq2tf1uxG4cBXnK0mb0pp7\npG2S04HfAopBgnt9VT2R5GPAvVV1U9fvOuBO4ADwwao6r4u/EXhvVb152bg+0laSjtoae6Rtkr3A\nySOargIuBy6vqluTvA24Hjhv9rM4A9g++2HH8pG2R9fXR9puvNg0bdO0t/ad5b7zGKfvYyx5jMHf\n5SvrJYEsnS2MkuQXq+p7u5efBK7rtheBbUNdtwIHu/jWZfFFRurzmeiStB5s5/l/XK+vZ6I/mmTp\nWtN3Aw9327cDFyc5Lsl2YAHYV1VPAs8k2dkV1S8Bblv1WUvSJtPLGcgK/hnw75O8GPiz7jVVtT/J\nLcB+4DCwu44UcHYDNwAnAHdU1V2rPmtJ2mTWXAKpqvuAnWPa9gB7RsTvB86c89QkSUPW4iUsSdI6\nYAKRJDUxgUiSmphAJElNTCCSpCYmEElSExOIJKmJCUSS1MQEIklqYgKRJDUxgUiSmphAJElNTCCS\npCYmEElSExOIJKmJCUSS1MQEIklqYgKRJDUxgUiSmqSq+p7DqkhScG7f05CkdeYeqiqjWo5d7an0\n6009Hvtzcz7+LMdvHeto9lup76T2UW3G1n5smrZp2lv7znLfeYzT9zHGuWdsi5ewJElNTCCSpCYm\nEElSExOIJKmJCUSS1MQEIklqYgKRJDUxgUiSmphAJElNTCCSpCYmEElSExOIJKlJLwkkyduS/O8k\nzyZ59bK2K5M8kuShJOcPxc9O8mDX9pGh+IuT/HIXvzfJ6av5vUjSZtXXGciDwFuB3xwOJtkBXATs\nAHYB1yRZuo3wtcBlVbUALCTZ1cUvA57q4j8PfGgV5i9Jm14vCaSqHqqqh0c0XQDcXFWHquoA8Ciw\nM8kpwJaq2tf1uxG4sNt+C/CJbvtTwPfMb+aSpCVrrQZyKnBw6PVB4LQR8cUuTvfvEwBVdRj4kyQv\nnf9UJWlzm9sDpZLsBU4e0fT+qvrMvI472eeGts8AtvczDUlasx4DDkzVc24JpKrOa9htEdg29Hor\ngzOPxW57eXxpn78O/EGSY4HvqKo/Hj18n08klKT1YDvP/+N6bT+RcPhZu7cDFyc5Lsl2YAHYV1VP\nAs8k2dkV1S8BPj20zzu67R8Afn2V5i1Jm1ovz0RP8lbgo8BJwK8meaCqvq+q9ie5BdgPHAZ2V1V1\nu+0GbgBOAO6oqru6+MeBX0jyCPAUcPEqfiuStGn1kkCq6lbg1jFte4A9I+L3A2eOiP858IOznqMk\nabK1cAlLkrQOmUAkSU1MIJKkJiYQSVITE4gkqYkJRJLUxASyKT3W9wT0bfyZrC3+PKZhAtmUDvQ9\nAX2bA31PQM9zoO8JrAsmEElSExOIJKlJjtxqamNLsjm+UUmasarKqPimSSCSpNnyEpYkqYkJRJLU\nxASygSW5OsnBJA90X9831HZlkkeSPJTk/KH42Uke7No+0s/MN48ku7qfwSNJ3tf3fDaLJAeSfLF7\nX+zrYi9NsjfJw0nuTvKSof4j3y+bnQlkYyvgw1V1Vvd1J0CSHcBFwA5gF3BN96RHgGuBy6pqAVhI\nsquPiW8GSY4B/h2Dn8EO4O1J/la/s9o0Cviu7n1xThf7cWBvVb2SwZNNfxzGvl/83YkJZDMYtXri\nAuDmqjpUVQeAR4GdSU4BtlTVvq7fjcCFqzPNTekc4NGqOlBVh4BfYvCz0epY/t54C/CJbvsTHPm/\nP+r9cg4ygWwCP5Lkd5N8fOiU/FTg4FCfg8BpI+KLXVzzcRrwxNDrpZ+D5q+Azya5L8k/7WIvr6qv\ndttfBV7ebY97v2x6vTzSVrOTZC9w8oimqxhcjvrJ7vVPAT8HXLZKU9PKXEPfnzdU1VeSvAzYm+Sh\n4caqqhU+O+bPDhPIuldV503TL8l1wGe6l4vAtqHmrQz+qlrstofjizOYpkZb/nPYxvP/0tWcVNVX\nun//MMmtDC5JfTXJyVX1ZHc592td91HvF98XeAlrQ+veBEveCjzYbd8OXJzkuCTbgQVgX1U9CTyT\nZGdXVL8EuG1VJ7253MdgocIZSY5jUKi9vec5bXhJ/nKSLd32XwHOZ/DeuB14R9ftHRz5vz/y/bK6\ns16bPAPZ2D6U5O8xON1+DPjnAFW1P8ktwH7gMLC7jtySYDdwA3ACcEdV3bXqs94kqupwkncDvwYc\nA3y8qn6/52ltBi8Hbu0WHh4L3FRVdye5D7glyWUMbsf7g7Di+2VT81YmkqQmXsKSJDUxgUiSmphA\nJElNTCCSpCYmEElSExOIJKmJCUSS1MQEIklqYgKRepLkA0muGHr9r5Nc3uecpKPhJ9GlniQ5HfiV\nqjq7e0DRw8Brq+rpnqcmTcV7YUk9qarHkzzV3a/sZOB/mjy0nphApH5dB1zK4AZ/1/c8F+moeAlL\n6lGSvwT8HoO78S54l1etJ56BSD2qqkNJfgN42uSh9cYEIvWoK56/DviBvuciHS2X8Uo9SbIDeAT4\nbFV9ue/5SEfLGogkqYlnIJKkJiYQSVITE4gkqYkJRJLUxAQiSWpiApEkNfn/ukF5nhLIpUMAAAAA\nSUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x15482e10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEZCAYAAABiu9n+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEl5JREFUeJzt3XvMZHddx/H3py21UEvaBuyNTbtIUTZeqGgLgUCJWlcT\nKHihNAEKFiRpKCQCkYvEArGmRggkCEZbbMEKbiTFom3pFrFBUdcqQrEUWrIr7NoLIqQlAvby9Y85\nm85un9/MnunMnGfneb+SJ535nu/vN7/t7Oz3me/vnJlUFZIkreWQoRcgSVq/LBKSpCaLhCSpySIh\nSWqySEiSmiwSkqQmi4TUQ5KXJ/nM2P17k5wy3IqkxbJISPtJ8qwkn03y7STfTPL3SX56rdyqOqqq\nds3xsQ9PclmSXUnuSfK5JFvHjj89yfZuXXcn2Zbk+Hk9vrQ/i4Q0Jsljgb8G3gscA5wEvB34/pKW\ncBjwNeDZVfVY4LeBbUlO7o4fDfwRcHL3cy/wp0tamzYgi4S0rycDVVV/USPfq6rtVXXzWslJHkzy\nxO72o5O8q3sX8O0kn0lyRHfs6d27k28l+fckz1lrvqr636p6e1V9rbv/N8BO4Ke6+9dV1ceq6jtV\n9V3gD4Fnzv3/gtSxSEj7+jLwQJLLk2xNckyPsX8AnAY8AzgWeCPwYJKTGL07eUdVHQO8AfhYksdN\nmzDJcYwK1380Up4NfLHHGqVeLBLSmKq6F3gWUMCfAHcn+askPzRpXJJDgFcAr6uqO6rqwar6p6r6\nP+AlwDVVdV33GDcANwG/NGXORwFXApdX1VfWOP4TwNsYFSNpISwS0n6q6taqekVVbQJ+DDgReM+U\nYY8DjgC+usaxk4Ff61pN30ryLUYtouaGc1d0Pgx8D3jNGsefBFwDvLaq/uEA/ljSTCwS0gRV9WXg\nCkbFYpL/ZvQP+pPWOPY14MNVdczYz1FV9ftrTZQkwGXA44FfqaoH9jt+MrCdUfvqyn5/Iqkfi4Q0\nJsmPJPnNbh+BJJuAc4F/nDSuqh4EPgi8O8kJSQ5N8owkhwN/BjwvyVld/IgkZ+59jDV8APhR4PlV\ntc9ZVd2YvwXeV1V//Ij+sNIBsEhI+7oXOAP45yTfYVQcvgC8vjte3Q9j9/d6A3Az8C/AN4HfAw6p\nqt3A2cBbgLsZvbN4PWu8/rp3Cb8B/CRwZ3ex3r1Jzu1SXglsBi4aO3bPI/9jS2uLXzokSWrxnYQk\nqckiIUlqskhIkposEpKkpsOGXsC8JXEnXpJmUFXZP7ZyRWLkooEe99PAc9fRY0zKfSRrndcaph1f\n69jBHBs/NvRaFhmbFO+bM0tun3HLWMOBWsa/H5NctGbUdpMkqckiIUlqskistFOGXoD2ccrQC9DD\nnDL0AtY9i8RK2zz0ArQPn4/1x+dkGouEJKnJIiFJarJISJKaLBKSpCaLhCSpySIhSWqySEiSmiwS\nkqQmi4QkqckiIUlqskhIkposEpKkJouEJKnJIiFJarJISJKaUlVDr2GukhQ8Z+hlSNJB5kaqKvtH\nDxtiKYs31JeJL+OLzOf1GI9knnl+efyk42sdM7b+Y5PifXNmyZ3HuEXPNcT809y4ZtR2kySpySIh\nSWqySEiSmiwSkqQmi4QkqckiIUlqskhIkposEpKkJouEJKnJIiFJarJISJKaLBKSpCaLhCSpySIh\nSWry+yQkSfh9Ekvh90n0z/X7JFYvNineN2eW3HmMW/RcQ8w/jd8nIUnqadAikWRXki8k+VySHV3s\n2CTbk3wlyfVJjh7Lf3OS25LcmuSs4VYuSRvD0O8kCjizqk6rqtO72JuA7VX1ZOBT3X2SbAHOAbYA\nW4H3Jxl6/ZK00tbDP7L7b5Q8H7iiu30F8ILu9tnAR6rqvqraBdwOnI4kaWGGLhIF3JDkpiSv6mLH\nVdVd3e27gOO62ycCu8fG7gZOWs4yJWljGvrspmdW1R1JHg9sT3Lr+MGqqtEprU2NY58eu30KsPkR\nLlOSVs1OYNfUrEGLRFXd0f33G0muYtQ+uivJ8VV1Z5ITgLu79D3AprHhT+hiaxjyNDJJOhhsZt9f\noNfZKbBJHpPkqO72kcBZwM3A1cB5Xdp5wMe721cDL05yeJLNwKnAjuWuWpI2liHfSRwHXJVk7zqu\nrKrrk9wEbEtyPqP3Qi8CqKpbkmwDbgHuBy6oVbtcXJLWmcGKRFXtBJ66Rvx/gJ9rjLkYuHjBS5Mk\ndfzsJkkSfnbTUvjZTf1z/eym1YtNivfNmSV3HuMWPdcQ80+zzjauJUnrn0VCktRkkZAkNVkkJElN\nFglJUpNFQpLUZJGQJDVZJCRJTV5xLUnCK66Xwiuu++d6xfXqxSbF++bMkjuPcYuea4j5p/GKa0lS\nTxYJSVKTRUKS1GSRkCQ1WSQkSU2eAitJwlNgl8JTYPvnegrs6sUmxfvmzJI7j3GLnmuI+afxFFhJ\nUk8WCUlSk0VCktRkkZAkNXl2kyQJz25aCs9u6p/r2U2rF5sU75szS+48xi16riHmn8azmyRJPVkk\nJElN7klIknBPYinck+if657E6sUmxfvmzJI7j3GLnmuI+adxT0KS1JPtJkkStpuWwnZT/1zbTasX\nmxTvmzNL7jzGLXquIeafxnaTJKkni4Qkqck9CUkS7kkshXsS/XPdk1i92KR435xZcucxbtFzDTH/\nNO5JSJJ6st0kScJ201LYbuqfa7tp9WKT4n1zZsmdx7hFzzXE/NPYbpIk9WS7SZKE7aalsN3UP9d2\n0+rFJsX75sySO49xi55riPmnsd0kSerJIiFJanJPQpKEexJL4Z5E/1z3JFYvNineN2eW3HmMW/Rc\nQ8w/jXsSkqSebDdJkrDdtBS2m/rn2m5avdikeN+cWXLnMW7Rcw0x/zS2myRJPdlukiRhu2kpbDf1\nz7XdtHqxSfG+ObPkzmPcoucaYv5pbDdJknqySEiSmtyTkGbwnMZb83m60b/HWir3JJbAPYn+uQfn\nnsRzuXHhj3ojz53zjMuKTYr3zZkldx7jFj3XEPNP456EJKkn203SDGw3afWsSLspyVbgPcChwKVV\ndcnDs2w3LXYe2022mybFJsX75sySO49xi55riPmnWYF2U5JDgfcBW4EtwLlJnjLsqiRpdU0tEkle\nm+SYZSzmAJwO3F5Vu6rqPuCjwNkDr0mSVtbUPYkkvwucA/wb8EHgkzXQRkaSXwV+oape1d1/CXBG\nVV04luOehBbOPQmtnhn3JKrqrUneBpwFvBx4X5JtwGVV9dW5r3PKcvoPOQXYPO91NLgn0T/XPYlW\nzD2JWXPnMW7Rcw0x//52ArumZh3QxnVVPZjkTuAu4AHgGOAvk9xQVW98BKvsaw+waez+JmD3w9OG\n3PyRpIPBZvb9BXrtd8cH0m56HfAy4JvApcBVVXVfkkOA26rqh+ey3gOQ5DDgy8DPAv8F7ADOraov\njeXYbtLC2W7S6pn9FNhjgV+uqv8cD3bvLp43r+UdiKq6P8lrgE8yOgX2svEC8RBPgV3sPLabbDdN\nik2K982ZJXce4xY91xDzT7P2Lz4HsifxOxOO3fIIVjSTqroWuHbZjytJG5FXXEszsN2k1bMiV1wf\nGNtNi53HdpPtpkmxSfG+ObPkzmPcoucaYv5pVuCKa0nSclkkJElN7klIknBPYinck+ife3DuSRib\nFJsU75szS+48xi16riHmn8Y9CUlST7abJEnYbloK2039c203rV5sUrxvziy58xi36LmGmH8a202S\npJ5sN0mSsN20FLab+ufablq92KR435xZcucxbtFzDTH/NLabJEk9WSQkSU3uSUiScE9iKdyT6J/r\nnsTqxSbF++bMkjuPcYuea4j5p3FPQpLUk+0mSRK2m5bCdlP/XNtNqxebFO+bM0vuPMYteq4h5p/G\ndpMkqSfbTZIkbDcthe2m/rm2m1YvNineN2eW3HmMW/RcQ8w/je0mSVJPFglJUpN7EpIk3JNYCvck\n+ue6J7F6sUnxvjmz5M5j3KLnGmL+adyTkCT1ZLtJkoTtpqWw3dQ/13bT6sUmxfvmzJI7j3GLnmuI\n+aex3SRJ6sl2kyQJ201LYbupf67tptWLTYr3zZkldx7jFj3XEPNPY7tJktSTRUKS1GSRkCQ1uXEt\nScKN66Vw47p/rhvXqxebFO+bM0vuPMYteq4h5p/GjWtJUk8WCUlSk0VCktTkxrUkCTeul8KN6/65\nblyvXmxSvG/OLLnzGLfouYaYfxo3riVJPVkkJElNFglJUpNFQpLUZJGQJDVZJCRJTV4nIUnC6ySW\nwusk+ud6ncTqxSbF++bMkjuPcYuea4j5p/E6CUlSTxYJSVKTRUKS1GSRkCQ1WSQkSU0WCUlSk0VC\nktTkxXSSJNbVxXRJLgJeCXyjC72lqq7tjr0Z+HXgAeC1VXV9F38acDlwBHBNVb2u/QheTLfYebyY\nztik2KR435xZcucxbtFzDTH/NOvrYroC3l1Vp3U/ewvEFuAcYAuwFXh/kr2V7QPA+VV1KnBqkq1D\nLFySNpIh9yQe9rYGOBv4SFXdV1W7gNuBM5KcABxVVTu6vA8BL1jOMiVp4xqySFyY5PNJLktydBc7\nEdg9lrMbOGmN+J4uLklaoIXtSSTZDhy/xqG3MmodvaO7/07gXcD583v0T4/dPgXYPL+pJWkl7AR2\nTc1aWJGoqp8/kLwklwKf6O7uATaNHX4Co3cQe7rb4/E97VmH3PyRpIPBZvb9BXodbVx3ewx7vRC4\nubt9NfDiJIcn2QycCuyoqjuBe5Kc0W1kvxT4+FIXLUkb0FDfJ3FJkqcyOstpJ/BqgKq6Jck24Bbg\nfuCCeuhCjgsYnQL7aEanwF639FVL0gYzSJGoqpdNOHYxcPEa8X8FfnyR65Ik7cuP5ZAkNVkkJElN\nFglJUpMf8CdJYl19wN/i+QF/i53HD/gzNik2Kd43Z5bceYxb9FxDzD/NOrpOQpJ0cLBISJKaLBKS\npCaLhCSpySIhSWqySEiSmiwSkqQmi4QkqckiIUlqskhIkposEpKkJouEJKnJIiFJarJISJKa/D4J\nSRJ+n8RS+H0S/XP9PonVi02K982ZJXce4xY91xDzT+P3SUiSerJISJKaLBKSpCaLhCSpySIhSWqy\nSEiSmiwSkqQmi4QkqckiIUlqskhIkposEpKkJouEJKnJIiFJarJISJKaLBKSpCaLxErbOfQCtA+f\nj/XH52Qai8RK2zX0ArSPXUMvQA+za+gFrHsWCUlSk0VCktSUqhp6DXOVZLX+QJK0JFWV/WMrVyQk\nSfNju0mS1GSRkCQ1WSRWQJKLkuxO8rnu5xfHjr05yW1Jbk1y1lj8aUlu7o69d5iVbxxJtnbPwW1J\nfmvo9WwUSXYl+UL3utjRxY5Nsj3JV5Jcn+Tosfw1Xy8bmUViNRTw7qo6rfu5FiDJFuAcYAuwFXh/\nkr0bUx8Azq+qU4FTk2wdYuEbQZJDgfcxeg62AOcmecqwq9owCjize12c3sXeBGyvqicDn+rut14v\nG/7fyA3/P2CFPOysBOBs4CNVdV9V7QJuB85IcgJwVFXt6PI+BLxgOcvckE4Hbq+qXVV1H/BRRs+N\nlmP/18bzgSu621fw0N/9tV4vp7PBWSRWx4VJPp/ksrG3zycCu8dydgMnrRHf08W1GCcBXx+7v/d5\n0OIVcEOSm5K8qosdV1V3dbfvAo7rbrdeLxvaYUMvQAcmyXbg+DUOvZVR6+gd3f13Au8Czl/S0jSd\n55kP55lVdUeSxwPbk9w6frCqasq1VRv+ubNIHCSq6ucPJC/JpcAnurt7gE1jh5/A6LejPd3t8fie\nOSxTa9v/edjEvr+xakGq6o7uv99IchWj9tFdSY6vqju71uvdXfpar5cN/7qw3bQCur/oe70QuLm7\nfTXw4iSHJ9kMnArsqKo7gXuSnNFtZL8U+PhSF72x3MTo5IBTkhzOaHP06oHXtPKSPCbJUd3tI4Gz\nGL02rgbO69LO46G/+2u+Xpa76vXHdxKr4ZIkT2X01ngn8GqAqrolyTbgFuB+4IJ66BL7C4DLgUcD\n11TVdUtf9QZRVfcneQ3wSeBQ4LKq+tLAy9oIjgOu6k7oOwy4sqquT3ITsC3J+Yw+BvZFMPX1smH5\nsRySpCbbTZKkJouEJKnJIiFJarJISJKaLBKSpCaLhCSpySIhSWqySEiSmiwS0oIl+ZnuE3p/IMmR\nSb7YfXeBtO55xbW0BEneCRzB6GNQvl5Vlwy8JOmAWCSkJUjyKEYf9Pdd4Bl+JpAOFrabpOV4HHAk\n8IOM3k1IBwXfSUhLkORq4M+BJwInVNWFAy9JOiB+VLi0YEleBny/qj6a5BDgs0nOrKq/G3hp0lS+\nk5AkNbknIUlqskhIkposEpKkJouEJKnJIiFJarJISJKaLBKSpCaLhCSp6f8Bi7dTcflLpucAAAAA\nSUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1577eba8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mesh.plotSlice(sigma, normal='X', grid=True)\n",
    "mesh.plotSlice(sigma, ind=22, normal='Z', grid=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xtemp = np.linspace(-150, 150, 21)\n",
    "ytemp = np.linspace(-150, 150, 21)\n",
    "xyz_rxP = Utils.ndgrid(xtemp-10., ytemp, np.r_[0.])\n",
    "xyz_rxN = Utils.ndgrid(xtemp+10., ytemp, np.r_[0.])\n",
    "xyz_rxM = Utils.ndgrid(xtemp, ytemp, np.r_[0.])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x15be3940>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVIAAAFRCAYAAAAmQSVBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGXFJREFUeJzt3X+QXWV9x/HPd3ezSUggG0cIgUSSpZvWTJ2BouGXlUiT\nmMaRH079wYyKljo4DOpMrYJCp5k6QwVLSjuMlhbQSBVL7QRjQ4BgU2ub1pVWkQ4/EpwsIWvICobF\nbEh27+63f9yzyc1y7r3n5jnn3HPvvl8zGe79nuc898n1+sm5z3POuebuAgCcuI5mDwAAWh1BCgCB\nCFIACESQAkAgghQAAhGkABCIIEVLM7OPmtkPK57/2syWNG9EmI4IUhSemb3dzHaY2Stm9rKZ/YeZ\nvTWurbuf7O4DKb/+9Wb2uJkdNrOvTdm2PNr2q2h8/2lmb0/z9VF8Xc0eAFCLmZ0i6V8kXSvpAUkz\nJf2upCM5DmNQ0hclvUvS7Jht75M0ED2/XtJ3JJ2e1+DQfByRouiWSXJ3/0cvO+zu29z9ybjGZjZh\nZr3R49lmdruZDURHiz80s1nRtguio9wDZvZTM7uk2gDcfZO7f1fSyzHbht19t5cvEeyUNCFpXwp/\nb7QQjkhRdM9KGjezr0v6tqQfufuBhPv+paQ3S7pQ0n5JKyRNmNmZKh/lfsjdHzazVZL+2cx+y91f\nqtGfVd1g9oqkOZJ+IenShONDm+CIFIXm7r+W9HZJLunvJQ2Z2XfN7LRa+5lZh6SPSfq0u+9z9wl3\n/293H5X0IUkPufvD0Ws8JulxSevqDafGOHskzVM57P/JzKqGLtoPQYrCc/dn3P1j7r5Y0m9LOkPS\nHXV2e6OkWZJ+HrPtLEnvi77WHzCzA5IuVv15zZrh6O6HJN2o8nTEW+r0hTZCkKKluPuzkjaqHKi1\nvCTpsKTfiNm2R9J97j6/4s/J7n5bvZdPMMROlf9/dShBW7QJghSFZma/aWZ/HM1ryswWS7pK0n/V\n2s/dJyTdK2mDmS00s04zu9DMuiX9g6T3mNmaqD7LzFZOvkbMGDqjRaouSZ1mNtPMOqNtq8zsnKjN\nKZI2SHrW3Z9L6z1A8RGkKLpfSzpf0o/M7KDKAfozSZ+JtruOP1KsfPwnkp6U9GOVV9z/QlKHu++V\ndLmkL0gaUvkI9TOq/v+HP1X5CPMGledXX5N0U7StR9L9kl5ReWHsVEmXndhfFa3KuLEzAIThiBQA\nAhGkABCIIAWAQAQpAARqu0tEzYzVMwCZcPfYizLaLkjL1ufwGtslvbMJ/VVr1+h4krSv1SZuW5Ja\n6PM0+5z631YZ91SN/u90Im2StM/itULGk7b1Vbfw1R4AAhGkABCIIG05S5o9gDaypNkDaDNLmj2A\npiFIW87SZg+gjfBepmv6vp8EKQAEIkgBIBBBCgCBCFIACESQAkAgghQAAhGkABCIIAWAQAQpAAQi\nSAEgEEEKAIEIUgAIRJACQCCCFAACEaQAEIggBYBABCkABCJIASAQQQoAgQhSAAhk7t7sMaTKzFy6\npNnDANB2fiB3t7gtXXkPJR/vzOE1tqf8OqH9Nbp/kva12sRtS1ILfZ5Fn6067nr1etsaaRPSPq19\n8+wzzg+qbuGrPQAEIkgBIBBBCgCBCFIACESQAkAgghQAAhGkABCIIAWAQAQpAAQiSAEgEEEKAIEI\nUgAIRJACQCCCFAACEaQAEIggBYBABCkABCJIASAQv9kEAInwm00Z4Deb+M2mIvRZq15vWyNtQtqn\ntW+efcbhN5sAIDMEKQAEIkgBIFBTg9TMBszsZ2b2EzPrj2pvMLNtZrbTzB41s56K9p83s11m9oyZ\nrWneyAHgmGYfkbqkle5+rruviGo3Strm7sskfT96LjNbLukDkpZLWivpK2bW7PEDQNODVJKmnk5w\nmaSN0eONkq6IHl8u6X53H3P3AUnPSVohAGiyZgepS3rMzB43s49HtQXuvj96vF/SgujxGZL2Vuy7\nV9KZ+QwTAKpr9nmkF7v7PjM7VdI2M3umcqO7e/kE+6ra62oCAC2pqUHq7vui//7SzDap/FV9v5md\n7u4vmtlCSUNR80FJiyt2XxTVYmyveLxE0tJ0Bw5gGtgtaSBRy6Z9tTezk8zs5OjxHElrJD0pabOk\nq6NmV0t6MHq8WdIHzazbzJZK6pPUH9/7Oyv+EKIATsRSHZ8l1TXziHSBpE1mNjmOb7r7o2b2uKQH\nzOwalf85eL8kuftTZvaApKcklSRd5+12owAALalpQeruuyWdE1P/laRVVfa5RdItGQ8NABrS7FV7\nAGh5BCkABOJ+pACQCPcjzQD3I+V+pEXos1a93rZG2oS0T2vfPPuMw/1IASAzBCkABCJIASAQQQoA\ngQhSAAhEkAJAIIIUAAIRpAAQiCAFgEAEKQAE4lp7AEiEa+0zwLX2XGtfhD5r1etta6RNSPu09s2z\nzzhcaw8AmSFIASAQQQoAgQhSAAhEkAJAIIIUAAIRpAAQiBPyASARTsjPACfkc0J+EfqsVa+3rZE2\nIe3T2jfPPuNwQj4AZIYgBYBABCkABCJIASAQQQoAgTj9CQAS4fSnDHD6E6c/FaHPWvV62xppE9I+\nrX3z7DMOpz8BQGYIUgAIRJACQCAWmwAgERabMsBiE4tNReizVr3etkbahLRPa988+4zDYhMAZIYg\nBYBABCkABGKxCQASYbEpAyw2sdhUhD5r1etta6RNSPu09s2zzzgsNgFAZghSAAjEHCkAJMIcaQaY\nI2WOtAh91qrX29ZIm5D2ae2bZ59xmCMFgMwQpAAQiDlSAEiEOdIMMEfKHGkR+qxVr7etkTYh7dPa\nN88+4zBHCgCZIUgBIBBzpACQCHOkGWCONIs50rvu+h0tW/YW9fbOV8+6S1UaH1VprEMzZ3ZpeOu/\nVqldotJ4h0pj41HtIfWsW6fSuEc1aXjrv0f7ujbfcIuu/PJNmjWrU0eOjKs0PqrNN9ynK798TY3a\nqErjHdp8w71RrUNHjkxE/d2rK798tWbNmhnte+w1xidcu3a9rPHxCV2ssxt6L5gjbWafcarPkbZp\nkKJVLVu2WCtXLik/OWWmpJkaHZtQ94wOzXvTvCq1kySporYgaqeYfaV1687XvL8rP545s0vSTK1b\n15egpik1VdTmTGl37DXO7p2vnbtezuLtQkEQpCiUQ4cOS5JeeeWwenbsUH//0xoe7tbq1WfXqQ1q\nePhwVDsYU6to98TPtXrHjnKtZ1a5vydeTFAbrKgdVE/P3Cm1ynbl1+jvH9SaNfdpWEdUPnJCO2rT\nIM3rA5v264T21+j+SdrXahO3LUmt+vOhoRUaGjqg0dGS5qx+h/pGDqtUGtfovDk6uGlLA7VLq7bb\nevMGjb5thWa7a8ysonZVyrUV6hs5rEcOvqBxTehindvQe5Hseb16vW2NtAlpn9a+efaZXJsGKXOk\n6bTPf450yZKFOu20+eUnM7o0v2fu0a/nixadmkpt1arz1H3nsY9+uXZhyrXya3T3zNXZvWdVfLVn\njrQd50g5/QmF0tu7UJI0Pj4hSSqVxjV5Zknt2kTi2tjYuCRp8oSVUmk8YS1u39q1UmlCA88Ph70p\nKDyCFIWyZ8+QJKmzs/zR7OrqVIdZglpH4lp3d/noMdqkrq5OdXd3Jqh1NFDrOvp4yVnzwt4UFB7n\nkaJQtmz5ktatu6BiMedpDQ+PaPXqt7Zkrb//aa1Z81kNDx9s9luLYJxHmgHmSLOYIx0aOqChoRGN\njo5rzuoL1DdySKVSh0bnzdLBTY9Uqa1Q30hJpdJEVPue5qxepb6R0ajWrYObtkX7jmrrzbdp9G0X\naLZPaMw61DdySFtv/madmtQ3UtLWmzdGtXGNWWfU38Zo8UrRvsdeo29kVI8cfDlabOI8UuZIgRyU\nF5vmaNGieZoxo0Pze+Zq3rxZ0YJRtVqX5vfMqqidGrWbrHVV7DurvBA0o0Mzu7uO9rdqVW+dWle0\n72RtRkV/veqe0VWx77HXmN8zS2f3zm/224qMEaQoFBab0IoIUhQKi01oRSw2oVBYbEJxsdiUARab\nWGxisalxLDYBmWOxCa2IIEWhsNiEVtSmX+25aUl67fO9acmePe/Vm9604LiFJR8rSVKdWkfiWvxi\n00SCWkcDtcrFppO0c9feht+LZM/r1etta6RNSPu09s2zz+TadLFpfQ6vxBxpFnOkW7a8K1psqrid\nXeUt86rWptxGb/IWd5W30Tvabl+0OJT9axy9jd7wkYbfC+ZIm9lnnPXTbbEJrerYYlNJc1avjFlY\niqu9I+Fi08oCLDahHTFHikJhsQmtqE2PSJkjTa99vnOkvb2fkFReHOpU5cJSR2q1ysUhs8naaynX\nJhebxjXw/F5JYw2/F8me16vX29ZIm5D2ae2bZ5/JtWmQch5pOu3znyPds2coZrGpvJiTVi1+sWl2\nyrWuo4+XnLWIGzunsm+efcapfh5pywWpma2VdIekTkl3u/utTR4SUvTqqyOSMv7Npp8+9/rfZ/pp\nzG82va42WFGrWGyKbfccv9k0jbRUkJpZp6Q7Ja2SNCjpx2a22d2fbu7IkBYWm9CK6p7+ZGafknSf\nux/IZ0g1x3KhpD9z97XR8xslyd2/VNHGuda+dW3ffodWrjyn/OSiiyRJo2Mldc/oknbsSKU29OAW\nnXbFu4973SxrL700rJ279kY/fofWFXat/QKVj/z+V9K9kh7x5p18eqakFyqe75V0/uubMUeaTvv8\n50grr2w6tmBkCWoTiWvxi02dCWoTDdQqr2w6JOmU6G/IHGk7zpHWPf3J3W+StEzlEP2opF1mdouZ\nNeN7SsIA317xZ3eGw0HauI0eimO3js+S6hLNkbr7hJm9KGm/pHFJ8yV9x8wec/fPhg22IYOSFlc8\nX6zyUekUefzrhCyw2ITiWBr9mVT9iDTJHOmnJX1E0suS7pa0yd3HzKxD0i53z+3I1My6JD0r6fck\n/UJSv6SrKhebmCNtbV/72g1at+4CjY6WtODKd+vgyGGVSuOaN2+OhjZtSaW29eYN+oO/+oLcXR1m\nmdZGRg5r584XND4+wRxpywubI32DpPe6+/OVxego9T1pDC8pdy+Z2fWSHlH59Kd74lfsmSNNp33+\nc6TlK5uiK4FmdGl+z1yNjk0cvWIpjdqqVeep+85jH/1y7cKUa+XX6O6Zq7N7z+I80lT2zbPPOAHn\nkbr7n9XY9tQJjuiEuftWSVvzfl3kI36xSQlqx98yr1YtfrEpSS1u39o1bqM3PXCtPQqFxSa0opY6\nIR/tj8UmtKI2DVJuWpJe+3xvWjI0tEJDQweiK5veob5owWh03hwd3LSlgdqlVdttvXlDdCWSa8ys\nonZVyrUV6hs5rEcOvhBd2XRuQ+9Fsuf16vW2NdImpH1a++bZZ3JtGqQsNqXTnsUmFpvSaJ/Wvnn2\nGSfghHwgT/xmE1oRQYpCYbEJrahNf7OJE/Jb1ZYtX4p+s2lyMedpDQ+PRL+x1Hq1/v6ntWbNZzU8\nfLDZby2ChZ2Q34KYI02nff5zpMduozeuOasviLllXlxtRcLb6F1QgNvoMUfKHCmQMX6zCa2IIEWh\nsNiEVtSmX+05jzS99vmeR7pnz3tjfrOpJGnqbzFNrXUkrsUvNk0kqHU0UKtcbDpJO3dN3qSM80jb\n8TzSNl1sWp/DKzFHmsUc6ZYt74oWmyquMKq8iqlqbcqVTZNXHVVe2XS03b5ocSj71zh6ZdPwkYbf\nC+ZIm9lnnPXTbbEJrYrfbEIrYo4UhcJiE1pRmx6RMkeaXvt850h7ez8haeot81xSR2q1+FvmvZZy\nbXKxaVwDz++VNNbwe5Hseb16vW2NtAlpn9a+efaZXJsGKeeRptM+/znSPXuGYhabyos5adXiF5tm\np1zrOvp4yVmLuNY+lX3z7DNO9fNI2zRI0aq4jR5aEUGKQmGxCa2oTU9/uqTZw8AJ2r79Dq1ceU75\nyUUXSZJGx0rqntEl7diRSm3owS067Yp3H/e6WdZeemlYO3ft5cfvWh7X2meAOdIs5kjjf7PJEtQm\nEtfiF5s6E9QmGqhVXtl0SNIp0d+QOdJ2nCPl9CcUCrfRQytq0yNStCoWm9CKmCNFocybN1d33fUZ\nfe5zf6vbbvuErr32dklq2dq1197OvUjbRvU50jYN0vU5vBJzpFnMkSabF0y7z1Ydd716vW2NtAlp\nn9a+efYZp/q19syRAkAgghQAAhGkABCoTedIWWwCkDZOyM8Ai00sNhWhz1r1etsaaRPSPq198+wz\nDifkA0BmCFIACMQcKQAkwhxpBpgjZY60CH3Wqtfb1kibkPZp7Ztnn3GYIwWAzBCkABCIOVIASIQ5\n0gwwR8ocaRH6rFWvt62RNiHt09o3zz7jMEcKAJkhSAEgEHOkAJAIc6QZYI6UOdIi9FmrXm9bI21C\n2qe1b559xmGOFAAyQ5ACQCCCFAACEaQAEIhVewBIhFX7DLBqz6p9EfqsVa+3rZE2Ie3T2jfPPuOw\nag8AmSFIASAQQQoAgQhSAAhEkAJAIE5/AoBEOP0pA5z+xOlPReizVr3etkbahLRPa988+4zD6U8A\nkBmCFAACEaQAEIggBYBABCkABCJIASAQQQoAgTghHwAS4YT8DHBCPifkF6HPWvV62xppE9I+rX3z\n7DMOJ+QDQGYIUgAIRJACQCCCFAACEaQAEIggBYBABCkABCJIASAQQQoAgQhSAAjEtfYAkEjBrrU3\ns/WS/kjSL6PSF9x9a7Tt85L+UNK4pE+5+6NR/TxJX5c0S9JD7v7p6q/AtfbptOda++L3Wateb1sj\nbULap7Vvnn3GKd619i5pg7ufG/2ZDNHlkj4gabmktZK+YmaT/wJ8VdI17t4nqc/M1jZj4AAwVTPn\nSOMOkS+XdL+7j7n7gKTnJJ1vZgslnezu/VG7b0i6Ip9hAkBtzQzST5rZE2Z2j5n1RLUzJO2taLNX\n0pkx9cGoDgBNl1mQmtk2M3sy5s9lKn9NXyrpHEn7JN2e1TgAIGuZLTa5++ok7czsbknfi54OSlpc\nsXmRykeig9Hjyvpg9V63VzxeonJmA0AjdksaSNSyWav2C919X/T0SklPRo83S/qWmW1Q+at7n6R+\nd3cze9XMzpfUL+nDkv6m+ivksYIHoL0t1fEHYdVX7Zv1UyO3mtk5Kq/e75Z0rSS5+1Nm9oCkpySV\nJF3nx050vU7l059mq3z608O5jxoAYjQlSN39IzW23SLplpj6/0h6S5bjAoATwSWiABCIIAWAQAQp\nAAQiSAEgEEEKAIG4jR4AJFKw2+hlj9vopdOe2+gVv89a9XrbGmkT0j6tffPsM07xbqMHAG2DIAWA\nQAQpAAQiSAEgEEEKAIEIUgAIRJACQCCCFAACEaQAEIggBYBABCkABCJIASAQQQoAgQhSAAhEkAJA\nIIIUAAIRpAAQiCAFgEAEKQAE4sfvACARfvwuA/z4HT9+V4Q+a9XrbWukTUj7tPbNs884/PgdAGSG\nIAWAQAQpAAQiSAEgEEEKAIEIUgAIRJACQCCCFAACEaQAEIggBYBABCkABCJIASAQQQoAgQhSAAhE\nkAJAIIIUAAIRpAAQiCAFgEAEKQAEIkgBIBBB2nJ2N3sAbYT3Ml3T9/0kSFvOQLMH0EYGmj2ANjPQ\n7AE0DUEKAIEIUgAIZO7e7DGkysza6y8EoDDc3eLqbRekAJA3vtoDQCCCFAACEaQFZWbrzWyvmf0k\n+vP7Fds+b2a7zOwZM1tTUT/PzJ6Mtv11c0beGsxsbfT+7TKzG5o9nlZgZgNm9rPo89gf1d5gZtvM\nbKeZPWpmPRXtYz+n7YggLS6XtMHdz43+bJUkM1su6QOSlktaK+krZjY5Af5VSde4e5+kPjNb24yB\nF52ZdUq6U+X3b7mkq8zszc0dVUtwSSujz+OKqHajpG3uvkzS96Pn1T6nbZs3bfsXaxNxK4SXS7rf\n3cfcfUDSc5LON7OFkk529/6o3TckXZHPMFvOCknPufuAu49J+rbK7yvqm/qZvEzSxujxRh37zMV9\nTleoTRGkxfZJM3vCzO6p+Mp0hqS9FW32Sjozpj4Y1fF6Z0p6oeL55HuI2lzSY2b2uJl9PKotcPf9\n0eP9khZEj6t9TttSV7MHMJ2Z2TZJp8dsuknlr+l/Hj3/oqTbJV2T09DaHef8nZiL3X2fmZ0qaZuZ\nPVO50d29znncbfu+E6RN5O6rk7Qzs7slfS96OihpccXmRSr/az8YPa6sD6YwzHY09T1crOOPnhDD\n3fdF//2lmW1S+av6fjM73d1fjKaXhqLmcZ/Ttv088tW+oKIP5aQrJT0ZPd4s6YNm1m1mSyX1Sep3\n9xclvWpm50eLTx+W9GCug24dj6u8GLfEzLpVXhTZ3OQxFZqZnWRmJ0eP50hao/JncrOkq6NmV+vY\nZy72c5rvqPPDEWlx3Wpm56j8dWi3pGslyd2fMrMHJD0lqSTpOj92edp1kr4uabakh9z94dxH3QLc\nvWRm10t6RFKnpHvc/ekmD6voFkjaFJ0g0iXpm+7+qJk9LukBM7tG5ds/vV+q+zltO1wiCgCB+GoP\nAIEIUgAIRJACQCCCFAACEaQAEIggBYBABCkABCJIASAQQYppx8zeFt1Va6aZzTGz/4vunwmcEK5s\nwrRkZl+UNEvly2lfcPdbmzwktDCCFNOSmc1Q+eYlr0m6sJ2vA0f2+GqP6eqNkuZImqvyUSlwwjgi\nxbRkZpslfUtSr6SF7v7JJg8JLYzb6GHaMbOPSDri7t+OfpBth5mtdPd/a/LQ0KI4IgWAQMyRAkAg\nghQAAhGkABCIIAWAQAQpAAQiSAEgEEEKAIEIUgAI9P/aWJGg1kIsuwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x155f4a20>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,1, figsize = (5,5))\n",
    "mesh.plotSlice(sigma, grid=True, ax = ax)\n",
    "ax.plot(xyz_rxP[:,0],xyz_rxP[:,1], 'w.')\n",
    "ax.plot(xyz_rxN[:,0],xyz_rxN[:,1], 'r.', ms = 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1323\n"
     ]
    }
   ],
   "source": [
    "rx = DC.RxDipole(xyz_rxP, xyz_rxN)\n",
    "tx = DC.SrcDipole([rx], [-200, 0, -12.5],[+200, 0, -12.5])\n",
    "print xyz_rxP.size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "survey = DC.SurveyDC([tx])\n",
    "problem = DC.ProblemDC_CC(mesh)\n",
    "problem.pair(survey)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    from pymatsolver import MumpsSolver\n",
    "    problem.Solver = MumpsSolver\n",
    "except Exception, e:\n",
    "    problem.Solver = SolverLU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "ImportError",
     "evalue": "No module named pymatsolver",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mImportError\u001b[0m                               Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-11-d6799536a06f>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[1;32mfrom\u001b[0m \u001b[0mpymatsolver\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mMumpsSolver\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mImportError\u001b[0m: No module named pymatsolver"
     ]
    }
   ],
   "source": [
    "from pymatsolver import MumpsSolver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "problem.Solver = SolverLU\n",
    "\n",
    "data = survey.dpred(sigmahomo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "IndentationError",
     "evalue": "expected an indented block (<ipython-input-13-ff092c0e869a>, line 3)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-13-ff092c0e869a>\"\u001b[1;36m, line \u001b[1;32m3\u001b[0m\n\u001b[1;33m    \u001b[0m\n\u001b[1;37m    ^\u001b[0m\n\u001b[1;31mIndentationError\u001b[0m\u001b[1;31m:\u001b[0m expected an indented block\n"
     ]
    }
   ],
   "source": [
    "# Plot pseudo section\n",
    "for ii in range(data):\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "u1 = problem.fields(sigma)\n",
    "u2 = problem.fields(sigmahomo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Msig1 = Utils.sdiag(1./(mesh.aveF2CC.T*(1./sigma)))\n",
    "Msig2 = Utils.sdiag(1./(mesh.aveF2CC.T*(1./sigmahomo)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "j1 = Msig1*mesh.cellGrad*u1[tx, 'phi_sol']\n",
    "j2 = Msig2*mesh.cellGrad*u2[tx, 'phi_sol']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# us = u1-u2\n",
    "# js = j1-j2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'mesh' 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-2-cb76a57fca1d>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mmesh\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplotSlice\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmesh\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0maveF2CCV\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mj1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvType\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'CCv'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnormal\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'Y'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mview\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'vec'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mstreamOpts\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m{\u001b[0m\u001b[1;34m\"density\"\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m\"color\"\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;34m'w'\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      2\u001b[0m \u001b[1;31m#xlim(-300, 300)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[1;31m#ylim(-300, 0)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'mesh' is not defined"
     ]
    }
   ],
   "source": [
    "mesh.plotSlice(mesh.aveF2CCV*j1, vType='CCv', normal='Y', view='vec', streamOpts={\"density\":3, \"color\":'w'})\n",
    "#xlim(-300, 300)\n",
    "#ylim(-300, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'js' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-23-575f23801c4a>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mmesh\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplotSlice\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmesh\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0maveF2CCV\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mjs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvType\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'CCv'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnormal\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'Y'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mview\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'vec'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstreamOpts\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m{\u001b[0m\u001b[0;34m\"density\"\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"color\"\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m'w'\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mxlim\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m300\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m300\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mylim\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m300\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'js' is not defined"
     ]
    }
   ],
   "source": [
    "mesh.plotSlice(mesh.aveF2CCV*js, vType='CCv', normal='Y', view='vec', streamOpts={\"density\":3, \"color\":'w'})\n",
    "xlim(-300, 300)\n",
    "ylim(-300, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = np.random.randn(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print (a.reshape([1,-1])).repeat(3, axis = 0)\n",
    "print (a.reshape([1,-1])).repeat(3, axis = 0).sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def DChalf(txlocP, txlocN, rxloc, sigma, I=1.):\n",
    "    rp = (txlocP.reshape([1,-1])).repeat(rxloc.shape[0], axis = 0)\n",
    "    rn = (txlocN.reshape([1,-1])).repeat(rxloc.shape[0], axis = 0)\n",
    "    rP = np.sqrt(((rxloc-rp)**2).sum(axis=1))\n",
    "    rN = np.sqrt(((rxloc-rn)**2).sum(axis=1))\n",
    "    return I/(sigma*2.*np.pi)*(1/rP-1/rN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data_analP = DChalf(np.r_[-200, 0, 0.],np.r_[+200, 0, 0.], xyz_rxP, sighalf)\n",
    "data_analN = DChalf(np.r_[-200, 0, 0.],np.r_[+200, 0, 0.], xyz_rxN, sighalf)\n",
    "data_anal = data_analP-data_analN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Data_anal = data_anal.reshape((21, 21), order = 'F')\n",
    "Data = data.reshape((21, 21), order = 'F')\n",
    "X = xyz_rxM[:,0].reshape((21, 21), order = 'F')\n",
    "Y = xyz_rxM[:,1].reshape((21, 21), order = 'F')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,2, figsize = (12, 5))\n",
    "vmin = np.r_[data, data_anal].min()\n",
    "vmax = np.r_[data, data_anal].max()\n",
    "dat0 = ax[0].contourf(X, Y, Data, 60, vmin = vmin, vmax = vmax)\n",
    "dat1 = ax[1].contourf(X, Y, Data_anal, 60, vmin = vmin, vmax = vmax)\n",
    "cb0 = plt.colorbar(dat1, orientation = 'horizontal', ax = ax[0])\n",
    "cb1 = plt.colorbar(dat1, orientation = 'horizontal', ax = ax[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
