examples/Creating new fields.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:05280da19004e75d9e31fd74082606bdf82b0045db4c82eb190399c9cb81a04e"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Adding new analytic field models"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Currently RAPT has a limited selection of field models. However, it is easy to add new field models to it. In this notebook we describe how to add analytic fields (i.e., fields that have an explicit formula in terms of t, x, y, and z)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "All new fields are subclasses of the `_Field` class. By subclassing, the new field class inherits all the auxiliary field-related functions, so we don't have to write them  again. All we need to provide is the new `E` and `B`. This is called _overriding_ in OOP parlance."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Overriding the static magnetic field: Double-dipole model"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This model is already built into the RAPT package, but it provides a good example for creating new fields. The double-dipole model consists of two identical dipoles, both moments in the negative-z direction, separated by 20$R_e$. This arrangement generates a planar magnetopause between them, and mimics the compression of the dayside magnetosphere by the solar wind."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import pylab as pl"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import rapt\n",
      "import numpy as np"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "class DoubleDipole(rapt.fields._Field):\n",
      "    def __init__(self, B0=rapt.B0, distance=20*rapt.Re, imagestrength=1):\n",
      "        rapt.fields._Field.__init__(self)\n",
      "        self.gradientstepsize = rapt.Re/1000\n",
      "        self._dd = distance  # distance between two dipoles\n",
      "        assert imagestrength >= 1\n",
      "        self._k = imagestrength   # >=1. Relative strength of the image dipole\n",
      "        self._coeff = -B0*rapt.Re**3 \n",
      "        \n",
      "    def B(self, tpos):\n",
      "        t,x,y,z = tpos\n",
      "        B1 =  np.array([3*x*z, 3*y*z, (2*z*z -x*x- y*y)]) / pow(x*x+y*y+z*z, 5.0/2.0)\n",
      "        x -= self._dd\n",
      "        B2 = self._k * np.array([3*x*z, 3*y*z, (2*z*z -x*x- y*y)]) / pow(x*x+y*y+z*z, 5.0/2.0)\n",
      "        return self._coeff*(B1+B2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Extending with an electric field: Electrically charged magnet"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As an entirely new model, consider a magnetic dipole coincident with an electric charge. To construct a new class, we override both `E` and `B`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "class ChargedDipole(rapt.fields._Field):\n",
      "    def __init__(self, B0=1, Q=1):\n",
      "        rapt.fields._Field.__init__(self)\n",
      "        self.B0=B0  # field strength (T) at unit distance on the x-y plane\n",
      "        self.Q=Q    # charge (C)\n",
      "        self._k = 8.9875517873681764e9 # Coulomb constant\n",
      "        self.static = False  # needs to be set to false if E and/or dB/dt is nonzero.\n",
      "    \n",
      "    def B(self, tpos):\n",
      "        t,x,y,z = tpos\n",
      "        return self.B0 * np.array([3*x*z, 3*y*z, (2*z*z -x*x- y*y)]) / pow(x*x+y*y+z*z, 5.0/2.0)\n",
      "    \n",
      "    def E(self, tpos):\n",
      "        t,x,y,z = tpos\n",
      "        return  self._k*self.Q * np.array([x,y,z]) / pow(x*x+y*y+z*z, 3.0/2.0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 152
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Try it out with a small charge."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f = ChargedDipole(Q=1e-6)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 153
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Follow a proton's trajectory"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "p = rapt.Particle([5,0,0], [0,1,0], t0=0, mass=rapt.m_pr, charge=rapt.e, field=f)\n",
      "p.setke(1)  # set kinetic energy to 1 eV"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 155
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Check the cyclotron period and cyclotron radius to gain an idea about the scale of the motion."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "p.cycper()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 160,
       "text": [
        "8.5710888990731629e-06"
       ]
      }
     ],
     "prompt_number": 160
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "p.cycrad()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 161,
       "text": [
        "0.098787549981589062"
       ]
      }
     ],
     "prompt_number": 161
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Follow the proton for some time."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "p.advance(4e-4)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 157
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Display the trajectory on the x-y plane."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pl.plot(p.getx(), p.gety())\n",
      "pl.xlabel(\"x\")\n",
      "pl.ylabel(\"y\");"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 159,
       "text": [
        "<matplotlib.text.Text at 0x7f5f6f52a630>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEPCAYAAABP1MOPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xuc1nP+//HHK1NUKlI5lEpClFPZCMtFTjnV7g9FWDnH\nhnW27Dbsrg2Lxa79bigrRcTqoGzIZYksSTqanNqcKkUHlZrm/fvjNbUdZ6bpuq73dXjeb7e5mcM1\n83lepvm8rvfZQgiIiEjhqhE7gIiIxKVCICJS4FQIREQKnAqBiEiBUyEQESlwKgQiIgUueiEwswZm\n9qyZzTCzaWZ2aOxMIiKFpCh2AOABYHQI4UwzKwLqxA4kIlJILOaCMjOrD0wKIewZLYSISIGL3TW0\nB/CtmQ00s/fNrL+Z1Y6cSUSkoMQuBEVAe+CvIYT2wDLg5riRREQKS+wxgi+AOSGE98o/HgbctOGD\nzEwbIomIVEMIwSp7TNQWQQhhLjDHzPYu/1RnYPpmHpu3b3379o2eQc9Pz03PL//eqip2iwDgKmCw\nmdUEPgV6Rc4jIlJQoheCEMJk4Cexc4iIFKrYg8UCJBKJ2BHSKp+fXz4/N9DzKxRR1xFUlZmFXMgp\nIpJNzIyQ7YPFIiISnwqBiEiBUyEQESlwKgQiIgVOhUBEpMCpEIiIFDgVAhGRAqdCICJS4FQIREQK\nnAqBiEiBUyEQESlwKgQiIgVOhUBEpMCpEIiIFLjoB9OY2efAIqAMWBVC6Bg3kYhIYYleCPACkAgh\nfBc7iMiGli+HEGC77aBGBe3nefNgwQIoK4PWrWHbbTOXUWRrZUMhMNRFJRkSAnz1FXz3HWy/Pey+\nO2yzzfqPmTgRHnoIkkn4+mv/elERdOwIv/gF9OgBNWvC4sXw4IPw5JMwdy7ssov//NmzoUMHuOIK\nf+yaArJkCYwdC9OmwbffQuPG0KYNdOniWURiyYYbcAD+ZWbvmtklscNIfpoxAy6/HBo1goMPhu7d\n4ac/9Ztxz57w4YewahVcdRV07Qr77w//+hesWAHLlsFnn8GVV8Jjj8Ghh8KQIf6YmTO9ECxY4NeY\nORO+/x6uvx4eeACOPx7mzIHbb4fmzaF/f1i5EvbYw1sbAwdCs2Zw3XXwww/rZ16yBKZMgalT/eeL\npEs2tAiOCCF8bWaNgZfNbEYI4c0NH1RcXLz2/UQiobNGZa2ZM/3V+2efgRnsthsccwy0a+cf//nP\n8Ic/QJ8+8MEH3gpY4+uvYfBgv2HXru3fM3Uq7LDD+tfYaSf42c+gWzf429/gj3/0G323bhvn2XZb\n//xpp8Ett3gB6N7dr92ixcaP/+YbLwSdOnmL4aWX/BpTp/rjzeDLL/3nnHOOF6s6dVL6v1DyRDKZ\nJJlMbvH3ZdWZxWbWF1gSQrhvg8/rzGLZyH/+A1df7V0xXbrAnnv6TXP2bBgzBpo2hfbt4ZVX/Oba\nsuXmf9a8eTB+vLcGKhoLqI7vv9+4sGwoBPjtb+HOO70g3HYbHHecd0kBrF4NEyZ48fnPf2D0aNhv\nv/99/7Jl3qKoV8/HM0Sg6mcWRy0EZlYHqBFCWGpmdYGxwO0hhLEbPE6FoACt6c+fNw/q1vVXxGtu\nco88Ar/5Ddx7r/fDb9jPv3o1PP203zD79Vu/FZDN5szxriKr4E930CC4+Wb/bzIJQ4fCf//rrYSl\nS2HvveGUU+BXv4Kdd85YdMlCuVII9gD+iY8TFAGDQwj9NvE4FYIC8tVXcNdd8NxzUFrqN7MffvAu\nlM6dvfvmH//wm2Dr1rHTxjFgABQXe3fVL34BBxzgrYfVq+G993wM48kn4Z574MILY6eVWHKiEFSV\nCkHhGDzYu3suvNDf2rT539cWLYJhw+D55+F3v/NuH9m8mTPh1FN94Pryy72ozpwJCxdCrVpeUDVb\nKb+pEEhWW7zY3+rWhR139M8NHeo3rTFj/CYlW2/WLDj8cB9vGDPGW1dNmviMpRkzfAbV9dfD6aen\nfmxE4lMhkKwzfz7cdx+MGAGff+4DqEuWQIMG/sr1ued8UFev9FPrtdd8xlL37j6jao0ff4RRo7x1\n1bIlPPEE1K8fLaakgQqBZJXBg33a49lnwwUX+CvRbbbxAeGSEnj2WdhrL79ZSWatXAmXXeaL3EaM\nqHigWnKLCoFkjQEDfEHVqFG+CEuyz8qVPm21a1dvob37rk97rVvXf2drFtlJblEhkIxbvty3Z5g7\n12ew7Lefv7rs1AneeGP9gV/JPmtWX++7r6+ebtTIx3EmTYKnnvLfZ//+vipacoMKgWTMtGm+EOqF\nF/xm0by59z9PngxffOHz+G+6KXZK2RqrVvn+S/36+e/58MNjJ5KqUCGQjHjkEbj1Vrj2Wujd27sV\n1gjBt0Zo2lT9zvlizBif1jt5ss8+kuymQiBpN3iwr+596SVfzSqF4corvYVw6qk+PXXFCm8FHnaY\nD/hL9lAhkLRatAhatYJx4+DAA2OnkUz68ksvAk2a+HjCdtv5Fhfjxvm/ibvvhiOPjJ1SQIVAUuTH\nH31q54gRvhvmihW+dXNRke+MOWRI7ISSLUpLfeX3tdf6WQy33RY7kagQyFb717+8P7hdOzj3XH/l\nX7eu7/kzfbpPKVQ/sWxo3jw46igvCJdeGjtNYVMhkK3y6KO+qdmgQb63v8iWmDrVz3j49FM/50Hi\nUCGQaps4EU4+Gd58U4N/Un1du/r5CLvu6ttjN2nix30ed5zOdM6UqhYCbTMlG+nb11sDKgKyNe67\nz89WaNzYVyWH4NuLN23qZz2XlsZOKGuoRVDAysrg3//2bZ0/+shXBu+wA7z6qm8Qp+MQJR2mT/cp\nqNtt5//21HWUPuoakgq9/74PBJeV+UZwBxzge9N/843PCe/UKXZCyWelpXDeeb7QUDPP0ienCoGZ\n1QDeA74IIZy+ia+rEKTQiy9Cr17edO/ZU6t+JY5ly3xG2hNPaN1BuuTaGMHVwPTYIQpBSYlvAz1y\npE8JVRGQWOrUgUsu8bOllyzxN73eiyN6ITCzZsDJwKOxsxSC227zE6kOPTR2EhE480zfr2rXXf2t\nQQPo1s1brSoKmRO9a8jMngX+ADQArlPXUHqsWOGvuFq3hjlzdBKVZI+ysv8dkzlvnu9dddddXhiG\nDoWddoqbL5dVtWuoKBNhNsfMTgHmhhA+MLMEsNnAxcXFa99PJBIkEol0x8tpZWX+RzR4sK8HWLbM\nB+guvlhFQLLLumclN2kC55/vY1e33QZHHAHjx6sYVFUymSSZTG7x90VtEZjZncC5QClQG6gHPB9C\nOH+Dx6lFsAVKSrz/v6gI+vTxBTyNGnlT20zjApI7rrnGDzp66qnYSXJTTs0aAjCzo1HX0Fb76CM4\n9lj49a994y/d9CWXLV8Ou+wCH3/sC9Nky+TarCFJgdWrvVl9yy2+YEdFQHJd7dqQSMDYsbGT5Les\nKQQhhNc31RqQqhszxscGrrgidhKR1OnZ06c8160Le+7pC9HGjdOsolTKmkIgW+/FF32VcA39ViWP\nnHUWrFzpYwUjR/pJaL17+9jX3Lmx0+WHrBkjqIjGCDZt/nwYPhzeeMOn3SWT8MorPtNCJJ+tXg13\n3AGPPw5vveUb2cnGcm6wuCIqBOsrLYU//MF3cDzhBB8c3m033zBORUAKye9/7+MHr7+uMbFNUSHI\nU8uXe1P5xx+hf39o2TJ2IpF4ysp8w8QHHoDOnWOnyT6aNZSnbrwRatb08QAVASl0NWr4TLnhw2Mn\nyW0qBDlkxgw/SP6xx7wYiAi0aQOvvebnaHzxhWYTVYcKQQ556ilfMbzjjrGTiGSPI4+EAw+E3/0O\n2reHffeFv/5VJ6BtCY0R5JDTToOLLvLdGUVkYyH43lq33+7jac88U9gzijRYnAdWrfJ/yKNGwdSp\nMGuWH/PXqlXsZCLZrazMp5f+85++ad3228dOFIcKQY4bP95XULZs6d1BHTr4K5tGjWInE8kNIfiK\n5J128tP4CpEKQQ574QW49FIYMABOPTV2GpHc9eWXfhzmnDmF2SrQ9NEc9cknXgRGj1YRENlaTZvC\nPvvApEmxk2Q3FYIs06+fbxp3yCGxk4jkh112gZkzfdxANk2FIIuE4IPDl14aO4lI/ujcGe6808cK\nzj7b9+aS9akQZJEFC3yl5G67xU4ikj/69IHPPvMFmUce6ZMvLrnEz/EWp0KQBUpLfdOsBx/Uxlki\n6bLLLn5g05QpsHChtw606MxFLQRmtq2ZvWNmk8xsipn1jZkn00KA//s/2GMPuPZa+O47GDQodiqR\n/Fa/vq/SX7gQHn44dprsEH36qJnVCSEsM7NtgPHAVSGE/2zwmLybPrpqFfToAbNnw9//7usERCRz\nJk+GE0/0/YmKimKnSY+cmT4aQlhW/u62QBGQX3f8zejTx09deustFQGRGA480AeQP/wwdpL4ohcC\nM6thZpOAb4CXQwjvxs6Ubu+959tIDxkCtWrFTiNSuFq2hK++ip0ivugNohBCGXCwmdUHXjCz/UII\n0zd8XHFx8dr3E4kEiUQiYxlT7eGH4ZproF692ElECtuSJX7E60EHQbNmsdNsvWQySTKZ3OLviz5G\nsC4z+w3wQwjhvg0+n1djBAce6GcKaNGYSFxjxvjf4muv+d/ln/7kW1nni5wYIzCzRmbWoPz92sDx\nwMyYmTJh2TJo0CB2ChHp0gWGDYOvv4aePf3jRx6JnSrzYncN7Qr8w8xq4EVpaAhhdORMafPNNzBy\nJHz7reYvi2STWrX8rI9EAo491g9/OuOM2KkyJ6u6hjYn17uG5s+Hq66Cl17y6Wo/+xmceaavIhaR\n7PL2214EZs7M/XG8qnYNxW4R5L05c+CnP/Ub/1dfQe3asROJSEU6dYKDD/bW+znnxE6TGXpNmkal\npXDKKb5m4J57VAREcsXxxxfW5nQqBGk0cCA0buzbR4hI7qhXz888LhQqBGk0ZAj86lfaSE4k13z2\nmf931aq4OTJFhSCNJk+Gjh1jpxCRLXXccfDxx77I7M47fcp3PlMhSCPNChLJTUcfDW++6auOJ06E\nQw/1zenylW5VafL99z44PGNG7CQiUl377usLzs491086y9eWgQpBis2Z4wdetGgBP/kJtGsXO5GI\nbA0zuOkm3xKmb56emKIFZSk0YQJ07Qq9e/tMofr1YycSkVSZPdu3jP/6a6hZM3aaqtGCsgybPRtO\nPx0efxxOPjl2GhFJtRYt/DzxDz/MvzNE1DWUItdd5wvHVARE8lfDhrBoUewUqaeuoRRYuNAPuJg7\nV6uHRfJVCLD77vD667DnnrHTVE1ObEOdL955x9cLqAiI5K/vvoMvv/Q9iBYvjp0mtVQIUmDFCqhb\nN3YKEUmnhg29NTBhAuy1l08rzRcaLE6BPfeEkpLYKUQk3Y46yt8mToSf/9y3mO/dO3aqrRf7hLJm\nZjbOzKaZ2RQzuypmnuooK4OPPvK3FStipxGRTOjQwY+37NvXi0Kui901VApcG0JoC3QCrjSzNpEz\nVdn48b5g7O674Z//hO22i51IRDKlVSv47W+hX7/YSbZeVs0aMrMXgIdCCK9u8PmsmzU0dqyfcfrI\nI76ITDuMihSeRYtg1119IHnbbWOn2VjOzRoys5bAQcA7cZNUbsECuOACHyzq1k1FQKRQNWgATZrk\n/oZ0WVEIzGx7YBhwdQhhaew8lfnrX+HUU32HQhEpbMuW5f7U8eizhsysCC8Cg0IIwzf3uOLi4rXv\nJxIJEolE2rNtznPPwcMPR7u8iGSJadOgVi3vHsoGyWSSZDK5xd8XfYzAzJ4Avg0hbPZAx2wbI9hx\nR/jkE59XLCKF68kn4cYb4W9/873Gsq2bOCfGCMzsCKAncKyZTTKz983spJiZqqJWrfzdl1xEqu7c\nc2HAALjlFujRA5YsiZ2oeqK3CKoi21oEp58OZ5wB558fO4mIZIPly+HSS32B2ahRUBS9093lRIsg\nFy1fDnXqwOjRsZOISLaoXRsGDoQff/TJJLlGhWAL/Pvf0Lo1rFwJd9wRO42IZJOiIvjTn/yttDR2\nmi2TJQ2Y7PfGG94d9OSTcMIJsdOISDbq0MFPL/vkE9hnn9hpqk4tgipYuhTOOQf+8Q8VARGpWMuW\nfmJhLlEhqIKnn4b27aFLl9hJRCTbzZsHjRvHTrFlVAiqYNQobxGIiFTk44/9pMJ27WIn2TIqBFUw\nZw7ssUfsFCKS7QYN8l1Jy8piJ9kylRYCM+tjZjtmIky2atjQN5oTEanIRRfBzjvDwQfD1Kmx01Rd\nVVoEOwPvmtkzZnaSWbYtok6/ww+Hl16KnUJEsl3z5jB8ONxwg08s+fjj2Imqpkori8tv/icAvYBD\ngGeAx0IIn6Q33trrR11Z/MkncMgh8PXXOnxGRKrm/vthxAgYNy7eHkQpXVlcfhf+pvytFNgRGGZm\nd29VyhwwdSp07w6HHRY7iYjkkj59/JyCt9+OnaRyVRkjuNrMJgJ3A+OB/UMIvYEOwP9Lc76o3nkH\nEgm47DLfUkKtARGpqqIi35esGrtCZ1xVVhY3BH4eQlhviUQIoczMTk1PrPiWLfOVxAMHwmmnxU4j\nIrmoVSuYMiV2ispVWghCCH0r+NqM1MbJHo8/7svFVQREpLoWLIAddoidonJaR7AZo0drm2kRqb4Q\nYORIOOaY2Ekqp0KwGSUl0KZN7BQikqtGjPBt6zt3jp2kctELgZk9ZmZzzezD2FnW1aABLF4cO4WI\n5KrFi327ibFjYyepXPRCAAwETowdYkOHHOLzf0VEquO882DoUD+57PHHY6epWFYcVWlmLYCRIYQD\nNvP1jC8omzQJTj4ZZszIjcEeEclOH30ERxwBb70Fe++d2WvrqMqt1LixnzL06quxk4hILttnH7jp\nJui72fmX8akQbMJ77/nU0WuugW7dYqcRkVx37rm+X1m2HmGZM0dVFhcXr30/kUiQSCTScp2FC301\nYP/+0LVrWi4hIgVm1119pfHChdCkSfquk0wmSVZjKXO2jBG0xMcI9t/M1zM2RnDzzfDdd/D3v2fk\nciJSAELwscZPP4WddsrcdXNmjMDMhgBvAXub2X/NrFfMPM88A1ddFTOBiOSbyZN93DGTRWBLZEWL\noDKZahF8/z3svjssWgQ1opdIEckHIcCZZ/phNbfemtlrV7VFoEKwjqVL/XShH35I+6VEpEC8/joc\ne6wvMKtbN7PXzpmuoWxSt66vKC4piZ1ERPJF27bQsSNccgmsWhU7zaapEKzDDHr0gAcfjJ1ERPJF\no0bw2ms+CeX662On2TR1DW1gwQKv4KNG+TYTIiKpsHAhtGvn6wkO2OQeCqmnrqFqWrgQatXKnUOn\nRSQ3NGwIvXrBk0/GTrIxFYJ1zJoFRx8Nv/2tdxGJiKRSIuE7F2QbFYJyIcAVV8B118HFF8dOIyL5\nqG7d7JyVmDNbTKRbSQlMn+4nk4mIpMO0ab4JXbZRi6DcuHFw4olQs2bsJCKSj0LwcwlOPjl2ko2p\nEJRbujR7l3+LSO578EE/uvLMM2Mn2Zi6hso1bAgffBA7hYjko9JS39Z+3DjYZpvYaTamFkG5447z\ns0WzcSBHRHJbURH85S9+LsFnn8VOszEtKFvHOef4vuH33pv2S4lIAbr3Xnj+eXjzTd/JIN206Vw1\nfPstHHSQV26dTCYiqVZWBvvuCwMG+DnG6aaVxdVQUgIrVvh4gYhIqtWo4YtVX3wxdpL1qRCUmz8f\nfv5zGDQIjjoqdhoRyVctWsAXX8ROsb7ohcDMTjKzmWZWYmY3xcpxzz0+ratLl1gJRKQQLFiQfVPV\noxYCM6sB/AU4EWgLnG1mbWJkefZZuPTSGFcWkUIybhx06BA7xfpitwg6ArNCCLNDCKuAp4GumQ5R\nVgb//S/st1+mrywiheTNN3290hlnxE6yvtiFoCkwZ52Pvyj/XEaVlvo839LSTF9ZRApFSQl07w6P\nPgrbbRc7zfpyZmVxcXHx2vcTiQSJRCJlP7tWLZ/SNXEiHH54yn6siMha99/vh9efckr6rpFMJkkm\nk1v8fVHXEZjZYUBxCOGk8o9vBkII4a4NHpf2dQT33AMTJsBzz6X1MiJSoMrKfBFZJhaSrZEr6wje\nBVqbWQszqwX0AEbECHLllTBjhs4rFpH0qFEjs0VgS0TtGgohrDazXwJj8aL0WAhhRowsder4WQSt\nWkH79nDkkTFSiIhkXvQxghDCS0BWHNXw9tvQrBm0bh07iYhI5mivoXILFkCbNvDyy77fkIhIrsuV\nMYKsMWiQrypWERCRQqNCUO6tt+Ckk2KnEBHJPBWCcvPnQ+PGsVOIiGSeCkG5Fi3g889jpxARyTwV\ngnLHHAPDhsVOISKSeSoE5Xr08AVlo0bFTiIiklkqBOW23RaGDoVevXxhmYhIodA6gg2MGwedO8OS\nJbD99hm5pIhIWmgdQTU99RScc46KgIgUjuhbTGSTkhIYPhw++SR2EhGRzFGLYB3Dh8NZZ0G9erGT\niIhkjgrBOkpKoG3b2ClERDJLhWAd9evD0qWxU4iIZJYKwTratfM9h0RECommj65jyRLfauKtt3xL\nahGRXJb100fN7Awzm2pmq82sfawc66pXD/r1gzPPhLlzY6cREcmMmF1DU4CfAa9HzLCRSy6BM86A\njh3hgw9ipxERSb/oXUNm9hpwXQjh/Qoek7GVxWucdRYUFcGQIRm9rIhIylS1a0gLyjZh4UJ49VUN\nHItIYUhrITCzl4Gd1/0UEIBbQwgjt+RnFRcXr30/kUiQSCRSkHDTnnkGTjgB9tknbZcQEUm5ZDJJ\nMpnc4u9T19AmXHwxdOgAvXtn7JIiIimX9bOGNlBp0EyrWTN2AhGRzIg5fbSbmc0BDgNGmdmYWFk2\n1Lw5TJ8eO4WISGZE7xqqikx3DU2b5mMEs2ZBnToZu6yISErlWtdQVmnbFo4+Gm68EXKgToqIbBW1\nCDZj4UI/0P7ww+Gee3RQjYjkHrUItlLDhpBMwg8/+L5DCxfGTiQikh5qEVTi44/hkEPg889hhx2i\nRBARqRa1CFLk4Yfhl79UERCR/KVCUIk33oAuXWKnEBFJHxWCSixfDg0axE4hIpI+KgSVaNrUxwlE\nRPKVCkElTjoJhg6NnUJEJH00a6gSS5bAvvvCoEG+rkBEJFdo1lCK1KsHAwdC9+7wyiux04iIpJ4K\nQRUcfzw89RT06gXnnafzjEUkv6gQVFHnzr4j6bvvwtNPx04jIpI6GiPYAitXQpMmvjtp06ax04iI\nVExjBGkwdSrsvruKgIjkFxWCLbB0KdSvHzuFiEhqxTyh7G4zm2FmH5jZc2aW9bfYvfaCmTN1RoGI\n5JeYLYKxQNsQwkHALOCWiFmqZNddoVkzGD06dhIRkdSJVghCCK+EEMrKP5wANIuVZUvcfjtcdx18\n/33sJCIiqZEtYwQXAllzeH1FunaFE0/0HUnnzImdRkRk6xWl84eb2cvAzut+CgjArSGEkeWPuRVY\nFUIYUtHPKi4uXvt+IpEgkUikOm6VmMH998Mf/wgHHQTnnw833AC77RYljojIWslkkmQyucXfF3Ud\ngZldAFwCHBtC+LGCx2XFOoINzZ4Nl10GzZtD//6x04iIrK+q6wjS2iKoiJmdBNwAHFVREchmLVr4\nofadOsVOIiJSfTHHCB4CtgdeNrP3zezhiFmqbfp0P9NYRCRXRWsRhBD2inXtVKpRA0pLY6cQEam+\nbJk1lLMOOADeeSd2ChGR6lMh2EoXXQQPPACrVsVOIiJSPSoEW+nYY6FVK7jmGm09ISK5SdtQp8Ci\nRV4Q9twTHnoIdt658u8REUk3bUOdQQ0awPjxvkV1mzZ+itlHH8VOJSJSNWoRpNj8+dCzJ7RrB/fd\nFzuNiBQytQgiadwYdtoJDj44dhIRkapRIUiDxYu9u0hEJBeoEKRBs2bw5ZexU4iIVI3GCNJg4UKo\nXdvfRERiqeoYgQqBiEie0mCxiIhUiQqBiEiBUyEQESlwKgQiIgUuWiEwszvMbLKZTTKzl8xsl1hZ\nREQKWcwWwd0hhANDCAcDLwJ9I2aJqjqHTeeSfH5++fzcQM+vUEQrBCGEpet8WBcoi5Ultnz/x5jP\nzy+fnxvo+RWKaEdVApjZ74Hzge+BY2JmEREpVGltEZjZy2b24TpvU8r/expACOG2EEJzYDDQJ51Z\nRERk07JiZbGZ7Q6MDiHsv5mvxw8pIpKDqrKyOFrXkJm1DiF8XP5hN2DG5h5blSciIiLVE61FYGbD\ngL3xQeLZwOUhhK+jhBERKWBZ0TUkIiLx5NTKYjPrY2Yzyged+8XOkw5mdp2ZlZlZw9hZUsXM7i7/\nvX1gZs+ZWf3YmVLBzE4ys5lmVmJmN8XOk0pm1szMxpnZtPK/t6tiZ0o1M6thZu+b2YjYWVLNzBqY\n2bPlf3fTzOzQih6fM4XAzBLAacD+5YPKf4qbKPXMrBlwPN5Vlk/GAm1DCAcBs4BbIufZamZWA/gL\ncCLQFjjbzNrETZVSpcC1IYS2QCfgyjx7fgBXA9Njh0iTB/AJOPsCB1LBGCzkUCEAegP9QgilACGE\nbyPnSYf7gRtih0i1EMIrIYQ1CwYnAM1i5kmRjsCsEMLsEMIq4Gmga+RMKRNC+CaE8EH5+0vxG0nT\nuKlSp/xF18nAo7GzpFp5i/unIYSBACGE0hDC4oq+J5cKwd7AUWY2wcxeM7NDYgdKJTM7HZgTQpgS\nO0uaXQiMiR0iBZoCc9b5+Avy6Ea5LjNrCRwEvBM3SUqtedGVj4OkewDfmtnA8q6v/mZW4XmJUVcW\nb8jMXgZ2XvdT+C/qNjzrjiGEw8zsJ8AzQKvMp6y+Sp7fr/FuoXW/ljMqeG63hhBGlj/mVmBVCGFI\nhIhSDWa2PTAMuHqDbWFylpmdAswNIXxQ3uWcU39rVVAEtAeuDCG8Z2Z/Bm6mgv3csqoQhBCO39zX\nzOxy4Pnyx71bPqC6UwhhQcYCbqXNPT8zawe0BCabmeFdJxPNrGMIYV4GI1ZbRb87ADO7AG+KH5uR\nQOn3JdB8nY+blX8ub5hZEV4EBoUQhsfOk0JHAKeb2clAbaCemT0RQjg/cq5U+QLvXXiv/ONhQIWT\nGXKpa+gS8FuFAAABr0lEQVQFym8iZrY3UDOXikBFQghTQwi7hBBahRD2wH+RB+dKEaiMmZ2EN8NP\nDyH8GDtPirwLtDazFmZWC+gB5NvskwHA9BDCA7GDpFII4dchhOYhhFb4721cHhUBQghzgTnl90mA\nzlQyKJ5VLYJKDAQGmNkU4Ed8s7p8Fciv5upDQC3gZW/wMCGEcEXcSFsnhLDazH6Jz4iqATwWQqhw\nZkYuMbMjgJ7AFDObhP+b/HUI4aW4yaSKrgIGm1lN4FOgV0UP1oIyEZECl0tdQyIikgYqBCIiBU6F\nQESkwKkQiIgUOBUCEZECp0IgIlLgVAhERAqcCoGISIFTIRCpBjM7xMwmm1ktM6trZlPNbL/YuUSq\nQyuLRarJzO7ANy2rjW/ydVfkSCLVokIgUk3l+7i8CywHDg/6Y5Icpa4hkeprBGwP1AO2i5xFpNrU\nIhCpJjMbDjyFnwi1WwihT+RIItWSS9tQi2QNMzsPWBlCeLr8IPvxZpYIISQjRxPZYmoRiIgUOI0R\niIgUOBUCEZECp0IgIlLgVAhERAqcCoGISIFTIRARKXAqBCIiBU6FQESkwP1/Kj1LWuBjKF0AAAAA\nSUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f5f6f5c0b38>"
       ]
      }
     ],
     "prompt_number": 159
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}