{
 "metadata": {
  "name": "",
  "signature": "sha256:790cfec8fbb91c33d8fc13b7a54a4a052916c1b24d03ebcc37b780d188dd16d4"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# grid\n",
      "nx=np.int32(5000)\n",
      "h=np.float32(0.005)\n",
      "\n",
      "# time\n",
      "nt=np.int32(5000)\n",
      "dt=np.float32(0.001)\n",
      "\n",
      "# velocity\n",
      "v=np.zeros(nx,dtype=np.float32)\n",
      "v[:]=4.0\n",
      "\n",
      "# source position\n",
      "delta=np.zeros_like(v)\n",
      "delta[nx/2]=1.\n",
      "\n",
      "# source wavelet\n",
      "w=np.fromfile(\"wavelet.bin\",dtype=np.float32)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(np.linspace(0.,nt*dt,nt),w,'k-')\n",
      "xlabel(\"Time (s)\",fontsize='large')\n",
      "ylabel('Amplitude',fontsize='large')\n",
      "xlim([0,5])\n",
      "savefig('wavelet.eps')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAESCAYAAADuVeJ5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFqtJREFUeJzt3X9w5HV9x/HnJrnkcpdcD44fx4+Tk5/1UIq/kKpIhl8V\ntKi0VZla1BkcbaXi4DgWtVyqTpXpFFuHVnHEFqjFGbUKFKliNRU7BQFBFLwCKeHHHRxc7k7JXXYv\n5LZ/fHYvm73d5Lu738139/t9PmYy2f3mm+w7O5DXfd6fz/fzBUmSJEmSJEmSJEmSJEmSJEkd5KvA\nVuAXC5zzBeAR4OfAy5eiKElS5zmNEAL1AuM84Lulx68B7lyKoiRJnWk99QPjS8A7Kp5vAg5td0GS\npKAn6QIacATwZMXzp4AjE6pFkjKnmwIDIFf1vJhIFZKUQX1JF9CAzcC6iudHlo7Nc8wxxxTHx8eX\nrChJSolx4NiFTuimEcbNwEWlx6cCOwmrquYZHx+nWCzu+7jzzjA3ftNNN807noWPjRs3Jl5Dp3z4\nXvhe+F4s/AEcs9gf4U4aYdwInA4cRJir2AgsK33tGsIKqfOAR4FdwHuj/NDnnnsOgB07dsRbrSRl\nTCcFxoURzrmk0R9aDgoDQ5Ja000tqaZMT0/P+5wlIyMjSZfQMXwv5vhezPG9aEzqA6NQKAAGRtb5\nXszxvZjje9GY1AdGPp+f91mS1JxMBMbw8HAmRxiSFKfUB0ahUGD16tWOMCSpRakPjHw+b2BIUgxS\nHxjlEYYtKUlqTeoDwxGGJMUjM4HhCEOSWpP6wLAlJUnxSH1g5PN5Vq1axczMTNKlSFJXy0RgDA0N\nGRiS1KLUB0ahUGB4eNjAkKQWpT4wyld679mzJ+lSJKmrZSYwHGFIUmtSHxi2pCQpHqkPDEcYkhSP\nTASGq6QkqXWpDwxbUpIUj9QHhi0pSYpHLukC2qBYLBbLD+jt7aVQKDAwMMDs7Cy5XBp/ZUlqTelv\n44J/IFM9wpiZmaGvr49ly5aRy+WYnZ1NuiRJ6lqpDox8Ps/AwAAAy5Ytsy0lSS1IfWAsX74cMDAk\nqVWpDozy3AUYGJLUqlQHhiMMSYpPqgOjUCgYGJIUk1QHhpPekhSf1AdGeYTR399vYEhSC1IdGNUt\nKe+JIUnNS3VgVLak+vr6vHBPklqQ+sAojzB6e3t54YUXEq5IkrpXqgOjsiXlCEOSWpPqwKhuSTnC\nkKTmpT4wbElJUjxSHRi2pCQpPqkOjMqWlCMMSWpN6gPDEYYkxSPVgVHZknKEIUmtSXVgeOGeJMUn\n9YFR2ZJyhCFJzUt1YNiSkqT4pDowbElJUnxSHRiOMCQpPqkODJfVSlJ8Uh8YXrgnSfFIdWC4NYgk\nxaeTAuONwCbgEeBjNb4+AvwauK/08cnFfqC71UpSfPqSLqCkF7gaOAvYDNwN3Az8quq8/wLOj/pD\n3a1WkuLTKSOMU4BHgQlgBvg68JYa5+Ua+aG2pCQpPp0SGEcAT1Y8f6p0rFIReC3wc+C7wIbFfqiT\n3pIUn05pSRUjnPMzYB2wGzgX+A5wfK0TR0dHAdi+fTv33nsv69atc4QhSRXGxsYYGxtr6HsaavG0\n0anAKGHiG+ByYC9w5QLf8xjwSmB71fFisRjyZ3h4mC1btjA8PMwVV1xBb28vGzdujLVwSUqDXC4H\ni2RCp7Sk7gGOA9YD/cA7CJPelQ5l7pc5pfS4OizmcWsQSYpPp7SkXgAuAb5HWDF1LWGF1PtLX78G\n+EPgT0vn7gbeudAPnJ2dZXZ2lmXLlgFhDiOfz7eleEnKgk4JDIDbSh+Vrql4/A+lj0gKhQIDAwPl\nYZbXYUhSizqlJRW7fD7P4ODgvue2pCSpNakNjOnp6X3XYIDLaiWpVakNjMqrvMERhiS1KjOB4QhD\nklqT6sBwDkOS4pPqwHCEIUnxyUxguKxWklqTqcCwJSVJzUttYLisVpLildrAcIQhSfHKTGA4wpCk\n1qQ6MFxWK0nxSXVgOMKQpPhkJjAcYUhSazIVGI4wJKl5qQ0Ml9VKUrxSGxi2pCQpXpkJDEcYktSa\nVAdG5bLa3t5eRxiS1IJUB4YtKUmKT2oDw0lvSYpXagNj165drFy5ct9zW1KS1JpGA6MHOKwdhcSt\nOjBsSUlSa6IGxgHAvwJ5YLx07HzgM+0oKg5TU1MMDQ3te25LSpJaEzUwvgT8BjgKKJSO/Q/wznYU\nFQdbUpIUr76I551JaEXNVBx7Djgk9opiYktKkuIVdYSxEzi46tiLgC3xlhOfqamp/UYYtqQkqXlR\nA+MrwDeBM0rf87vAdcA1baqrJTMzM+zdu5eBgYF9x2xJSVJrorakrgSmgauBZcA/EeY1/r5NdbWk\n3I7K5XL7jhkYktSaqIFRJIRDRwZEtep2FLi9uSS1aqHAOJMQFIv5YUy1xGbXrl3zltSCIwxJatVC\ngXEt8wPjSGAvMAmsIcxlPAkc3bbqmlS9QgoMDElq1UKBsb7i8ccJIfGXwG5gBfApYHvbKmtBvZaU\ngSFJzYs6h3EZcDiwp/R8NyFEtgB/3Ya6WrJz505Wr14975jLaiWpNVGX1e4CTqk69urS8Y6zc+dO\nDjjggHnHenrCr7p3794kSpKkrhd1hPFJ4DbgFuApYB3wZuCDbaqrJTt27NhvhAFzbalyeEiSoov6\nl/MG4DXAJmAV8CvgVOD6NtXVklojDLAtJUmtiDrCAHiIMNHd8Xbu3Mm6dev2O+5KKUlqXtTAuKHq\neeVy24tiqiU2O3bs4KSTTtrvuCulJKl5UQNjnBAS5b021gJ/AHytHUW1qt4chi0pSWpe1MAYrXHs\nK3WOJ+7ZZ5/lkEP233ndlpQkNa+V5UL3A6fHVUicnn76aQ47bP87ydqSkqTmNXIDpcp5i5WEu+09\nGHtFMXjmmWdYu3btfsdtSUlS86IGRvW+UrsII4wLY68oBgMDA6xYsWK/47akJKl5UQNjfTuLiFut\ndhQYGJLUiqhzGPfVOX5PXIXE6dhjj6153DkMSWpe1MCo9Rc4RwdubQ5wwgkn1DzuHIYkNW+xllT5\ngr0BwjYguYqvrSfeSe83An8H9BKW7F5Z45wvAOcSdst9D3VGPieeeGLNF7AlJUnNWywwxkufi6XH\nuYrnPwG+EVMdvYT7hZ8FbAbuBm4m7FlVdh5hpHMcYV+rLxL2s9rPaaedVvNFbElJUvMWC4zR0uc7\ngf9oYx2nAI8CE6XnXwfewvzAOB+4rvT4LmA1cCiwtfqH1ZvDsCUlSc1bKDDeAPy49HgGOKPOeXHc\n0/sIwu1ey54ijCIWO+dIagRGLperPgREa0mNj48zMTHB9PQ0xWKUW5pLUjYsFBj/CLy09Lj6OoxK\nL46hjqh/mauToOb3jY6O7ns8MjLCyMgIsHBLqlgscvHFF3PrrbeyYcMGBgcHvW+GpNTatm0bk5OT\nDX3PQoHx0orH65spqAGbCTdlKltHGEEsdM6RpWP7qQyMSgu1pG655RbuvvtuHnvsMQYHB6NVLUkp\nUa8zU6lT/gl9D2Eyez3QD7yDMOld6WbmtlI/FdhJjXbUQhZqSV1//fVceumlhoUk1bHQCOPJBb5W\nVgReFEMdLwCXAN8jrJi6ljDh/f7S168BvktYKfUoYWuS9zb6IvVaUsVikTvuuIOrrrqqqeIlKQsW\nCow/WbIqgttKH5WuqXp+SSsvUK8l9cQTT9DT01PzLn2SpGChwBhbqiKWSr2W1MMPP8yGDRsi9fAk\nKauizmEMAJ8mtIN2lz5/Bljeprraol5gTExMsH79+qUvSJK6SNTdar8IHA/8OfAEYd7iE4RrIxqe\nS0hKX19fzZbU448/zlFHHZVARZLUPaIGxluBY4AdpecPEq62HqeLAmOhEcY555yTQEWS1D2itqSe\nBqrvSDQIbIm3nPaqFxibN2/miCOOSKAiSeoeUUcYNxBWMF1NWG77IuDPCDvYVm4ZEsc2IW1Tb1nt\n5OQkBx10UAIVSVL3iBoYHyh9vrziWK50/AMVx+LYJqRt6i2rnZycZM2aNQlUJEndI5W3aK2nVkuq\nWCwyOTnJgQcemFBVktQdOmVrkCVRqyW1e/ducrkcK1ZUT9FIkipFDYyTCfMTOwhbnZc/9rSprrao\n1ZKyHSVJ0URtSd0IfBO4FJhuXzntVasltX37dgNDkiKIGhhrgSuIft+KjlSrJbV9+3bnLyQpgqgt\nqeuBP25nIUuhVkvq+eefZ9WqVQlVJEndI+oI47OE+3pfDjxbcbxI/Vu3dpxaLampqSmGhoYSqkiS\nukfUwPgGYRuQbwP5iuNd1aKq1ZIyMCQpmqiBcTJwEFBoYy1tV6slZWBIUjRR5zDuADa0s5ClYEtK\nkpoXdYQxAXwf+Df2n8O4Iuaa2qa3t5eZmZl5x6ampjjkkEMSqkiSukfUwFgB3Eq4kdI65uYuettR\nVLv09fWRz+fnHZuamuLoo49OqCJJ6h5RA+M9Vc9PAi6iy5ba1ltWa0tKkhbXyF5SBwMfBu4rfbwa\n+FA7imoX5zAkqXmLjTD6gfOBdwO/BzxEWGJ7FPB2YGtbq4uZy2olqXmLBcYzhEnuG4DLgEdKxy+h\ny67BAJfVSlIrFmtJPUC4u95rgFOA4bZX1Eb1WlLDw139a0nSklgsMEYI11/cA/wVYbRxEzBEaFd1\nlVotqV27dnkvDEmKIMqk9wTwKeBY4GxCaOwFfg78Tdsqa4NaLalCocDy5csTqkiSukejd9z7CfA+\nwnbnlwAvjb2iNqrVksrn8wwMDCRUkSR1j2Zv0TpNuKnSuTHW0na1WlKFQsHAkKQIMnVP71otqXw+\nb0tKkiLIXGBUjjDK4dHXF/WCd0nKrkwHhhPekhRdpgKjr69vXkvKCW9Jii5TgVFrhGFgSFI0mQ4M\nJ7wlKbpMBUb1slpHGJIUXaYCo3pZrSMMSYouc4HhCEOSmpOpwLAlJUnNy1Rg2JKSpOZlLjAcYUhS\nczIVGNUtKUcYkhRdpgKjuiXlCEOSostcYLiXlCQ1J9OB4V5SkhRdpgKjevNBW1KSFF2mAsO9pCSp\neZkODEcYkhRdJwTGgcDtwMPA94HVdc6bAB4A7gN+2swL1brS2xGGJEXTCYHxF4TAOB74z9LzWorA\nCPBy4JRmXqjWld6OMCQpmk4IjPOB60qPrwPeusC5uVZeyJaUJDWvEwLjUGBr6fHW0vNaisAPgHuA\n9zXzQl7pLUnN61ui17kdWFvj+CeqnhdLH7W8DngaOLj08zYBd9Q6cXR0dN/jkZERRkZGAK/0lqSy\nsbExxsbGGvqello8MdlEmJt4BjgM+BHw24t8z0ZgCvjbGl8rFou1M2f37t2sWbOG6elpAC644ALe\n9a53ccEFFzRXuSSlRC6Xg0UyoRNaUjcD7y49fjfwnRrnrACGS49XAucAv2j0hWq1pBxhSFI0nRAY\nnwPOJiyrPaP0HOBw4NbS47WE9tP9wF3AvxOW4DbElpQkNW+p5jAWsh04q8bxLcCbSo//Dzi51Rfq\n6emhWCxSLBbJ5XJOektSAzphhLFkcrncvKW1jjAkKbpMBQbMb0t5pbckRZfJwCiPMJz0lqToMh0Y\ntqQkKbrMBUblPTGc9Jak6DIXGI4wJKk5mQ4MRxiSFF3mAqN8tXexWGTPnj309/cnXZIkdYXMBUZ5\nWW05LHp6MvcWSFJTMvfXstySckmtJDUmc4FRbkk54S1JjclcYJRbUk54S1JjMhkYjjAkqXGZC4zK\nlpQjDEmKrhO2N19S5ZbUzMyMIwxJakAmA2N2dtbAkKQGZa4lVbms1paUJEWXucAobz7opLckNSZz\ngVG5SsoRhiRFl9nA8EpvSWpM5gLDK70lqTmZCwyv9Jak5mQyMBxhSFLjMhcYtqQkqTmZDIyZmRny\n+TyDg4NJlyNJXSNzgdHf378vMJzDkKToMhkYhULBwJCkBmUuMAYGBtizZ4+BIUkNylxg9Pf3GxiS\n1ITMBoarpCSpMZkMDOcwJKlxmQsM5zAkqTmZCwznMCSpOQaGJCmSTAaGcxiS1LhMBoYjDElqXOYC\nw0lvSWpO5gLDEYYkNSeTgVEoFLyntyQ1KJOB4QhDkhqXucConMNwaxBJii5zgdHf38/U1BQ9PT30\n9fUlXY4kdY1MBsbk5CTDw8NJlyJJXSVzgTEwMMC2bdsYGhpKuhRJ6iqZC4zh4WG2bdvGypUrky5F\nkrpKJgMDcIQhSQ3qhMD4I+BBYBZ4xQLnvRHYBDwCfKzZF1u1ahVgYEhSozohMH4BvA348QLn9AJX\nE0JjA3Ah8JJmXqzcihocHGzm27vK2NhY0iV0DN+LOb4Xc3wvGtMJgbEJeHiRc04BHgUmgBng68Bb\nmnmxnp7wK2dhSa3/M8zxvZjjezHH96IxnRAYURwBPFnx/KnSsaatWLGipYIkKWuW6p/ZtwNraxz/\nOHBLhO8vxlsOHHjggXH/SElKtVzSBVT4EfAR4Gc1vnYqMEqYwwC4HNgLXFnj3EeBY9pQnySl2Thw\nbNJFRPUj4JV1vtZH+GXWA/3A/TQ56S1J6l5vI8xPTAPPALeVjh8O3Fpx3rnA/xJGEJcvZYGSJEmS\nMiiWC/tS4qvAVsI1Llm2jtDqfBD4JfChZMtJ1HLgLkI79yHgs8mW0xF6gfuItvAmzSaABwjvxU+T\nLWVp9BJaVeuBZTjHcRrwcgyMtcDJpcdDhJZmlv+7KK8l7wPuBF6fYC2d4DLga8DNSReSsMeASMtG\nu+U6jMXEdmFfStwB7Ei6iA7wDOEfDwBTwK8Ic2NZtbv0uZ/wj6ztCdaStCOB84Cv0FmrRZMS6T1I\nS2DEfmGfUmc9YdR1V8J1JKmHEKBbCa26h5ItJ1GfBz5KWJ6fdUXgB8A9wPsWOjEtgRH7hX1KlSHg\nm8ClhJFGVu0ltOiOBN4AjCRaTXLeDDxL6Nk7uoDXEf4xdS7wQUJLu6a0BMZmwgRn2TrCKENaBnwL\n+BfgOwnX0il+TViy/qqkC0nIa4HzCb37G4EzgOsTrShZT5c+Pwd8m9DiTzUv7Nvfepz0zhH+EHw+\n6UI6wEHA6tLjQcLu0GcmV07HOJ1sr5JaAZTvV70S+G/gnOTKWTpe2DfnRmALUCDM7bw32XIS83pC\nG+Z+QvvhPua2l8malxG23bmfsITyo8mW0zFOJ9urpF5M+G/ifsLS86z/7ZQkSZIkSZIkSZIkSZIk\nSZIkxeeXhC01lsIG4O4I510CfK7NtUiSqkwBz5c+9hJ2dy0/v3CJa/kW8PYI5w0QLs48uL3lSJLq\neYywx1ASDgMmCdvcRPFl4CPtK0cK0rL5oLQUJpgLkVHgG8ANwG8I220cR9haYSvwOHB2xff+FnAt\nYcuWp4BPU///v7OBe4E9Fcc+Vvq+3xDuLFkZZmPAm5r5haRGGBhSdNXb6L+ZsLnhAYR9qm4vHT+c\nEAjXVJz7z4QAOIawlfQ5wMV1XudlhH3Ryk4gbDv9KmBV6XsnKr6+CfidRn4RSVJ8arWkKo+NAt+r\n+NrvE+Y6yvdXGCbMg6wCDgXyhPtql10I/LDOa3+Z+ffcPpYwajmTsF17teOAF+r+JlJM+pIuQOpi\nz1Y8nga2MTcKmS59HiLcsGgZc/cdgDC6f6LOz93B3JbTEHZg/jAhpE4kBNVlFT9vmHCPC6mtbElJ\n7fckYav5NYT21QGEOY2X1Tn/AeD4qmM3Eu6EdhQhlK6s+NpLmLt3udQ2BobUfk8D3weuIowGeghz\nGfWu6fgB8ArmVkkdT2iFDRCCJw/MVpx/OnBb7FVLVQwMqTlF9p8EX+j5RYQAeAjYTlhhtbbOz95K\nmN94a+n5AGFO4zlC+BzE3I1ulhNuHnZdw7+BJCkVXgL8NMJ5XuktSZIkSZIkSZIkSZIkSZIkSZIk\nSZLUjf4frhlRg80ad4kAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10907e490>"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Python, Numpy"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#import numpy as np\n",
      "def model1d(nx,nt,h,dt,v,delta,w):\n",
      "    um=np.zeros_like(v)\n",
      "    uo=np.zeros_like(um)\n",
      "    ud=np.zeros_like(um)\n",
      "    up=np.zeros_like(um)\n",
      "\n",
      "    h2=h*h\n",
      "    dt2=dt*dt\n",
      "    seismo=np.zeros((nt,nx),dtype=np.float32)\n",
      "    \n",
      "    for it in xrange(nt):\n",
      "        for ix in xrange(1,nx-1):\n",
      "            ud[ix]=(uo[ix-1]+uo[ix+1]-2.*uo[ix])/h2\n",
      "        for ix in xrange(nx):\n",
      "            up[ix]=(ud[ix]+delta[ix]*w[it])*v[ix]*v[ix]*dt2+2.*uo[ix]-um[ix]\n",
      "        um,uo,up=uo,up,um\n",
      "        for ix in xrange(nx):\n",
      "            seismo[it,ix]=uo[ix]\n",
      "    return seismo\n",
      "\n",
      "def model1d_vec(nx,nt,h,dt,v,delta,w):\n",
      "    um=np.zeros_like(v)\n",
      "    uo=np.zeros_like(um)\n",
      "    ud=np.zeros_like(um)\n",
      "    \n",
      "    h2=h*h\n",
      "    dt2=dt*dt\n",
      "    seismo=np.zeros((nt,nx),dtype=np.float32)\n",
      "    \n",
      "    for it in xrange(nt):\n",
      "        ud[1:-1]=(uo[:-2]+uo[2:]-2.*uo[1:-1])/h2\n",
      "        up=(ud+delta*w[it])*v*v*dt2 + 2.*uo - um\n",
      "        um,uo=uo,up\n",
      "        seismo[it,:]=uo[:]\n",
      "    return seismo"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "seismo=model1d_vec(nx,nt,h,dt,v,delta,w)\n",
      "seismo=seismo.clip(seismo.min()*0.01,seismo.max()*0.01)\n",
      "plt.imshow(seismo,aspect=3,extent=[0,(nx-1)*h,(nt-1)*dt,0],cmap=cm.gray_r)\n",
      "xlabel('Distance (km)',fontsize='large')\n",
      "ylabel('Time (s)',fontsize='large')\n",
      "savefig('seismo.eps')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8k1X///FX05a2zBtQ9qggsoeMsvdUUVT0VhRB0FuQ\nDYIKiMiP4a0UkI2oDNnIxj1u5Ct7dEChlFlkL1EECtIkvz+uXiWUtDlpM64kn+fjwYM2PSc5hPST\nk/d1rnOBEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhAgQHYBDwBHgHS+PRQghhJsFA0eBSCAUiAMq\ne3NAQggRqEweepwotMKfDNwBlgOdPPTYQgghbHiq8JcETtl8fzrtNiGEEB7mqcJv9dDjCCGEcCDE\nQ49zBiht831ptFl/uoIFC1qvXr3qoeEIIYR/KFq0KBcuXAhypo9TjXMgBEgCWgNngV1AFyDRpo11\n9OjRHhqOsf3666+0aNHC28MwBG8+F7t37+bbb7+957aWLVvSrFkzr4xHXhd3yXNx15gxY8DJWu6p\nGX8q0A/4AW2FzxfcW/SFMIw7d+4wd+5cLl++fN/PNm3axLZt2xg0aBDh4eFeGJ0QOeepjB/gO6Ai\n8DDwoQcfVwhlly5dYsKECXaLvu727dt89NFHnDlzxoMjE8J1PFn4haLIyEhvD8EwPPlc7Nmzh1mz\nZim3//zzz9myZYsbR3QveV3cJc9Fzngq6hFOkBf1XZ54LlJTU5k7dy6XLl1yuu8vv/zCli1bGDx4\nMGFhYW4Y3V3yurhLnouckRm/CGiXLl1i/Pjx2Sr6utu3b/Pf//6Xs2fPunBkQriPFH4RsPbu3etU\ntOPIZ5995tHoR4jskqhHBJzU1FQ+++wzLl68qNS+a9euhISEsGDBAodtf/nlF7Zt28bAgQPdHv0I\nkV0y4xcBRY92VIv+kCFDqFatGpUrV+btt9/GZHL8K5OSkiLRjzC0YG8PwMYHckKGcKeYmBiWLFmi\n1LZWrVr06tWLQoUKYTKZCAoKInfu3DRq1Ig7d+5w6tQph/cRExNDaGgoZcqUyenQhcjU5s2bAcY4\n00eiHuH3UlNT+eKLLzh//rxS+5dffpkqVaqQK1eue24PCgoiT548PPHEEzzyyCPMnz/f4X39/PPP\nbN26lUGDBt13f0J4i0Q9wq/p0Y5q0R8yZAjVq1fPskiHhoZSqVIlhg0bpnSfKSkpfPjhh5w7d06p\nvRDuJlGP8FvORDs1a9a8J9pxRI9+GjduzD///KMU/ezdu1eiH+FyEvUIgRbtzJs3T3mG/dJLL1G1\nalWnoxg9+unYsSOPPPKI0qqfn3/+me3btzNgwACJfoTXSNQj/Mrly5cZP368ctEfPHgwNWrUyFER\nDg0NpXLlygwbNoygIMebJN64cYMPP/xQOX4SwtUk6hF+IzY2lsWLFyu1rVGjBr169aJw4cJK0Y4j\nttHPrVu3OH36tMM+e/fuJSwsjNKlSztsK0RmJOoRAclsNjN//nzl3TK7dOmSHu2ozNBV6dHPU089\nRcWKFVm4cKHDPj/++CPbt2+nX79+Ev0Ij5GoR/i0K1euMG7cOOWiP3jwYGrWrElYWJhLi76t0NBQ\nqlSpwtChQ5Xa//333xL9CI+SqEf4rLi4OBYtWqTUtlq1avTu3dtl0Y4j+uy/cePGpKSkKL0x7d27\nl/DwcEqVKuX28Qn/IVGPCAipqaksXLhQKUcHLdqpVq0aoaGhbpvl26Pn/p06daJixYp8+eWXDvv8\n8MMPbN++nb59+0r0I9xGoh7hU65cucL48eOVi/6gQYOoWbOmy/N8VUFBQfdEPypjuHbtGh9++CEX\nLlzwwAhFIJLCL3xGfHw8M2bMUGpbtWpV3n//fUqUKEFwsPcTzeDgYIoUKcLo0aOJiopS6jNnzhx2\n7tzp5pGJQCRRjzA8s9nMwoULlc6OBXjxxRepXr26x6MdR2yjn0qVKilFP99//z07duygT58+hIaG\nemCUIhDIjF8Y2h9//MG4ceOUi/6gQYOoVauW16IdR4KCgsiVK5dTq37+/PNPJkyYINGPcBkp/MKw\n9u3bx/Tp05XaVqlSxVDRjiN69PPBBx9Qr149pT5z5sxh165dbh6ZCAQS9QjDMZvNLF68mOTkZKX2\nL7zwAjVq1DBctOOIHv0888wzVKpUSWlp6nfffcfOnTvp3bu3RD8i22TGLwxFj3ZUi/6gQYN49NFH\nDRvtOKKv+qlWrZpy9PPHH38wYcIE5auICZGRFH5hGM5EO5UrV/apaMcRk8mUHv3UrVtXqc/s2bMl\n+hHZIlGP8Dpno51///vf6Ttq+uIsPzN69PPss89SqVIlpQ3nvvvuO3bt2kWvXr0k+hHKZMYvvOrq\n1atORzu1a9d261473qRHP9WrV+ett95S6nPlyhWJfoRTpPALr9m3bx/Tpk1TalupUiVGjRrlN9GO\nIyaTiaJFizJ69Gjq1Kmj1Gf27Nns2bPHzSMT/kCiHuFxZrOZpUuXcvz4caX2zz//vFe3XfAWfaO3\nZ599lsqVKytFP9988w27du3ijTfeICREfr2FfTLjFx6lRzuqRX/gwIHUqVPHb6MdR/QTvqpVq8aQ\nIUOU+ugXmL906ZKbRyd8lRR+4TEJCQnK0U6FChUYNWoUJUuWDIhox5Hg4GCKFSvG6NGjqV27tlKf\nWbNmSfQj7PJk4Z8HXAD2e/AxhQGYzWaWLFnC6tWrldo/99xzdOvWjXz58gXkLD8zevTTuXNnXn75\nZaU+33zzDbNnz+bOnTtuHp3wJZ4s/POBDh58PGEAf/75J+PGjePo0aNK7QcOHEjdunUJDw+Xom+H\nHv1Ur15dOfq5ePEiEyZM4PLly24enfAVniz8vwFXPfh4wssSEhKYOnWqUluJdpxjG/3UqlVLqc/M\nmTPZu3evm0cmfIEc9hcuZzabWb58ufIs/7nnnnP7dXD9kR79PP/881StWpUlS5Y47PP111+zZ88e\nevbsKSd8BTA5uCtc6q+//pJox4Nso5/Bgwcr9Tl//rxEPwHO079pkcBGoLqdn1mbN29+t2FkJJGR\nkZ4ZlXCJAwcOsGrVKqW25cuX56WXXpIDuC5ktVq5ceMG69atIz4+XqnPk08+qbxKSBhDcnLyPWe6\np11s3alfIkMV/tGjR3t2NMIlzGYzK1eu5PDhw0rtO3fuTK1atSTacQOr1cqdO3c4cOAAS5cuVepT\nokQJXn31VYl+fNSYMWPAyVruyahnGbANeAQ4BfTw4GMLN9GjHdWiP2DAAOrVqyfRjpvo0U+NGjWU\no5+zZ88yYcIErly54ubRCaPwZOHvApQAwoDSaMs7hQ87ePAgn3zyiVLbcuXKyaodDwoODqZ48eKM\nHj2a6tXtfcC+34wZM4iJiXHzyIQRyKoe4TSz2cyqVas4dOiQUvvOnTv79MVSfJW+6ufFF1+kWrVq\nLFu2zGGfjRs3EhMTw6uvvip7/fgxWdUjnKJHO6pFX492JM/3Dj36qVmzpnL0c+bMGcaPHy/Rjx+T\nwi+UORPtPPTQQ4waNYpSpUpJtGMAevTz/vvvU61aNaU+M2bMIDY21s0jE94gn+WEQ2azmdWrV5OY\nmKjU/tlnn+XRRx+VWb7BBAUFkTdvXrp06UJCQoJS9LNhwwZiY2Pp1q2bRD9+RGb8IkvXrl1j3Lhx\nykV/wIABREVFyaodg8pO9HPq1CnGjx/PH3/84ebRCU+Rwi8ylZiYyJQpU5Tali1blvfee0+iHR+R\nnehn+vTpxMXFuXlkwhPks5u4j8ViYfXq1Rw8eFCp/dNPPx3QF0vxVbbRz/79+1m+fLnDPuvXrycm\nJkaiHx8nM35xj2vXrjF27Fjlot+/f38aNGgg0Y6P0qOfWrVqSfQTQKTwi3SHDh1SjnbKlCnDe++9\nR+nSpSXa8QO20U+VKlWU+kj047vks5rAYrGwZs0aDhw4oNReoh3/pEc/L730Evv372fFihUO+6xf\nv564uDi6du0q0Y8PkRl/gNOjHdWi369fP+rXry/Rjp8KCgoiLCyMRx99lEGDBin1OXnyJOPHj+fq\nVbnOkq+Qwh/AkpKSlKOdUqVKpUc7MrPzf8HBwZQoUYL333+fSpUqKfWZNm0a+/btc/PIhCvIb3AA\nslgsrFu3jv371a5736lTJ+rWrSvRToDRo5+uXbuyb98+Vq5c6bDP2rVriYuL46WXXpIJgoHJjD/A\n/P3334wdO1a56Pfr14+GDRtKtBOg9Oindu3aDBw4UKnPiRMnJPoxOCn8AeTw4cNMnjxZqW3JkiUZ\nOXKkrNoRgBb9lCxZklGjRjkV/ahOMIRnyWexAGA2m1m/fr3yL+FTTz0lO2qK+wQFBZEvXz5eeeUV\n4uPjlaKfNWvWEBcXR5cuXST6MRCZ8fu569evM27cOOWi37dvXxo1aiTRjrBLP+HLmejn+PHjjB8/\nnj///NPNoxOqpPD7scOHDzNp0iSltiVKlGDkyJGUKVNGoh3hkG3088gjjyj1mTp1qkQ/BiGfvfyQ\nxWJh/fr1ykvrJNoR2aFHP926dSM+Pp6vvvrKYZ81a9YQHx9Ply5dZILhRTLj9zPXr19n7NixykVf\noh2RE/qqnzp16ihHP8eOHWPcuHES/XiRFH4/cuTIEeVop1ixYhLtCJexjX4qVKig1Gfq1KkkJCS4\neWTCHol6/IDFYmHDhg3Ex8crtX/yySeJioqSaEe4lB79dO/enbi4OFatWuWwz+rVq4mPj+eFF16Q\nVT8eJDN+H6dHO6pFv0+fPhLtCLfRo5+6desqRz9Hjx6VVT8eJoXfhzkT7RQpUoQRI0ZQtmxZmVkJ\nt9P3+hk1ahTly5dX6iPRj+dIBfBBFouFr7/+mtjYWKX2HTt2pH79+hLtCI8ymUzky5ePHj16OBX9\n7Nu3jxdeeEGOPbmRzPh9jB7tqBb9Pn360LhxY4l2hFfYRj8DBgxQ6nPkyBHGjRvHX3/95ebRBS4p\n/D7k2LFjytHOAw88INGOMIzg4GBKlSrFqFGjKFeunFKfTz75RPkSoMI5UhF8gNls5ttvvyUmJkap\n/eOPP07Dhg0l2hGGoq/66dGjB7GxsaxZs8Zhn6+++oqKFSvy3HPPyQTGhWTGb3D6XjuqRf/NN9+k\nadOmEu0IQwoKCiI8PJyoqCjl6CcpKYnx48dL9ONCUvgN7Pjx48rRTuHChRkxYgSRkZEyMxKGZxv9\nREZGKvX55JNPSExMdO/AAoRUCAOyWCx8++237N27V6n9Y489RqNGjSTaET5Fj35ee+015ehn5cqV\nVKxYkeeff15W/eSAzPgN5saNG4wdO1a56Pfu3ZtmzZpJ0Rc+yTb66d+/v1KfpKQkxo0bx7Vr19w8\nOv/lycJfGtgEHAASALWAL4AcP36c6OhopbaFChVixIgRPPTQQ4SEhEjRFz5Nj37ee+89ypYtq9Rn\nypQpHDp0yM0j80+ejHruAIOBOCAvsBf4CQj40M5isfDdd9+xZ88epfaPPfaYXAdX+B2TyUT+/Pl5\n/fXXiYmJYe3atQ77rFixgsqVK9O5c2eJfpzgaMb/IPAW8D/gCpCa9vcvwNC0n6s6j1b0Aa6jFfwS\nzgzWH928eZOxY8cqF/0333yTZs2aERERIUVf+J3sRD+JiYkS/Tgpq8L/XyAGqAh8DrQFKqf9PQ94\nJO3n/83G40YCjwI7s9HXb5w4cYKJEycqtS1QoICs2hEBIyQkJD36KVOmjFKfKVOmkJSU5OaR+Yes\npoz9gM+A21m0CQdeB2Y48Zh5gV+BccA6m9utzZs3T/8mMjJSeZmXr7FYLHz//ffs3r1bqX2HDh1k\nR00RkKxWK7dv32bv3r2sW7fOcQegatWqPPPMM34b/SQnJ5OcnJz+/ebNmyHrWn4fT1eRUOBr4Dvg\nkww/s44ePdrDw/G8mzdvKs/yQVu1U7ZsWYKDg6Xoi4BlNps5ffo0M2aozzEHDx5M/vz53TgqYxgz\nZgw4WctVV/W0AvQNNooDXwLzgWJOPFYQ8AVwkPuLfkBITk52KtoZPny4rNoRAm3VT+nSpXnvvfco\nVaqUUh+JfjKnWvhnoR3YBZiMthrICsx14rEaA12BlkBs2p8OTvT3WfqqnYULFyq1b9++PUOHDqVg\nwYKYTHKqhRBwd9VPr1696NSpk1Kf5cuXs3r1asxms5tH51tUjxKWAH5Hi2raA2XRsv9zTjzWFgLw\nhLEbN24wadIkrFarUns92pEDuELcT1/107BhQ0qXLq0U/SQkJJCQkBAw0Y8K1UJ8DS3WaYZ2Atbf\naNFNqJvG5ReSk5OJjo5WKvr58uW7J9oRQmROj35GjhxJyZIllfpMmTKFw4cPu3lkvkG1wkwHdgFh\nwKC02xojJ1/ZZbFY+Omnn9ixY4dS+3bt2tG0aVPZdkEIJ5hMJgoUKEDv3r3Zs2cP69evd9hn2bJl\nVK9enU6dOvntqh8VqoX/I7Sll2bgaNptp9GWcgobN2/eVJ7lA/Tq1UvW5guRTXr006BBA+XoZ//+\n/ezfv58hQ4aQL18+D4zSeJzJ3JO4W/QBDgP7XTsc33by5EkmTpyoVPTz5MnDu+++S7ly5aToC5FD\nISEh6dFP8eLFlfpMnjw5YKOfrAr/buDfQK5Mfp4r7ee7XD0oX2OxWPjxxx9ZsGCBUvu2bdvyzjvv\nUKhQIVm1I4SL6NFPnz59eOqpp5T6LFu2jLVr1wbcqp+spprdgbHAbLQN1ZLQ9tjJh7ZdQ220PXy6\nu3mMhibRjhDGkXHVz8yZMx322bdvH/v27Quo6Cer6eZBoDNQDVgM3AIKAzeAhUBVtBl/wB7g/f33\n35WjnYiICIl2hPCQkJAQypQpw4gRIyhatKhSn8mTJ3P06FHHDf2ASgU6h3amrkhjsVj45Zdf2LZt\nm1L7Nm3a0KxZM9lrRwgPMplM/Otf/6Jfv37s2rWLjRs3OuyzZMkSatWqRceOHf161Y9MPZ2UkpKi\nPMsHeOONN2RtvhBeokc/jRo1onTp0syaNcthn7i4OOLi4vw6+pEji044deoUH3/8sXK0M3z4cMqX\nLy9FXwgvCwkJoWzZsowcOZIiRYoo9Zk8eTLHjh1z88i8w0ifZT5o0aKFt8dglx7tbNiwQal969at\n6dq1K/ny5ZNVO0IYRFBQEGFhYdSpU4c8efIoLeXct28f165do0KFCoaNadO2ZR7jTB+ZijrgbLTz\nn//8Rw7gCmFQevTTuHFjypQpoxT9xMbGEhsby1tvvUXevHk9MEr3U52OmoA30JZv6idtNUNb1eO3\nnIl2wsLCePfdd3n44Yel6AthcHr0M2LECB58UO0KspMmTfKb6Ee18I8BXkO7Ipd+HbQzwLvuGJS3\n6dHOvHnzlNq3bt2akSNHyglZQvgQfdXPgAEDeOKJJ5T6LF68mA0bNmCxWNw8OvdSnZr2QLtG7iW0\nvfkBTnD34ix+IyUlhejoaOX/2Ndffz092jFqBiiEsE+Pfpo0aULZsmWdin6GDh1Knjx5PDBK13Mm\n6rme4bY8aNsz+43Tp0/z8ccfKxV9PdqpUKECoaGhUvSF8GG20U/hwoWV+kRHR3P8+HE3j8w9VAv/\nd2hX3gq36TcWcHxGhA+wWCz873//44svvlBq36pVK0aMGCHRjhB+RI9+Bg0axOOPP67UZ9GiRWzc\nuNHn9vpRjXqGAAuAP9EuvnId+BHo5p5heU5KSgqTJ08mNTXVcWPuRjuhoXINGiH8jR79NG3alLJl\nyzJ79myHfWJiYoiJifGp6Ed1uvoX8AzaJRcbAuWBp9GuzOWz9GhHpeiHhITwzjvvpEc7Qgj/FRIS\nQmRkZPonexXR0dGcOHHCzSNzDWdP4ApG26wtCG2Xzny4Luf32AlcFouFX3/9lXXr1im1b9myJd27\ndyd//vwS7QgRIPTZf7169QgPD1fawC0+Pp4bN25Qvnx5j9UKd57A1Rb4FIjMcLsVY53969CtW7eY\nPHkyd+7cUWr/2muvpW+7IAdwhQgsevFv1qwZkZGRStHPnj172Lt3L2+99ZZhox/Vt6TPgQlAAbQL\nsOh/wtw0Lrc4c+YMH330kVLR16OdRx55RFbtCBHg9Ohn+PDhFChQwGF7q9Vq6OhHdcYfDsxHu+au\nz7FarWzevFn/SORQixYtaNmyJREREVLwhRCAtuqnYMGCDB06lK1bt/L999877PPll18SFRVF+/bt\nDRUTq47kE+BttGzfp9y6dYsPP/xQuei/9tprtG/fXoq+EOI++kZvzZs3p3fv3kp9du3axdixY7lx\n44abR6dOtfCvAnqhreI5YfPH0GcvnD17VjnaMZlM6dGO5PlCiMwEBQUREhLCQw89xPDhw8mfP79S\nv+joaJKTk907OEWqUc9qYDPaG0CK+4bjGhaLhS1btrBp0yal9s2aNaN169YyyxdCKMsY/fzwww8O\n+yxcuJD69evTrl07r0Y/qoU/Em2vHsNn/M6u2unZs2f6jppS9IUQztBX/bRo0YLIyEg+/fRTh312\n7tzJzp07vXrCl+pbznqglTsH4grnzp1TjnaCg4N5++23qVixoqzaEUJkmx79lCtXjuHDhytfrjE6\nOpqTJ0+6eXT2ObOqZwPwf8BFm9utGGDbBovFwtatW/nf//6n1L5p06a0adNGoh0hhMvo0c+wYcPY\nsmULP/74o8M+CxYsoEGDBrRt29aj0Y9q4T+Q9icjtctSudGtW7eYMmUK//zzj1J7iXaEEO6iRz8t\nW7YkMjKSuXPnOuyzY8cOduzYwbBhw8idO7cHRmms5ZnW0aNHO9Xh3LlzSk8saO/GQ4cOpXDhwoZa\nTyuE8E8Wi4U///yTadOmKS/l7NGjB2XKlHHc0MaYMWPAyVqe1Yy/GVq0A1nn+2r5ihYXbUY72zcX\n2nGD4Yp972GxWNi2bRu//PKLUvsmTZrQtm1biXaEEB6jRz/vvPOOcvQzf/58GjZsSOvWrQkOdt9u\nOFkV/llAtbSv55F5rPOQ4mPdAloCN9MedwvQJO1vZbdu3WLq1KncunVLqX2PHj2oUKGCRDtCCI/L\nTvSzfft2duzYwdChQ90W/WSVeVQDuqR9HYlW4O39ccbNtL9zoW3u9ocznc+fP89HH32kVPRNJhNv\nv/02lSpVklU7Qgiv0Vf9lC9fnuHDhysVc6vVysSJE/n999/dMiZHYbdagO7c48UBF4BNwEGVTlar\nla1btyqtkQVo3Lgx77//Pg888IDk+UIIQ9Cjn3fffZc2bdoo9Zk/fz4//fSTyy/u7mga/Dfanvuu\nVgD4AXgX+DXtNmvz5s3TG0RGRhIZGcnt27eZOnUqKSlqJwy/+uqrsu2CEMLQUlNTOXHihPLiFCB9\n1U9ycvI9Wz+k7UPmVLFz1Pgm0NFBG9WDuxmNQtv+ITrt+/tW9Zw/f155lg/w9ttvy6odIYRPsFgs\nXL161amJrb1VP65e1QPaShxHVyBXzfkfAFLRrtsbgXZxF7tXjbFarWzbto2ff/5Z6Y4bNWpEu3bt\nyJ07t8zyhRA+wWQyUahQId59911+++03pXo3f/58GjVqROvWrXM0wXVU+G/g/AHczBQHFqLl/CZg\nEXDfeszbt28zffp05XWvEu0IIXxVUFAQERERtGrVinLlyilFP9u2bWPHjh289dZb2V71462M3x5r\n7969mTNnjnKHYcOGyQFcIYRf0KOfKVOmcPv2baU+PXv2ZN68eeDijP86kNeZO8wB5e0fGjZsSPv2\n7SXaEUL4FavVyq1bt/j111+V9x5L49LC70lKhV+iHSGEv0tNTeX48eN89tlnql38t/BLtCOECBR6\n9BMdHU1qaqqj5v5X+Bs0aECHDh0k2hFCBBQ9+tm0aZOjKwr6V+Hv3r07FStWlGhHCBGwFKIfp4qj\noTITq9V635+2bdvKXjtCiIAWEhJC27Zt7dbI7DBSNbVm9o/Yu3cvmzdvdvl+FUIIYXTBwcF06NCB\nypUr2/152qTYd6MeR+9eM2bMUN6OWQghfF2ePHl48803s2zj94UfYNmyZZw5c8YDwxFCCO956KGH\n6Ny5s8N2AVH4AXbv3s1vv/0m0Y8Qwu8EBwfTvn17qlSpotQ+O4XfUAd3ly9frtSuXr16DBkyhPDw\ncDePSAghPCdPnjwMHjxYuegvWrQoW49jqMJ/+PBhZsyYody+X79+lChRwo0jEkIIz4iMjHSY5+vu\n3LnD1KlTOXHiRLYey31X83XeB82bNyclJYVdu3YRHh6uVNSrV69OSEgIp06dyvbSJiGE8BY92mnR\nooVS+y1btrBy5UquX7+O1WrVL8Rid4v7zDjaltkrUlJS+P777zly5Agvvviiw/ZRUVFERUXJqh8h\nhE9RWbVja9GiRZw8eRKz2Zyjia6hoh64exJXamqqRD9CCL9VtmxZ5aJvNpuZNm0aJ06cIDU1Ncfp\nhqGinhYtWqSfoav/nZ3oJzg4mNOnT0v0I4QwHJPJRLt27WjZsqVS+y1btrB8+fL0aCcjv4h6rFbr\nfdszpKSk8N1333H06FFeeOEFh/dRv3596tWrx+zZs5WvZSmEEO7mbLSzePFikpOTcxztZGS4qAe4\n7x9otVoxm80kJSUxc+ZMpfswmUz07duX4sWLu2OIQgjhFGeindTUVKZNm8bx48ddXvTBgFGPrYwz\nf6vVmh79REREKBX1GjVqYDKZJPoRQnhFcHAwbdu2pVWrVkrtt27dyooVKzKNdjLyi6jHlr1/tF78\nv/32W44ePcq///1vh/fToEED6tevz8yZM2XVjxDCY3Lnzk2fPn2U2zsb7WR312JDRj2O6NHPoUOH\nlKOfoKAg+vXrR7Fixdw8OiGEgDJlyigXfX3VjrPRTnZTDENHPY7os//du3eTJ08epaJeo0YNgoKC\nOHPmjEQ/QgiXCw4Opk2bNrRu3Vqp/bZt27JcteOI30U9KqxWKzdv3uTrr7/m8OHDStFPw4YNiYqK\nYvbs2RL9CCFcxtloZ8mSJZw4ccItB3Cz4pNRT0b6CV/ORD/BwcES/QghXKZ06dJORTvTp0/n2LFj\nLjkhy1k+HfVkZBv95M2bVzn6AST6EUJkS3BwMK1bt6ZNmzZK7bdv386yZcv4+++/s6w5qgduAzLq\nycg2+kmOO7JvAAASwUlEQVRKSlKKfho1akT9+vWZNWsWt2/f9sAohRD+ICIigr59+yq3dybacedE\n1C+inoyyG/30799foh8hhJLSpUsrF32LxcK0adNcGu1kdykn+FnUk5HtCV/58uWjaNGiDvvUqFED\nq9Uql3cUQthlMplo1aqVcrSzY8cOli5dmu1VOzrbQh8UFJR+XxL12KEX/40bN5KUlMTzzz/vsE/j\nxo2Jiopizpw5Ev0IIdJ5c9WOvo+ZvoNxTvhl1JOR1Wrlzp07JCYmMmvWLKU+oaGh9O/fX+lTghDC\n/5UqVUq56FsslvRVO65aqmk7y88pv456MtIP/O7evZv8+fMrFfWaNWtisVg4e/asrPoRIgCZTCZa\ntGhBu3btlNrv3LmTpUuXOly14yrZiXo8PeMPBmKBjR5+3HR68d+wYQOrVq1S6tOkSRMGDhxIWFiY\nm0cnhDCS3LlzM2TIEOrUqaPUfunSpfz000/cunUrx0U/JwdvHfF04R8IHAS8OnXOTvSjr/qR6EeI\nwFCyZEmn8vzp06dz9OhRl56Q5a7i78nCXwp4HPgccN9bmRPMZjOXLl1i4sSJJCQkKPV55ZVXqF+/\nPiZTQBweESLg6NFOly5dlNrv2rWLjz/+mD/++AOLxeKycbjiIG5mPLmqZwowDMjvwcd0SI9+1q1b\nx6FDh3juuecc9mnatCn169fn008/lVU/QvgRZ0/IWrZsmUsP4HqKp6atHYGLaPm+IWb7tvQTvg4e\nPMjs2bOV+uTKlYv+/ftTpEgRN49OCOEJJUuWVC76VquVGTNmcOTIEZeekOXOXP+ex/LIo8AE4BUg\nFQhHm/WvBrrZtLE2b948/ZvIyEgiIyM9NLy7goKCyJ07N48//jhVq1ZV6vN///d/7Nmzx6Uf84QQ\nnhEcHEyTJk2oV6+eUvvdu3ezadMmUlJSXJrlq95XcnIyycnJ6d+nrepxqpZ7Y/bdHBgKPJnhduvo\n0aO9MJz7BQUFERISQqVKlejcubNSn9u3bzN37lyJfoTwIREREfTp00d5pr18+XKOHj1qqGhnzJgx\n4GQt99YRSmM8Y5nQV/0cOHBAOfoJCwuT6EcIH1KiRAn69u2rXPRnzJjB4cOHvbKNsqt5o/BvBp7y\nwuM6zWKxcPHiRSZOnMjBgweV+nTr1o169ep5LKsTQjjHZDLRvHlzXnrpJaX2e/bs4eOPP+bKlSt+\nE+f6/V49OaWv+lm7di2JiYlK0U/z5s2Jiori888/l+hHCAPxh2jHFWQxuoLsrPqJiIigf//+PPjg\ng24enRBCRXajHX8r+iCFX5nVasVsNqdHP4mJiUr9unfvTt26deWELyG8xGQy0bRpU+VoJyYm5p5o\nx9+KPkjU4zQ9+lmzZg2VK1fm2WefddinRYsWREVF8cUXX0j0I4QHhYeH06dPH+WJ14oVKzhy5Ihf\nzvJtyTQ0G7Kz6id37twS/QjhQcWLF6dfv37KRX/mzJn3RTv6SVX+tlhDCn8O6NFPdHQ0SUlJSn26\nd+9OnTp1JPoRwk1MJhNNmjTh5ZdfVmqvRzuXL19OL/r+WOxtSdSTQ1arlRs3brB69WoqV67MM888\n47BPy5YtiYqKYt68eRL9COFCOY129GJvG/Pot7nyQijeJtNOF7Barfzzzz8kJCQwZ84cpT558uSh\nf//+PPDAA24enRCBwdloZ8aMGSQlJaWfkOXosoa28Y+vk8LvQmazmQsXLhAdHc3hw4eV+rz66qvU\nrl1boh8hsik4OJjGjRsrRzuxsbF2T8jKajafsej7evYvUY+L6dHPqlWrqFKlCk8//bTDPq1atUqP\nfv755x8PjFII/xAREcGbb76pPHFauXKlUztq2hZ2/RNBxmLvixGQTDPdQF/1s3//fuXoJ2/evAwY\nMECiHyEUFS9enL59+zq1aufQoUNOF33bGCizGb6vzfyl8LuJfsKXHv0cOXJEqZ8e/fjaC0kITzGZ\nTDRs2FA52omLi0tftePMVa30trZFP6u+vvQ7K1GPm+nRz8qVK6latapy9FO3bl0WLFgg0Y8QNpxd\ntfPVV1+RlJSUoxOy7M34M5v9+0rsIzN+D8hO9JM/f34GDBhA4cKF3Tw6IXxDsWLFnD4h69ChQy45\nCzdjoXc0uzf67F8KvwfZRj/Hjh1T6tOjRw9q1apl+BeSEO6iRztdu3ZVah8fH8/EiRO5fPlyjvfa\nybiGP+N92btv2z5GJVGPh+nRz4oVK6hatSqdOnVy2KdNmzbUq1ePhQsXSvQjAkp4eDi9e/cmJESt\nVLki2rGV1QFdUJ/5Gy3+kRm/F+gnfO3bt49PP/1UqU+BAgUk+hEBpWjRovTr10+56M+aNeueaCcn\n6+wzztptC7c+8zdaMXeGFH4vso1+jh8/rtSnR48e1KxZ09AfI4XICZPJRIMGDXjllVeU2u/fv5+J\nEydy6dKlLLdecIbez1GBV31DMNrvq0Q9XmaxWLhx4wbLly+nWrVqPPWU46tStm3blrp167Jo0SKJ\nfoRfCQ8Pp1evXoSGhiq1X7Vq1X2zfE/OxjOu4nEmFvLmJwaZ8RuAHv3Ex8crRz8FCxaU6Ef4FT3a\nUS36s2bNIjExEbPZDDhXhO3Jzqw8q8czchQkhd9AzGYz58+fJzo6mhMnTij1kehH+DqTyUT9+vWV\no52EhIR7oh24N5LJKtvP6nZXFeqMJ35lxpu/s34T9dh7kn3lZApb+qqfZcuWORX91KlTh8WLF0v0\nI3xKWFgYvXr1IleuXErtV61adc+OmpD5rpkZi6+9k7Bsb88pe+NwZr2/J2uV38z47S278rWir9NP\n+IqPj2fu3LlKfQoVKsSAAQMoVKiQm0cnhGsULVqU/v37Kxf9WbNmcfDgQbt77Tia+NlbW69/MnB1\nncjpeQOe+CTgU4Xf3n9axndMf4k89L1+zp07x6RJkzh58qRSv549e1KjRg2/eR6E/zGZTERFRSlH\nOwcOHEiPdjI7ISvj7N/2a5XIxdW/L5kdYDZK9GPYqMfeu7Xt1xnzPNv/eF+MeDJjtVq5fv06S5Ys\noXr16jz55JMO+7Rr147atWuzdOlSiX6Eobgi2lGR1Uw/q/auYq9OqYzDtr87a5lhZ/yOTsCwfWIz\nzvRt+/rDzFePfuLi4pSjnwceeECiH2EoRYoUcTraSUxMzHbR9xZ7xTo7Y7It/hn/5JRhC7+qzA6k\n2J5MYa+Nr70h5CT6qV69us/9e4X/MJlM1KtXj27duim1P3jw4H0nZGWH6mw5p2f5OrrfnHDXcUvD\nRT2ZHZnPrK29wp7xI5K9Ezu8PSvILtvop0aNGnTs2NFhn/bt26dHP3fu3PHAKIXQhIWF8Z///Ifw\n8HCl9qtXr05fm5/Tg6Qq2yjrbW3buGInz4xfZxZdO8tVb06Gn/Fnd/tTXzhtOrv06Cc2NpbPPvtM\nqc+DDz7IwIEDJfoRHqNHO6pFf/bs2Rw4cMDpaMeejJO9rM7ozWzdvTP1Iqs3lYyRjSvk9H4MV/hV\n98VQYZuJueLd1kgyRj+nTp1S6tezZ0+qVavmN8+DMB6TyUTdunWVo53ExEQmTpzIxYsX77n4eU7Z\nm3ln/Brurwm2BVr19ySz2uTsyh5PMXzUo/ozuP9jlb2DvrY/9wcWi4Xr16+zaNEiatasyRNPPOGw\nT4cOHXj00UdZsWKFrPoRLhUWFsbrr79ORESEUvs1a9Zk6wBuTtmrBZnFQfY2fLN3W2Z9MjvOmFM5\nqWOenvEnA/uAWGCXvQY5yeDtvVNn9Y7rL7Nefa+fmJgY5einaNGiDBgwgIIFC7p5dCJQPPjgg/Tv\n31+56M+ePZuEhATu3LnjtYmYyjYPOts2mcVGtvfriX9TdmuYpwu/FWgBPApEZdrIRQdZHMU8Rp31\nJycnZ6uf2Wzm7NmzTJo0idOnTyv1ee2116hataph3wRVr1QWCIz6XJhMJurUqUP37t2V2iclJeU4\n2snu74g99l77RjnRyhFfmfEDKD1rrirKmR1Uye6BG3tnD7taTl7U+qqfL7/8km+++Uapz2OPPcbL\nL7+svCuiJxm12HmDEZ8L/YSsli1bKrVfu3Ytq1at4ubNmzn6HXdl4bc3O3d2UUnGWuCpNwdfmvH/\nDOwB/qPUwQWzftv/2MwioKzeGDIeGLKX3xltFuBs9FOsWDEGDhwo0Y9Qpkc7efLkUWo/Z84cr0c7\n9jjK67PqlzH+yZgyGJWnC39jtJjnMaAv0FSlkyuKv737cnR7Zkf2bYu/kf+TbaOfs2fPKvUxevQj\nvC+70c6FCxfSt1E2EmeP+Tn6fbe3c4CrMn/bY6A5uU9v/naPBq4Dk9K+jwNqem84QgjhkzajHTs1\npNxAvrSv8wBbgXbeG44QQgQmT67jLwqstXncJcCPHnx8IYQQQgghhLd0AA4BR4B3vDwWb0vGwUlu\nfmwecAHYb3NbIeAn4DDaJ8R/eWFc3mDvufgAOI322ohF+70JBKWBTcABIAEYkHZ7IL42MnsuPsDH\nXhvBwFEgEghFO8hb2ZsD8rITaC/oQNQUbdWXbbH7GHg77et3gP96elBeYu+5GA0M8c5wvKoYUCvt\n67xAElqNCMTXRmbPhVOvDSNs0haFVviTgTvAcqCTNwdkAIG6lvI34GqG254CFqZ9vRB42qMj8h57\nzwUE5mvjPNqEELSVgIlASQLztZHZcwFOvDaMUPhLArZbS57m7j8kEDl9kpufK4oWeZD2d1EvjsUI\n+gPxwBcERrSRUSTaJ6GdyGsjEu252JH2vfJrwwiF35hnP3lPtk5yCxBWAvv1Mht4CO2j/jnungMT\nKPICq4GBwN8ZfhZor428wCq05+I6Tr42jFD4z6AdsNCVRpv1B6pzaX9fQlv+mulmdgHiAlquCVAc\nuOjFsXjbRe4WuM8JrNdGKFrRXwSsS7stUF8b+nOxmLvPhVOvDSMU/j1ABbSPLbmAF4AN3hyQF2U8\nya0d9x7cC0QbAH1vgO7cfaEHouI2Xz9D4Lw2gtDii4PAJza3B+JrI7PnwidfG4+hHZ0+Cgz38li8\n6SG0AzdxaEu1Au25WAacBf5BO+7TA22F088E1pI9uP+56Al8ibbUNx6tyAVKpt0EsKD9XtguVwzE\n14a95+IxAve1IYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEcJ/ZwHveHkQOhKFtpauvr14AjHXR\nfe8EqrjovoQQwiOSgZvANbQdKbcCvcjeTpTJQCtXDcyF+gOzbL6fD/w/F93382j7tgghhM84wd1i\nnQ94EjiOdjGS7NxXaxeNy5USgIY238/HdTP+cOAKcramEMKH2BZ+XT3AzN0IYwF3C+UDwNdonw6u\nAP+H9ulgUVqfm2i7Nw5Na/8V2uZ3fwKbuTcWWQDMTLu/a2hb3Jaz+XlVtKs8XUHbC13fRsMEvIu2\nzchlYAVQMJN/X5m0MdnujWVb+POhXVVJ33tlAdqng2/T/h2/oW1GNjXt35zI3Ytx6H4EumXy+ELc\nwwibtAlhz260XVr1baltt919C23/mgeAImjF2Aq8AvwOdEQrptFp7b8BHgYeBGKAJRke6wW0S9cV\nRCvk49Nuz4e2F8y3aJtgPQz8kvaz/mgXAmmW9rOraG8g9lRH+wRjyXC7FSicdp+/AYNsfvY8MDLt\n3/gP2hvSbrT9aVYBkzPcVyJQM5PHF+IeUviFkZ3F/mUo/0ErtpFoM/ytDu5nAXAD7QpvY9AKpL4L\nqhVYg7ZLrBntTUGfTXdMG8OUtMe8zt3rIPdCO9h81uZ+n8P+79S/uH//eNAuOPQr2qeF921u18cU\nC9xG2577Bto2vFZgJdo1G2z9TWBsUiZcQAq/MLJSwB823+sHeyeizcx/BI6hXW81Mya0a7EeBf5C\ni5VAm0nrLth8nYJ2kQvQrg1xPJP7jUQryFfT/hwEUrGfs1/l7huN7b/lCbR8/lM7fWz3lr+V4Xvb\nMeryY/9SjULcRwq/MKp6QAlgi52fXUfL78ujxS1DgJZpP8t4FaaX09q0BgqgbX0NaiuGfufevD/j\nzzqgxUP6n9zcvZCOrX1pj2v7+2YFPgN+QIuSciuMJyuV0bbkFcIhKfzCKPRCnB8tYlmGdrD2QIaf\nk/bzh9Nuu4YW0ej5+QW0NwRdXrS45A+0i9tMyORx7fkGLVIaiLYOPx93r2w0J+2+yqR9/yDaG4w9\np9E+cdS387j90K5FsRFt9u9oTPaEA7XRDkIL4ZAUfmEUG9GK+O9oB2snoV2IRWd7cPdhtCL3N7AN\n7aDq5rSffYiWvV9F+yTwJXAS7RKfCcB27v1UYO9arfr3fwNt0ZaXnkO74EeLtJ9NRbsC1I9p495O\n1pe7+xTt4LO9x30D7c1hHdobTMYxZTVG0sa3CW3VkRBCCIPIxb1n7rrSDuTMXSGEEEIIIYQQQggh\nhBBCCCGEEEIIIYQQQgghhBBCCCGEEML7/j9KbIrtnHPY9wAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x161b94cd0>"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "seismo=model1d(nx,nt,h,dt,v,delta,w)\n",
      "seismo=seismo.clip(seismo.min()*0.01,seismo.max()*0.01)\n",
      "plt.imshow(seismo,aspect=1.0,cmap=cm.gray_r)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "<matplotlib.image.AxesImage at 0x10af756d0>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAARUAAAEACAYAAACd9eLKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH81JREFUeJztnXuQFOXd7z+zM7vgAsZbwuIlrFWAiubU+8ZUMF7yLhc9\nqxj1n8QYIZQxCQRSKhATkWAqf0TBVB1fL+gxoAZShiRVJpYpI1F4s+QkMUJy1MUYhCVuDKssWUR9\nT0FgF+b80d1sM/RtZrp7+vL9VE1Nz9PPdP92p+c7/XyfX/8ahBBCCCGEEEIIIYQQQgghhBBCiMTQ\nCWwDdgDfanAsQoiUUwR6gHagGXgFOK+RAQkhwqUp5v19EkNUeoFB4CfAtTHHIISIkLhF5QzgH7bX\nu8w2IURGiFtUyjHvTwgRM6WY99cHnGV7fRbG2cpRTj755PK+fftiDUoIAWPHjqW/v79Q73bq3kCV\nlIA3gOnA28Bm4Abgr7Y+5e985zsxh1U7XV1ddHR0NDqMqogi5tdee42nnnoKgJkzZ/KJT3witG3r\nfxwP3/3udyEETYh7+DMEfB34NfA68FOOFRSRQvbu3XtUUACeffZZdu/e3cCIRCOJW1QAngPOASYA\n9zRg/yJEhoaGeOihh45rf/TRRzl48GADIhKNphGikina29sbHULVhBnz9773Pdd1y5cvD2Ufef8f\npw2JSp2k8eAJK+af/vSnvn1Wr15d937y/D9OIxIVURNbt25l27Ztvv36+vr405/+FENEIilIVETV\nDAwM8POf/zxwfxm3+UKiIqpicHCQlStXVv0+Gbf5QaIiquLuu+92XXfRRRcxbdo01/VhGbci2UhU\nRGC8jNmmpiZmzpzJtGnTGD16tGu/xx57LIrQRIKQqIhA+BmzS5cuZeTIkYwYMYLFixe79tu1axdb\ntmyJIkSRECQqwhc/Y3bhwoWMGTPm6OtRo0bxjW98w7X/r371Kxm3GUaiIjwZGhryNGZnzZpFW1sb\nhcLwJSOFQoEPf/jD3HTTTa7vk3GbXSQqwhOvjNkpU6Zw/vnn09R0/GHU1NTEOeecI+M2h0hUhCt+\nxuzVV19NqeRePaNYLDJ9+nRP4zaMjFuRLCQqwpHXXnvN05j99re/zciRI32309LS4mncKuM2e0hU\nxHEMDAwcU8qgkoULF3qefVTiZ9wq4zZbSFTEMfhlzDoZs35Yxu2XvvQl1z4ybrODREUcg1fGrJcx\n60dTUxOTJk1i6tSprn1k3GYDiYo4ip8x+5nPfMbTmPWjWCwyY8YMz6HTqlWrat6+SAYSFQEEy5ht\naWmpez9+xu3bb78t4zblSFSEb8bsokWLGDNmTFU+ihejRo3i9ttvd10v4zbdSFRyThBjduzYsaEJ\nChjG7WmnnSbjNqNIVHKOXymDCy64oCZj1g/LuPXLuC2Xdf+5tCFRyTFBShkUi8XI9l8sFn1LJSjj\nNn1IVHJK0FIGUeNXKuHtt99WqYSUIVHJIXv37q2qlEHU+Bm3KpWQLiQqOePQoUOON/+yqCVjtl5k\n3GYLiUrOuOce95tCRmnM+iHjNjtIVHLET37yE9d1cRizflilEryGXjJuk49EJSds3bqVN954w3V9\n0FIGURMk41bGbbKRqOSAIMZsNaUMoqa1tVXGbYqRqGScwcHBxBmzfgQ1bg8dOhRjVCIoEpWME1Up\ng6gJYtzec889Mm4TSPKOJhEafsasX43ZRqMat+lEopJR/IzZuDJm60XGbfqQqGSQam/+lXR0c7J0\nIVHJGH6lDGbPnp04Y9YP1bhNFxKVjOFXymDy5MmJNGb9UMZteghydD0O9ANbbW2nAC8A24HngZNs\n65YAO4BtwBW29gvNbewA7q89ZOGGVymDYrGYeGPWD2XcpoMgovIE0FnRdgeGqEwCNpqvASYD15vP\nncDDgHWe/QhwMzDRfFRuU9RBd3e3bymDESNGxBhRNLS0tLBo0SLX9TJuG08QUfk/wL6KtmuANeby\nGuA6c/laYB0wCPQCPcAUYBwwBths9ltre4+ok7179/KLX/zCdX3SMmbrRaUSkk2tg+uxGEMizOex\n5vLpwC5bv13AGQ7tfWa7qBO/UgZpNGb9kHGbbMIYYJfNR2h0dXUdXW5vb6e9vT3MzWcKr1IGn/rU\npxKbMVsvhUKBSZMmMX36dDZu3OjYZ/ny5dx1112ZEtQw6e3tpbe3N/Tt1ioq/UAbsBtjaLPHbO8D\nzrL1OxPjDKXPXLa397ltvKOjo8aw8oVfxuxVV13V0FIGUVMsFpk6dSpbtmzhgw8+cOyzevVqvvKV\nr8QcWTqo/MHetGlTKNut9SfsGWCOuTwHeNrW/nmgBTgbw5DdjCE+H2D4KwVgtu09oga6u7tTUcog\nakaMGOFr3G7evNl1vQifIKKyDvgDcA7wD+AmYDlwOcaU8jTzNcDrwM/M5+eA+QwPjeYDqzGmlHuA\n9aH8BTkkb8asH36lEp577jkZtzESRFRuwDBaWzCGNk8A7wIzMKaUrwDes/W/G5gAnAv82tb+Z+Bj\n5rpb6g08r6SxlEHUqMZtssieg5dxvDJms2zM+mFl3E6fPt21jzJu4yF/R1+K8TJmC4UCV111Vaoz\nZuvFujmZMm4bi0QlJfgZs8uWLcuFMetHkIxbGbfRIlFJAQMDA57G7KJFi3JlzPrhl3Er4zZaJCoJ\nJ0gpg7Fjx+bKmPVDxm1jkagkHD9jNq2lDKJGpRIah47GBONlzJZKJWbOnJlrY9YPq1TCiSee6NpH\nxm34SFQSip8xe+edd2ailEHUtLS0sHDhQtf1Mm7DR6KSQPwyZmXMVoeM23iRqCSMIKUMZMxWh0ol\nxItEJUGUy+XcljKIGnupBDdk3IaDjs4E4VdjNuulDKLGKpUg4zZaJCoJwc+YTcvNv5LOiBEjfI3b\nl156KcaIsodEJQGolEG8+Bm369evl3FbBxKVBpPHGrONxjJub775Ztc+Mm5rR6LSQGTMNo5CocDE\niRNl3EaAjtYGkvcas43GKpUg4zZcJCoNoru7m+3bt7uuz0uN2UajjNvwkag0AL9SBjJm40UZt+Ei\nUYmZIKUMZMzGi0olhItEJUbK5bJKGSQUlUoIDx29MeJlzBaLRZUyaDBWqYQPfehDrn1WrVoVY0Tp\nRKISE37G7NKlS1XKIAG0tLRw2223ua5/5513ZNz6IFGJAZUySBcybutDohIxKmWQPoJm3P7rX/+K\nMar0IFGJEGXMppcgGbcrVqyQceuAjuYICWLMKmM2uQTJuJVxezwSlYiQMZsN/DJu33nnHZVKqECi\nEgEqZZAtVCqhOiQqIaNSBtlDpRKqQ6ISIjJms4tl3M6YMcO1jzJuDXR0h4hKGWQbq8atMm69kaiE\nhEoZ5IMgxm3eM24lKiGgjNl80draqoxbDyQqdaKM2fyhjFtvgojKWcBvgL8ArwG3mO2nAC8A24Hn\ngZNs71kC7AC2AVfY2i8Etprr7q8n8CQQxJhVKYNsEsS4zWvGbZCjfRBYCJwPXAQsAM4D7sAQlUnA\nRvM1wGTgevO5E3gYsH6mHwFuBiaaj84w/ohGsW7dOtd1KmWQfYIYtz/4wQ9ijCgZBBGV3cAr5vL/\nA/4KnAFcA6wx29cA15nL1wLrMMSoF+gBpgDjgDGA5WKttb0ndXR3d7Njxw7X9cqYzQd+xu3u3btz\nl3Fb7Xl5O/DvwEvAWKDfbO83XwOcDuyyvWcXhghVtveZ7alDxqyw42fc5i3jthpRGQ08BdwK/HfF\nurL5yDx+xuwXv/hFGbM5Q8btsQQd8DdjCMqPgKfNtn6gDWN4NA7YY7b3YZi7FmdinKH0mcv29j6n\nnXV1dR1dbm9vp729PWCY0eJnzF588cUyZnOK3bjdsGGDY58VK1awbNmyxBwfvb299Pb2hr7dIKJS\nAB4DXgf+09b+DDAHWGE+P21r/zHwvzCGNxMxfJQy8AGGv7IZmA084LTDjo6O6v6KmPAyZkulkjJm\nc45l3G7ZsoX333/fsc+qVauYO3duzJE5U/mDvWnTplC2G0QyLwFmAVOBl81HJ7AcuBxjSnma+RoM\n8fmZ+fwcMJ/hodF8YDXGlHIPsD6MPyIO/IzZO++8U8askHFLsDOV3+EuPm6T9Hebj0r+DHwswD4T\nhYxZUQ2Wcfv973/fcf369esZP348bW1tMUcWD8kY3CUYGbOiWvJu3EpUPJAxK2olaMbtkSNHYowq\nHvRt8ODHP/6x67pSqcSVV14pY1a4ktdSCRIVF1599VV6enpc18uYFUHIo3ErUXFgYGCAp59+2nW9\nVWNWPooIQmtrK9/85jdd12ct41aiUsGhQ4dYuXKl63rVmBXVUigUOPXUU3Nj3EpUbAQxZlVjVtRC\nU1NTboxbfTtseBmzhUJBGbOiLvJSKkGiYtLd3e1pzC5btkzGrKgbP+O2v78/9catRAVlzIp4ybpx\nm3tRUcasiJtCocBpp52WWeM216KijFnRKKyM2+nTp7v2Satxm+tvi1/GrIxZESXFYpFp06ZlzrjN\nraj4GbPKmBVxkEXjNpeiImNWJImsGbe5E5WDBw/KmBWJImvGba5EpVwus3z5ctf1l1xyiYxZ0RAs\n4/byyy937ZMW4zZX3x4vY7a5uZnOzk4Zs6JhFItFOjo6OOmkk1z7pMG4zY2o+JUyWLJkCSNHjowx\nIiGOp6Wlhdtuu811fRqM21yIil8pAxmzIkmk3bjNvKj4lTKQMSuSRtqN20yLil/GrIxZkVTSXOM2\n098mL2MWUI1ZkWisUglpM24zKyp+GbN33XWXMmZF4kmjcZtJUVHGrMgSaTNuMycqfhmzs2fPljEr\nUkXajNtMiUqQjFnVmBVpJE3Gbaa+XU8++aTrOt38S6Qdq1SCl3H76KOPxhiRM5kRlVdffZWdO3e6\nrlcpA5EFmpubPY3bPXv2NNy4zYSo+GXMLl68WMasyAxJN25TLyoHDx70zZj9yEc+ImNWZIagxu2B\nAwdijGqYVIuKShmIvBKkVMK9997bEOM21d82v1IGqjErskyQUgmNMG5TKypBShk0NzfHGJEQ8eOX\ncdsI4zaVohK0lIF8FJEHkmbcpk5UghizypgVeaJQKHDqqafy5S9/2bVPnMatn6iMBF4CXgFeB6w6\nAqcALwDbgecB+6BuCbAD2AZcYWu/ENhqrru/lmCDGLPnnXeejFmRO5qampgwYUIijFu/b9+/gKnA\nvwH/w1y+FLgDQ1QmARvN1wCTgevN507gYcA6ZXgEuBmYaD46qw3WK2O2WCxy5ZVXUiqVqt2sEJkg\nKcZtkJ/0/eZzC1AE9gHXAGvM9jXAdebytcA6YBDoBXqAKcA4YAyw2ey31vaeQHR3d3tmzC5dulQZ\nsyL3JMG4DSIqTRjDn37gN8BfgLHma8znseby6cAu23t3AWc4tPeZ7YEYGBhQKQMhAtJo4zaIqBzB\nGP6cCXwaYwhkp2w+IkHGrBDVYWXc+hm3UZVKqMaAeB94FsNw7QfagN0YQ5s9Zp8+4Czbe87EOEPp\nM5ft7X1uO+rq6gIMY/a3v/2ta0AXX3yxMmaFcKBQKDBhwgRmzJjBhg0bHPusWLGCT3/606H/IPt9\nG09jeGbnBOBy4GXgGWCO2T4HsJJGngE+j+G/nI1hyG7GEJ8PMPyVAjDb9p7j6OjooKOjg74+V92h\nVCopY1YID6xSCSeffLJrn23bth39voWFn6iMA/4Lw1N5CfglxmzPcgyB2Q5MM1+DMe38M/P5OWA+\nw0Oj+cBqjCnlHmC9146DlDJoaWnxCV+IfNPc3Mytt97quj4K4zaJRkR5wYIFnj7K4sWL5aMIEZBy\nuczAwAD33nuva5+5c+da0811f6mS+K30NH3nzJkjH0WIKimXy2zfvp3Vq1f7da1bE1L1zVTGrBC1\nYRm3Xhm3YZGab2dzczOdnZ0yZoWoESvj1su4DYPUiMqSJUuUMStEnbS0tHgat2GQClFRKQMhwsMv\n47ZeEi8qypgVIlyClEqoh0SKSrlcplwu8/7773PBBRfImBUiZJqamjj33HOPftfK5fCutEn0t/XE\nE09k/vz5jQ5DiMxRKBRYuHBhJNtOtKgAjBo1iilTpjQ6DCEyRUdHR2QjgMSLCsBll13GuHHjGh2G\nEJngox/9KBdeeGFk20+FqADceOONnHDCCY0OQ4hUM2rUKD73uc9Fuo/UiArAggULlPwmRI2USiW+\n9rWvRb6fRIrKwYMHXdd99atfjTESIbJBoVDglltucV0/MDAQ2r4SKSoPPPCA6zoZt0JUj58x+9hj\nj4W2r0SKyv79+z1LH1x22WW0tbXFGJEQ6cXPmH3ggQdCLS2ZSFEpl8vs3bvX817Js2bNYuTIkTFG\nJUT6GD16tKcxu3btWt577718JL8dOXKEnTt38oc//MG1z9e//nUZt0K4UCqVmDdvnuv6jRs38ve/\n/z30G4wlVlTAEJauri7PU7O5c+fGGJEQ6SCIMfviiy9y+PDh0PedaFEpl8sMDg7y4IMPuvZpbW2V\ncStEBUGM2SgEBRIuKmAISxDjVhm3QhgENWbD9FHsJF5UIJhxq4xbIfyN2TVr1jgas2GWFkmFqIAh\nLH7G7YIFC1QmQeSWIMbsW2+95XiGkovZn0rK5bKMWyFc8DNm//nPf/Liiy9y5MiRyIY9FqkRFQhm\n3CrjVuQRP2P2iSee4PDhw66Ckvnhj9cfaBm3Dz30kGsfGbciT4wfP97TmH3wwQc5cOCA5xlK5oc/\n9j/QSWDK5TLvvvsu69atc92GjFuRB0aPHs1nP/tZ1/Vr167l3XffjXzIYyeRomKnXC47CsuRI0fo\n6enh97//vet7VSpBZJkgxmxvb2+sggIpEBVwPzU7fPgwXV1dnqUS5s6dq0r8InMUCgXP+/dYxmzc\nggIpERUvhoaGPEsltLa28slPfjLGiISIno6ODs8fSz9jNkpSLypBM25VKkFkhSAZs37GbJSkXlRg\nOOPWy7idNWsWra2tMUYlRPgEKWWwb9++hgkKZERUwBAWP+N2/vz5Mm5FaglayqCRggIZE5UjR46w\nadMm34xbGbcibQTJmP3jH/8Y2ZXH1ZAZUYHgpRJk3Iq04Zcx+/jjjzM0NBRjRO5kSlRApRJE9mh0\nKYNqyZyogEoliOwQpJRBo43ZSoKKShF4Gfil+foU4AVgO/A8cJKt7xJgB7ANuMLWfiGw1Vx3f+0h\nB8OqcauMW5FW/IzZDRs28NZbbx3X3mjPMKio3Aq8DlhyeAeGqEwCNpqvASYD15vPncDDgPUXPgLc\nDEw0H511xu7L4cOH2bRpEwcOHHDto1IJIokENWadShkEOWuJUniCiMqZwFXAaoYF4hpgjbm8BrjO\nXL4WWAcMAr1ADzAFGAeMATab/dba3nMcYf7Bg4ODnlc0t7a2ctFFF4W2PyHCoN5SBn5EOVwKIir3\nAbcD9jr+Y4F+c7nffA1wOrDL1m8XcIZDe5/Z7kiYf3CQUgmXXnqpjFuRGPxKGTQ6Y9aPks/6q4E9\nGH5Kh0ufMsPDolDo6uo6utze3k57e3td27OXSrjhhhsc+9x4442sXLnSc6gkRNSMGTPGs5RBmMZs\nb28vvb29dW+nEj9RuRhjqHMVMBI4EfgRxtlJG7AbY2izx+zfB5xle/+ZGGcofeayvb3PbacdHR1B\n4w+MvVTCJZdc4thnwYIF3HfffYlIIBL5o1QqeXp8XjVma8H+g10oFNi0aVMo2/Ub/tyJIRJnA58H\n/guYDTwDzDH7zAGeNpefMfu1mO+ZiOGj7AY+wPBXCuY2rPfEhlUqQRm3ImkELWUQxt0EC4XCccd4\nIyu/WXteDlyOMaU8zXwNxgzRz8zn54D5tvfMxzB7d2AYuOtrjroOhoaGlHErEkc9pQwskQj6Y2ht\nI6ofz2pEZRPGUAjgXWAGxpTyFcB7tn53AxOAc4Ff29r/DHzMXOc+VxYxyrgVSSOoMetGuVw++vDC\nLj5B+tdKJjNq/QiacatSCSJq/GrMWjf/gvqHKFHeQMxOLkUFgmXcqlSCiJJqbv5Vi6BUikblGYrO\nVCIgiHE7b948GbcidPyM2T179hyXMVvtcejmv0RNrkUF/I3bE044QTcnE6ETxJgdGho6ejcJ6ywj\nKNa2Kw3c3N5MLE6CGLfKuBVhEsSYtUoZWF92+9lKEAFw8k/s22v0tT+ZR6USRFxUY8zaBaWaKWML\n+9mN27NF5u9Q2AhUKkFETTXGrB0nMfETGPsZidsN+aIi16Ji/0cXCoWjNW79SiXIuBXVEsSYrSZj\n1u/Mwj7L4yQsGv5EhP2DsT6EIKUSZNyKapk6darnF/mHP/xh1decuZ3B2IdKlcJi0ejSB7nB+iCC\nlEo4/fTT4wpLpJzx48fz8Y9/3HW9PWM26JfdbXbH2obTD2YcJi1IVI5iTwyyl0pw4wtf+IIyboUv\nfqUM1q5de9SYhdq/8NUMbaKuwyJRwXnMaZVK+N3vfuf6PmXcCi+am5t9SxmEefOvuDwTPyQqJm7C\nIuNW1IJfjdlqjdmgVA59KhPflFEbE24ueblcZmhoyNe4VakEUUkUxmwQKnNbosxHcSPXolKp4vYP\nwm6ABSmVIONWWFRjzEZBpbDETa5FxQ0nhffLuJVxKyB4jdksk3lRqUapK42uysvEd+7cKeNWuOJX\nY7by5l9RDkVquQgxLDIvKl6Xf7sZWm79LONWpRJEJdWUMoDGDU3iIPOi4kQ1H6hdWCzjVqUSRCXV\n1JiN8+xBnkpEOGUfVl5K7vVBVxq6KpUg7FRTygCcj8MskQtRcbvmIegHW3lmoxq3wiKIMWsvZdCI\nvJG4z1YyJypu/8DK6x68alXY0/W95vtV4zbfBDVm46i25odKH9SIl6C43cYgaBUtpyFUkIxbGbfZ\npFpjNk9kSlTcZnoqx7JOV3UG2bbTGYtfqQQZt9nEL2PWMmb9iDPjtZbqcbWQKVGB488q3OpJ2KnG\nNKvMXwGVSsgbQTJmrbSDILM9fqkNaSNzogLOBpifKeZ2JuJUg6JShKxSCcq4zT5BM2YrK685UWvt\n2XrQtT8x4nXZeKWJa2+zJ8Yp4zbb1GLMeglGI85K3CyAMMm8qAT9J7rN9LiduTid+VjG7f79+133\no1IJ6SRIKYNKY9Z+LCUpH0VFmmqglmseKpPgvE5Z7esrE5mGhoY8E+NUKiGdTJ06laYm96+LUymD\nyqveG03lj2VUMWVOVKr5EP0KAVcKhtOvjpN7v3//fs9UfpVKSBfVljJwSnBLwpmK0wWzUQhLpkSl\n2g/PaSjjNhPkZPx6HTD79u2TcZsBaill4DYDGfd1P244TT6ESepFpV7jyekDr+baDLeZJpVKSD9+\nNWY3bNhwtMasUwa212xikgg7plSIiluKvVvfoP8kt0xZp/Ve+3ATHtW4TS9+GbP9/f2uGbNJGvL4\nYRfCsEi8qFSO/Spf1/rBOXkhTr8qTu1BruWwl0rQzcnSR0dHh+d6rxqzQQz/pFA5QREGiRcVpyuE\n7cJS64fm9M/0OitxisMpVqe2IKUSZNwmh2pLGYB3LlPShSVsgopKL9ANvAxsNttOAV4AtgPPAyfZ\n+i8BdgDbgCts7RcCW8119wfZsZeghEXYyUBOw6O9e/fy5JNPur5Hxm0yGD16dFWlDKznymPSzU+p\n16xNw5AqqKiUgQ7g3wEryeIODFGZBGw0XwNMBq43nzuBhwHrv/0IcDMw0Xx0eu3ULbPVyfsISq1X\nKbvR29vrur3Kg+5vf/tbIozbnTt3Rr6PMIkr3lKpxLx581zXu2XMOv3Yvfnmm67bqecHLIqznrCF\nqprhT+Vfcw2wxlxeA1xnLl8LrAMGMc5weoApwDhgDMNnOmtt7zl2R1X+46rpX82sTpD3VIpK5fvt\nsSWlVIJE5XhqNWbdZnvsx4VTykKWqeZMZQPwJ+ArZttYoN9c7jdfA5wO7LK9dxdwhkN7n9l+/M4i\ncs+DXpPhNesTdB9u8fsZtyqV0BiC3vzLKSvV63j1OoaTMgwKe/gfVFQuwRj6XAksAC6rWF82H6ER\ntrBUsz2nM5MgIhM01V+lEpKFX8bs/ffff4wxWzkN63Vc+aU/1EqYIpCEpLzvAIsxTNg2s22c+RoM\nb+UOW//1GMOfNuCvtvYbgP/tsP1XGBYpPfTQI75HFzHRiuGFAIwCfo8xo3Mv8C2z/Q5gubk8GUMY\nWoCzgZ0M+zEvYQhMAfgVPkatECKbnI0hEq8Ar2FMF4MxpbwB5ynlOzEM2m3A/7S1W1PKPcADkUYt\nhBBCCBEmnRhnNjsYHlY1gscxZrO22tpiSfSrg7OA3wB/wTibtKoJJTXukRhD4VeA14F7Eh6vnSJG\nEugvzddJj7mXBiWuNpoixpCoHWjGONjOa1Asl2HMdNlF5V7gm+bytzjeP2rGiL2HYf9oM8OJglH7\nR23Av5nLo4E3MP5/SY7bSh8uAX8ELk14vBaLgCeBZ8zXSY/5TQwRsZP0mEPhUxizRBaVM0hx086x\norKN4TycNoZnupZw7FnVeuAijNkw+0zX53Ge6YqKp4EZpCPuVmALcD7Jj/dMDB9xKsNnKkmP+U3g\n1Iq2SGNOygWFZwD/sL22EuaSQmSJfhHQjnGm9RLJjrsJ41exn+GhW5LjBbgPuB2wp9UmPeYyMSau\ngnHqmQTKjQ6gCqw5/SQyGngKuBX474p1SYv7CMaQ7UPArzF+/e0kLd6rgT0Y3kSHS5+kxQxG4uo7\nwIcxfJRtFetDjzkpZyp9GGajxVkcq4yNpp9jE/32mMuVcZ+JEXefuWxv74s4xmYMQfkRxvAH0hH3\n+8CzGEZgkuO9GON6tzcxrm2bhvG/TnLMYAgKwD+BX2D4IkmPORRKGEly7RhJc400auF4TyXpiX4F\njAs076toT2rcpzE843AC8FtgeoLjreQ/GPZUkhxz7hNXr8SYtehhOMGuEawD3gYOYfg8N5H8RL9L\nMYYTr2Ccnr+M8aEnNe6PAf/XjLcbw6cgwfFW8h8Mz/4kOWYlrgohhBBCCCGEEEIIIYQQQgghhBBC\nCCGEEEH4/ywwU8IaSQ5YAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10af50290>"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit -r 3 -n 1 model1d_vec(nx,nt,h,dt,v,delta,w)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 loops, best of 3: 273 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit -r 3 -n 1 model1d(nx,nt,h,dt,v,delta,w)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 loops, best of 3: 6min 26s per loop\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Numba"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numba\n",
      "model1d_nb=numba.autojit(model1d)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit -r 3 -n 1 model1d_nb(nx,nt,h,dt,v,delta,w)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 loops, best of 3: 6min 21s per loop\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Cython"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%load_ext cythonmagic"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%cython\n",
      "import cython\n",
      "import numpy as np\n",
      "cimport numpy as np\n",
      "@cython.boundscheck(False)\n",
      "@cython.wraparound(False)\n",
      "@cython.nonecheck(False)\n",
      "def model1d_cy(int nx,int nt,float h,float dt,np.ndarray[float] v,np.ndarray[float] delta,np.ndarray[float] w):\n",
      "    cdef np.ndarray[float] um=np.zeros_like(v)\n",
      "    cdef np.ndarray[float] uo=np.zeros_like(um)\n",
      "    cdef np.ndarray[float] ud=np.zeros_like(um)\n",
      "    cdef np.ndarray[float] up=np.zeros_like(um)\n",
      "\n",
      "    cdef float h2=h*h\n",
      "    cdef float dt2=dt*dt\n",
      "    cdef np.ndarray[float,ndim=2] seismo=np.zeros((nt,nx),dtype=np.float32)\n",
      "    cdef int it,ix\n",
      "    \n",
      "    for it in xrange(nt):\n",
      "        for ix in xrange(1,nx-1):\n",
      "            ud[ix]=(uo[ix-1]+uo[ix+1]-2.*uo[ix])/h2\n",
      "        for ix in xrange(nx):\n",
      "            up[ix]=(ud[ix]+delta[ix]*w[it])*v[ix]*v[ix]*dt2+2.*uo[ix]-um[ix]\n",
      "        um,uo,up=uo,up,um\n",
      "        for ix in xrange(nx):\n",
      "            seismo[it,ix]=uo[ix]\n",
      "    return seismo"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "seismo=model1d_cy(nx,nt,h,dt,v,delta,w)\n",
      "seismo=seismo.clip(seismo.min()*0.01,seismo.max()*0.01)\n",
      "plt.imshow(seismo,aspect=1,cmap=cm.gray_r)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 16,
       "text": [
        "<matplotlib.image.AxesImage at 0x1426ad690>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAARUAAAEACAYAAACd9eLKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH81JREFUeJztnXuQFOXd7z+zM7vgAsZbwuIlrFWAiubU+8ZUMF7yLhc9\nqxj1n8QYIZQxCQRSKhATkWAqf0TBVB1fL+gxoAZShiRVJpYpI1F4s+QkMUJy1MUYhCVuDKssWUR9\nT0FgF+b80d1sM/RtZrp7+vL9VE1Nz9PPdP92p+c7/XyfX/8ahBBCCCGEEEIIIYQQQgghhBBCiMTQ\nCWwDdgDfanAsQoiUUwR6gHagGXgFOK+RAQkhwqUp5v19EkNUeoFB4CfAtTHHIISIkLhF5QzgH7bX\nu8w2IURGiFtUyjHvTwgRM6WY99cHnGV7fRbG2cpRTj755PK+fftiDUoIAWPHjqW/v79Q73bq3kCV\nlIA3gOnA28Bm4Abgr7Y+5e985zsxh1U7XV1ddHR0NDqMqogi5tdee42nnnoKgJkzZ/KJT3witG3r\nfxwP3/3udyEETYh7+DMEfB34NfA68FOOFRSRQvbu3XtUUACeffZZdu/e3cCIRCOJW1QAngPOASYA\n9zRg/yJEhoaGeOihh45rf/TRRzl48GADIhKNphGikina29sbHULVhBnz9773Pdd1y5cvD2Ufef8f\npw2JSp2k8eAJK+af/vSnvn1Wr15d937y/D9OIxIVURNbt25l27Ztvv36+vr405/+FENEIilIVETV\nDAwM8POf/zxwfxm3+UKiIqpicHCQlStXVv0+Gbf5QaIiquLuu+92XXfRRRcxbdo01/VhGbci2UhU\nRGC8jNmmpiZmzpzJtGnTGD16tGu/xx57LIrQRIKQqIhA+BmzS5cuZeTIkYwYMYLFixe79tu1axdb\ntmyJIkSRECQqwhc/Y3bhwoWMGTPm6OtRo0bxjW98w7X/r371Kxm3GUaiIjwZGhryNGZnzZpFW1sb\nhcLwJSOFQoEPf/jD3HTTTa7vk3GbXSQqwhOvjNkpU6Zw/vnn09R0/GHU1NTEOeecI+M2h0hUhCt+\nxuzVV19NqeRePaNYLDJ9+nRP4zaMjFuRLCQqwpHXXnvN05j99re/zciRI32309LS4mncKuM2e0hU\nxHEMDAwcU8qgkoULF3qefVTiZ9wq4zZbSFTEMfhlzDoZs35Yxu2XvvQl1z4ybrODREUcg1fGrJcx\n60dTUxOTJk1i6tSprn1k3GYDiYo4ip8x+5nPfMbTmPWjWCwyY8YMz6HTqlWrat6+SAYSFQEEy5ht\naWmpez9+xu3bb78t4zblSFSEb8bsokWLGDNmTFU+ihejRo3i9ttvd10v4zbdSFRyThBjduzYsaEJ\nChjG7WmnnSbjNqNIVHKOXymDCy64oCZj1g/LuPXLuC2Xdf+5tCFRyTFBShkUi8XI9l8sFn1LJSjj\nNn1IVHJK0FIGUeNXKuHtt99WqYSUIVHJIXv37q2qlEHU+Bm3KpWQLiQqOePQoUOON/+yqCVjtl5k\n3GYLiUrOuOce95tCRmnM+iHjNjtIVHLET37yE9d1cRizflilEryGXjJuk49EJSds3bqVN954w3V9\n0FIGURMk41bGbbKRqOSAIMZsNaUMoqa1tVXGbYqRqGScwcHBxBmzfgQ1bg8dOhRjVCIoEpWME1Up\ng6gJYtzec889Mm4TSPKOJhEafsasX43ZRqMat+lEopJR/IzZuDJm60XGbfqQqGSQam/+lXR0c7J0\nIVHJGH6lDGbPnp04Y9YP1bhNFxKVjOFXymDy5MmJNGb9UMZteghydD0O9ANbbW2nAC8A24HngZNs\n65YAO4BtwBW29gvNbewA7q89ZOGGVymDYrGYeGPWD2XcpoMgovIE0FnRdgeGqEwCNpqvASYD15vP\nncDDgHWe/QhwMzDRfFRuU9RBd3e3bymDESNGxBhRNLS0tLBo0SLX9TJuG08QUfk/wL6KtmuANeby\nGuA6c/laYB0wCPQCPcAUYBwwBths9ltre4+ok7179/KLX/zCdX3SMmbrRaUSkk2tg+uxGEMizOex\n5vLpwC5bv13AGQ7tfWa7qBO/UgZpNGb9kHGbbMIYYJfNR2h0dXUdXW5vb6e9vT3MzWcKr1IGn/rU\npxKbMVsvhUKBSZMmMX36dDZu3OjYZ/ny5dx1112ZEtQw6e3tpbe3N/Tt1ioq/UAbsBtjaLPHbO8D\nzrL1OxPjDKXPXLa397ltvKOjo8aw8oVfxuxVV13V0FIGUVMsFpk6dSpbtmzhgw8+cOyzevVqvvKV\nr8QcWTqo/MHetGlTKNut9SfsGWCOuTwHeNrW/nmgBTgbw5DdjCE+H2D4KwVgtu09oga6u7tTUcog\nakaMGOFr3G7evNl1vQifIKKyDvgDcA7wD+AmYDlwOcaU8jTzNcDrwM/M5+eA+QwPjeYDqzGmlHuA\n9aH8BTkkb8asH36lEp577jkZtzESRFRuwDBaWzCGNk8A7wIzMKaUrwDes/W/G5gAnAv82tb+Z+Bj\n5rpb6g08r6SxlEHUqMZtssieg5dxvDJms2zM+mFl3E6fPt21jzJu4yF/R1+K8TJmC4UCV111Vaoz\nZuvFujmZMm4bi0QlJfgZs8uWLcuFMetHkIxbGbfRIlFJAQMDA57G7KJFi3JlzPrhl3Er4zZaJCoJ\nJ0gpg7Fjx+bKmPVDxm1jkagkHD9jNq2lDKJGpRIah47GBONlzJZKJWbOnJlrY9YPq1TCiSee6NpH\nxm34SFQSip8xe+edd2ailEHUtLS0sHDhQtf1Mm7DR6KSQPwyZmXMVoeM23iRqCSMIKUMZMxWh0ol\nxItEJUGUy+XcljKIGnupBDdk3IaDjs4E4VdjNuulDKLGKpUg4zZaJCoJwc+YTcvNv5LOiBEjfI3b\nl156KcaIsodEJQGolEG8+Bm369evl3FbBxKVBpPHGrONxjJub775Ztc+Mm5rR6LSQGTMNo5CocDE\niRNl3EaAjtYGkvcas43GKpUg4zZcJCoNoru7m+3bt7uuz0uN2UajjNvwkag0AL9SBjJm40UZt+Ei\nUYmZIKUMZMzGi0olhItEJUbK5bJKGSQUlUoIDx29MeJlzBaLRZUyaDBWqYQPfehDrn1WrVoVY0Tp\nRKISE37G7NKlS1XKIAG0tLRw2223ua5/5513ZNz6IFGJAZUySBcybutDohIxKmWQPoJm3P7rX/+K\nMar0IFGJEGXMppcgGbcrVqyQceuAjuYICWLMKmM2uQTJuJVxezwSlYiQMZsN/DJu33nnHZVKqECi\nEgEqZZAtVCqhOiQqIaNSBtlDpRKqQ6ISIjJms4tl3M6YMcO1jzJuDXR0h4hKGWQbq8atMm69kaiE\nhEoZ5IMgxm3eM24lKiGgjNl80draqoxbDyQqdaKM2fyhjFtvgojKWcBvgL8ArwG3mO2nAC8A24Hn\ngZNs71kC7AC2AVfY2i8Etprr7q8n8CQQxJhVKYNsEsS4zWvGbZCjfRBYCJwPXAQsAM4D7sAQlUnA\nRvM1wGTgevO5E3gYsH6mHwFuBiaaj84w/ohGsW7dOtd1KmWQfYIYtz/4wQ9ijCgZBBGV3cAr5vL/\nA/4KnAFcA6wx29cA15nL1wLrMMSoF+gBpgDjgDGA5WKttb0ndXR3d7Njxw7X9cqYzQd+xu3u3btz\nl3Fb7Xl5O/DvwEvAWKDfbO83XwOcDuyyvWcXhghVtveZ7alDxqyw42fc5i3jthpRGQ08BdwK/HfF\nurL5yDx+xuwXv/hFGbM5Q8btsQQd8DdjCMqPgKfNtn6gDWN4NA7YY7b3YZi7FmdinKH0mcv29j6n\nnXV1dR1dbm9vp729PWCY0eJnzF588cUyZnOK3bjdsGGDY58VK1awbNmyxBwfvb299Pb2hr7dIKJS\nAB4DXgf+09b+DDAHWGE+P21r/zHwvzCGNxMxfJQy8AGGv7IZmA084LTDjo6O6v6KmPAyZkulkjJm\nc45l3G7ZsoX333/fsc+qVauYO3duzJE5U/mDvWnTplC2G0QyLwFmAVOBl81HJ7AcuBxjSnma+RoM\n8fmZ+fwcMJ/hodF8YDXGlHIPsD6MPyIO/IzZO++8U8askHFLsDOV3+EuPm6T9Hebj0r+DHwswD4T\nhYxZUQ2Wcfv973/fcf369esZP348bW1tMUcWD8kY3CUYGbOiWvJu3EpUPJAxK2olaMbtkSNHYowq\nHvRt8ODHP/6x67pSqcSVV14pY1a4ktdSCRIVF1599VV6enpc18uYFUHIo3ErUXFgYGCAp59+2nW9\nVWNWPooIQmtrK9/85jdd12ct41aiUsGhQ4dYuXKl63rVmBXVUigUOPXUU3Nj3EpUbAQxZlVjVtRC\nU1NTboxbfTtseBmzhUJBGbOiLvJSKkGiYtLd3e1pzC5btkzGrKgbP+O2v78/9catRAVlzIp4ybpx\nm3tRUcasiJtCocBpp52WWeM216KijFnRKKyM2+nTp7v2Satxm+tvi1/GrIxZESXFYpFp06ZlzrjN\nraj4GbPKmBVxkEXjNpeiImNWJImsGbe5E5WDBw/KmBWJImvGba5EpVwus3z5ctf1l1xyiYxZ0RAs\n4/byyy937ZMW4zZX3x4vY7a5uZnOzk4Zs6JhFItFOjo6OOmkk1z7pMG4zY2o+JUyWLJkCSNHjowx\nIiGOp6Wlhdtuu811fRqM21yIil8pAxmzIkmk3bjNvKj4lTKQMSuSRtqN20yLil/GrIxZkVTSXOM2\n098mL2MWUI1ZkWisUglpM24zKyp+GbN33XWXMmZF4kmjcZtJUVHGrMgSaTNuMycqfhmzs2fPljEr\nUkXajNtMiUqQjFnVmBVpJE3Gbaa+XU8++aTrOt38S6Qdq1SCl3H76KOPxhiRM5kRlVdffZWdO3e6\nrlcpA5EFmpubPY3bPXv2NNy4zYSo+GXMLl68WMasyAxJN25TLyoHDx70zZj9yEc+ImNWZIagxu2B\nAwdijGqYVIuKShmIvBKkVMK9997bEOM21d82v1IGqjErskyQUgmNMG5TKypBShk0NzfHGJEQ8eOX\ncdsI4zaVohK0lIF8FJEHkmbcpk5UghizypgVeaJQKHDqqafy5S9/2bVPnMatn6iMBF4CXgFeB6w6\nAqcALwDbgecB+6BuCbAD2AZcYWu/ENhqrru/lmCDGLPnnXeejFmRO5qampgwYUIijFu/b9+/gKnA\nvwH/w1y+FLgDQ1QmARvN1wCTgevN507gYcA6ZXgEuBmYaD46qw3WK2O2WCxy5ZVXUiqVqt2sEJkg\nKcZtkJ/0/eZzC1AE9gHXAGvM9jXAdebytcA6YBDoBXqAKcA4YAyw2ey31vaeQHR3d3tmzC5dulQZ\nsyL3JMG4DSIqTRjDn37gN8BfgLHma8znseby6cAu23t3AWc4tPeZ7YEYGBhQKQMhAtJo4zaIqBzB\nGP6cCXwaYwhkp2w+IkHGrBDVYWXc+hm3UZVKqMaAeB94FsNw7QfagN0YQ5s9Zp8+4Czbe87EOEPp\nM5ft7X1uO+rq6gIMY/a3v/2ta0AXX3yxMmaFcKBQKDBhwgRmzJjBhg0bHPusWLGCT3/606H/IPt9\nG09jeGbnBOBy4GXgGWCO2T4HsJJGngE+j+G/nI1hyG7GEJ8PMPyVAjDb9p7j6OjooKOjg74+V92h\nVCopY1YID6xSCSeffLJrn23bth39voWFn6iMA/4Lw1N5CfglxmzPcgyB2Q5MM1+DMe38M/P5OWA+\nw0Oj+cBqjCnlHmC9146DlDJoaWnxCV+IfNPc3Mytt97quj4K4zaJRkR5wYIFnj7K4sWL5aMIEZBy\nuczAwAD33nuva5+5c+da0811f6mS+K30NH3nzJkjH0WIKimXy2zfvp3Vq1f7da1bE1L1zVTGrBC1\nYRm3Xhm3YZGab2dzczOdnZ0yZoWoESvj1su4DYPUiMqSJUuUMStEnbS0tHgat2GQClFRKQMhwsMv\n47ZeEi8qypgVIlyClEqoh0SKSrlcplwu8/7773PBBRfImBUiZJqamjj33HOPftfK5fCutEn0t/XE\nE09k/vz5jQ5DiMxRKBRYuHBhJNtOtKgAjBo1iilTpjQ6DCEyRUdHR2QjgMSLCsBll13GuHHjGh2G\nEJngox/9KBdeeGFk20+FqADceOONnHDCCY0OQ4hUM2rUKD73uc9Fuo/UiArAggULlPwmRI2USiW+\n9rWvRb6fRIrKwYMHXdd99atfjTESIbJBoVDglltucV0/MDAQ2r4SKSoPPPCA6zoZt0JUj58x+9hj\nj4W2r0SKyv79+z1LH1x22WW0tbXFGJEQ6cXPmH3ggQdCLS2ZSFEpl8vs3bvX817Js2bNYuTIkTFG\nJUT6GD16tKcxu3btWt577718JL8dOXKEnTt38oc//MG1z9e//nUZt0K4UCqVmDdvnuv6jRs38ve/\n/z30G4wlVlTAEJauri7PU7O5c+fGGJEQ6SCIMfviiy9y+PDh0PedaFEpl8sMDg7y4IMPuvZpbW2V\ncStEBUGM2SgEBRIuKmAISxDjVhm3QhgENWbD9FHsJF5UIJhxq4xbIfyN2TVr1jgas2GWFkmFqIAh\nLH7G7YIFC1QmQeSWIMbsW2+95XiGkovZn0rK5bKMWyFc8DNm//nPf/Liiy9y5MiRyIY9FqkRFQhm\n3CrjVuQRP2P2iSee4PDhw66Ckvnhj9cfaBm3Dz30kGsfGbciT4wfP97TmH3wwQc5cOCA5xlK5oc/\n9j/QSWDK5TLvvvsu69atc92GjFuRB0aPHs1nP/tZ1/Vr167l3XffjXzIYyeRomKnXC47CsuRI0fo\n6enh97//vet7VSpBZJkgxmxvb2+sggIpEBVwPzU7fPgwXV1dnqUS5s6dq0r8InMUCgXP+/dYxmzc\nggIpERUvhoaGPEsltLa28slPfjLGiISIno6ODs8fSz9jNkpSLypBM25VKkFkhSAZs37GbJSkXlRg\nOOPWy7idNWsWra2tMUYlRPgEKWWwb9++hgkKZERUwBAWP+N2/vz5Mm5FaglayqCRggIZE5UjR46w\nadMm34xbGbcibQTJmP3jH/8Y2ZXH1ZAZUYHgpRJk3Iq04Zcx+/jjjzM0NBRjRO5kSlRApRJE9mh0\nKYNqyZyogEoliOwQpJRBo43ZSoKKShF4Gfil+foU4AVgO/A8cJKt7xJgB7ANuMLWfiGw1Vx3f+0h\nB8OqcauMW5FW/IzZDRs28NZbbx3X3mjPMKio3Aq8DlhyeAeGqEwCNpqvASYD15vPncDDgPUXPgLc\nDEw0H511xu7L4cOH2bRpEwcOHHDto1IJIokENWadShkEOWuJUniCiMqZwFXAaoYF4hpgjbm8BrjO\nXL4WWAcMAr1ADzAFGAeMATab/dba3nMcYf7Bg4ODnlc0t7a2ctFFF4W2PyHCoN5SBn5EOVwKIir3\nAbcD9jr+Y4F+c7nffA1wOrDL1m8XcIZDe5/Z7kiYf3CQUgmXXnqpjFuRGPxKGTQ6Y9aPks/6q4E9\nGH5Kh0ufMsPDolDo6uo6utze3k57e3td27OXSrjhhhsc+9x4442sXLnSc6gkRNSMGTPGs5RBmMZs\nb28vvb29dW+nEj9RuRhjqHMVMBI4EfgRxtlJG7AbY2izx+zfB5xle/+ZGGcofeayvb3PbacdHR1B\n4w+MvVTCJZdc4thnwYIF3HfffYlIIBL5o1QqeXp8XjVma8H+g10oFNi0aVMo2/Ub/tyJIRJnA58H\n/guYDTwDzDH7zAGeNpefMfu1mO+ZiOGj7AY+wPBXCuY2rPfEhlUqQRm3ImkELWUQxt0EC4XCccd4\nIyu/WXteDlyOMaU8zXwNxgzRz8zn54D5tvfMxzB7d2AYuOtrjroOhoaGlHErEkc9pQwskQj6Y2ht\nI6ofz2pEZRPGUAjgXWAGxpTyFcB7tn53AxOAc4Ff29r/DHzMXOc+VxYxyrgVSSOoMetGuVw++vDC\nLj5B+tdKJjNq/QiacatSCSJq/GrMWjf/gvqHKFHeQMxOLkUFgmXcqlSCiJJqbv5Vi6BUikblGYrO\nVCIgiHE7b948GbcidPyM2T179hyXMVvtcejmv0RNrkUF/I3bE044QTcnE6ETxJgdGho6ejcJ6ywj\nKNa2Kw3c3N5MLE6CGLfKuBVhEsSYtUoZWF92+9lKEAFw8k/s22v0tT+ZR6USRFxUY8zaBaWaKWML\n+9mN27NF5u9Q2AhUKkFETTXGrB0nMfETGPsZidsN+aIi16Ji/0cXCoWjNW79SiXIuBXVEsSYrSZj\n1u/Mwj7L4yQsGv5EhP2DsT6EIKUSZNyKapk6darnF/mHP/xh1decuZ3B2IdKlcJi0ejSB7nB+iCC\nlEo4/fTT4wpLpJzx48fz8Y9/3HW9PWM26JfdbXbH2obTD2YcJi1IVI5iTwyyl0pw4wtf+IIyboUv\nfqUM1q5de9SYhdq/8NUMbaKuwyJRwXnMaZVK+N3vfuf6PmXcCi+am5t9SxmEefOvuDwTPyQqJm7C\nIuNW1IJfjdlqjdmgVA59KhPflFEbE24ueblcZmhoyNe4VakEUUkUxmwQKnNbosxHcSPXolKp4vYP\nwm6ABSmVIONWWFRjzEZBpbDETa5FxQ0nhffLuJVxKyB4jdksk3lRqUapK42uysvEd+7cKeNWuOJX\nY7by5l9RDkVquQgxLDIvKl6Xf7sZWm79LONWpRJEJdWUMoDGDU3iIPOi4kQ1H6hdWCzjVqUSRCXV\n1JiN8+xBnkpEOGUfVl5K7vVBVxq6KpUg7FRTygCcj8MskQtRcbvmIegHW3lmoxq3wiKIMWsvZdCI\nvJG4z1YyJypu/8DK6x68alXY0/W95vtV4zbfBDVm46i25odKH9SIl6C43cYgaBUtpyFUkIxbGbfZ\npFpjNk9kSlTcZnoqx7JOV3UG2bbTGYtfqQQZt9nEL2PWMmb9iDPjtZbqcbWQKVGB488q3OpJ2KnG\nNKvMXwGVSsgbQTJmrbSDILM9fqkNaSNzogLOBpifKeZ2JuJUg6JShKxSCcq4zT5BM2YrK685UWvt\n2XrQtT8x4nXZeKWJa2+zJ8Yp4zbb1GLMeglGI85K3CyAMMm8qAT9J7rN9LiduTid+VjG7f79+133\no1IJ6SRIKYNKY9Z+LCUpH0VFmmqglmseKpPgvE5Z7esrE5mGhoY8E+NUKiGdTJ06laYm96+LUymD\nyqveG03lj2VUMWVOVKr5EP0KAVcKhtOvjpN7v3//fs9UfpVKSBfVljJwSnBLwpmK0wWzUQhLpkSl\n2g/PaSjjNhPkZPx6HTD79u2TcZsBaill4DYDGfd1P244TT6ESepFpV7jyekDr+baDLeZJpVKSD9+\nNWY3bNhwtMasUwa212xikgg7plSIiluKvVvfoP8kt0xZp/Ve+3ATHtW4TS9+GbP9/f2uGbNJGvL4\nYRfCsEi8qFSO/Spf1/rBOXkhTr8qTu1BruWwl0rQzcnSR0dHh+d6rxqzQQz/pFA5QREGiRcVpyuE\n7cJS64fm9M/0OitxisMpVqe2IKUSZNwmh2pLGYB3LlPShSVsgopKL9ANvAxsNttOAV4AtgPPAyfZ\n+i8BdgDbgCts7RcCW8119wfZsZeghEXYyUBOw6O9e/fy5JNPur5Hxm0yGD16dFWlDKznymPSzU+p\n16xNw5AqqKiUgQ7g3wEryeIODFGZBGw0XwNMBq43nzuBhwHrv/0IcDMw0Xx0eu3ULbPVyfsISq1X\nKbvR29vrur3Kg+5vf/tbIozbnTt3Rr6PMIkr3lKpxLx581zXu2XMOv3Yvfnmm67bqecHLIqznrCF\nqprhT+Vfcw2wxlxeA1xnLl8LrAMGMc5weoApwDhgDMNnOmtt7zl2R1X+46rpX82sTpD3VIpK5fvt\nsSWlVIJE5XhqNWbdZnvsx4VTykKWqeZMZQPwJ+ArZttYoN9c7jdfA5wO7LK9dxdwhkN7n9l+/M4i\ncs+DXpPhNesTdB9u8fsZtyqV0BiC3vzLKSvV63j1OoaTMgwKe/gfVFQuwRj6XAksAC6rWF82H6ER\ntrBUsz2nM5MgIhM01V+lEpKFX8bs/ffff4wxWzkN63Vc+aU/1EqYIpCEpLzvAIsxTNg2s22c+RoM\nb+UOW//1GMOfNuCvtvYbgP/tsP1XGBYpPfTQI75HFzHRiuGFAIwCfo8xo3Mv8C2z/Q5gubk8GUMY\nWoCzgZ0M+zEvYQhMAfgVPkatECKbnI0hEq8Ar2FMF4MxpbwB5ynlOzEM2m3A/7S1W1PKPcADkUYt\nhBBCCBEmnRhnNjsYHlY1gscxZrO22tpiSfSrg7OA3wB/wTibtKoJJTXukRhD4VeA14F7Eh6vnSJG\nEugvzddJj7mXBiWuNpoixpCoHWjGONjOa1Asl2HMdNlF5V7gm+bytzjeP2rGiL2HYf9oM8OJglH7\nR23Av5nLo4E3MP5/SY7bSh8uAX8ELk14vBaLgCeBZ8zXSY/5TQwRsZP0mEPhUxizRBaVM0hx086x\norKN4TycNoZnupZw7FnVeuAijNkw+0zX53Ge6YqKp4EZpCPuVmALcD7Jj/dMDB9xKsNnKkmP+U3g\n1Iq2SGNOygWFZwD/sL22EuaSQmSJfhHQjnGm9RLJjrsJ41exn+GhW5LjBbgPuB2wp9UmPeYyMSau\ngnHqmQTKjQ6gCqw5/SQyGngKuBX474p1SYv7CMaQ7UPArzF+/e0kLd6rgT0Y3kSHS5+kxQxG4uo7\nwIcxfJRtFetDjzkpZyp9GGajxVkcq4yNpp9jE/32mMuVcZ+JEXefuWxv74s4xmYMQfkRxvAH0hH3\n+8CzGEZgkuO9GON6tzcxrm2bhvG/TnLMYAgKwD+BX2D4IkmPORRKGEly7RhJc400auF4TyXpiX4F\njAs076toT2rcpzE843AC8FtgeoLjreQ/GPZUkhxz7hNXr8SYtehhOMGuEawD3gYOYfg8N5H8RL9L\nMYYTr2Ccnr+M8aEnNe6PAf/XjLcbw6cgwfFW8h8Mz/4kOWYlrgohhBBCCCGEEEIIIYQQQgghhBBC\nCCGEEEH4/ywwU8IaSQ5YAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x143245a90>"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit -r 3 -n 1 model1d_cy(nx,nt,h,dt,v,delta,w)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 loops, best of 3: 517 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Fortran"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%file model1d_f.f90\n",
      "\n",
      "subroutine model1d_vec(nx,nt,h,dt,v,delta,w,seismo)\n",
      "implicit none\n",
      "integer,intent(in):: nx,nt\n",
      "real,intent(in):: h,dt,v(nx),delta(nx),w(nt)\n",
      "real,intent(out):: seismo(nx,nt)\n",
      "real,dimension(nx):: um,uo,up,ud\n",
      "real:: h2,dt2\n",
      "integer:: it\n",
      "\n",
      "um=0.\n",
      "uo=0.\n",
      "ud=0.\n",
      "seismo=0.\n",
      "\n",
      "h2=h*h\n",
      "dt2=dt*dt\n",
      "\n",
      "do it=1,nt\n",
      "    ud(2:nx-1)=(uo(1:nx-2)+uo(3:nx)-2.*uo(2:nx-1))/h2\n",
      "    up=(ud+delta*w(it))*v**2*dt2 + 2.*uo-um\n",
      "    um=uo\n",
      "    uo=up\n",
      "    seismo(:,it)=uo(:)\n",
      "enddo\n",
      "end subroutine\n",
      "                                                \n",
      "subroutine model1d(nx,nt,h,dt,v,delta,w,seismo)\n",
      "implicit none\n",
      "integer,intent(in):: nx,nt\n",
      "real,intent(in):: h,dt,v(nx),delta(nx),w(nt)\n",
      "real,intent(out):: seismo(nx,nt)\n",
      "real,dimension(nx):: um,uo,up,ud\n",
      "real:: h2,dt2\n",
      "integer it,ix\n",
      "um=0.\n",
      "uo=0.\n",
      "ud=0.\n",
      "seismo=0.\n",
      "\n",
      "h2=h*h\n",
      "dt2=dt*dt\n",
      "\n",
      "do it=1,nt\n",
      "    do ix=2,nx-1\n",
      "        ud(ix)=(uo(ix-1)+uo(ix+1)-2.*uo(ix))/h2\n",
      "    enddo\n",
      "    do ix=1,nx\n",
      "        up(ix)=(ud(ix)+delta(ix)*w(it))*v(ix)**2*dt2&\n",
      "        + 2.*uo(ix)-um(ix)\n",
      "    enddo\n",
      "    um=uo\n",
      "    uo=up\n",
      "    do ix=1,nx\n",
      "        seismo(ix,it)=uo(ix)\n",
      "    enddo\n",
      "enddo\n",
      "end subroutine"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Overwriting model1d_f.f90\n"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!f2py -c -m model1d_f model1d_f.f90 --f90exec=/opt/local/bin/gfortran-mp-4.9 > log.txt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import model1d_f\n",
      "print model1d_f.__doc__"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "This module 'model1d_f' is auto-generated with f2py (version:2).\n",
        "Functions:\n",
        "  seismo = model1d_vec(h,dt,v,delta,w,nx=len(v),nt=len(w))\n",
        "  seismo = model1d(h,dt,v,delta,w,nx=len(v),nt=len(w))\n",
        ".\n"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "vtr=v.transpose()\n",
      "dtr=delta.transpose()\n",
      "seismo=model1d_f.model1d(h,dt,vtr,dtr,w)\n",
      "seismo=seismo.clip(seismo.min()*0.01,seismo.max()*0.01)\n",
      "plt.imshow(seismo.transpose(),aspect=1,cmap=cm.gray_r)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 21,
       "text": [
        "<matplotlib.image.AxesImage at 0x143210c90>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAARUAAAEACAYAAACd9eLKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHuJJREFUeJztnX2QFOWdxz+zsyzIi1FjBHwJa5V4iubqLqaCUcwtL3oI\nRv0nMUYNdTEJBlIqEBOBYCp/RMFUnecLegZUIGVIUmVimTISgcuSM1EgOXUxBmFJ1oQ1LFl8uysI\n7MvcH93NNsN09zMz3T398v1UTU3P0890Pwvd3+n+Pr/5DgghhBBCCCGEEEIIIYQQQgghhBCJYSaw\nE9gNfKPBYxFCpJwi0Am0AsOAV4DzGzkgIUS4NMW8v49jiUoX0Af8ELgm5jEIISIkblE5A/iL6/Ve\nu00IkRHiFpVSzPsTQsRMc8z76wbOcr0+C+tq5Sgnn3xy6Z133ol1UEIIGDt2LD09PYV6t1P3Bqqk\nGXgDmA68BWwDrgf+4OpT+ta3vhXzsGqnvb2dtra2Rg+jKqIY82uvvcZTTz0FwOzZs/nYxz4W2rb1\nbxwP3/72tyEETYj79qcf+CrwC+B14EccKygihRw4cOCooAA8++yz7Nu3r4EjEo0kblEBeA74B+Ac\n4J4G7F+ESH9/Pw899NBx7Y8++iiHDx9uwIhEo2mEqGSK1tbWRg+hasIc83e+8x3PdcuXLw9lH3n/\nN04bEpU6SePBE9aYf/SjHwX2Wb16dd37yfO/cRqRqIia2LFjBzt37gzs193dzW9/+9sYRiSSgkRF\nVE1vby8/+clPjPvLuM0XEhVRFX19faxcubLq98m4zQ8SFVEVd999t+e6iy++mGnTpnmuD8u4FclG\noiKM8TNmm5qamD17NtOmTWP06NGe/R577LEohiYShERFGBFkzH7zm99kxIgRDB8+nEWLFnn227t3\nL9u3b49iiCIhSFREIEHG7IIFC465Ohk1ahRf+9rXPPv//Oc/l3GbYSQqwpf+/n5fY/bGG29k3Lhx\nFApDXxkpFAqcdtppfOELX/B8n4zb7CJREb74VcxefPHFXHjhhTQ1HX8YFQoFzj33XBm3OUSiIjwx\nMWaLxaJnn2KxyPTp032N2zAqbkWykKiIirz22mtGxmwQLS0tvv6KKm6zh0RFHEdvb+8xUQbllBuz\nQYwcOdJXWFRxmy0kKuIYgipmKxmzQRQKBT70oQ/JuM0JEhVxDH4Vs5MnT+aCCy6oaMwG0dTUJOM2\nJ0hUxFGCjNlPfepTNDfXHmtsYtyuWrWq5u2LZCBREUBwxezSpUsZPnx43ftpaWnxrbh96623ZNym\nHImKCKyYXbhwIWPGjAltf6NGjeKOO+7wXC/jNt1IVHKOiTE7duzYqozZIAqFAqeeeqqM24wiUck5\nQVEGXhWz9WJq3JZK+v25tCFRyTFBxuxVV13lWzFbL45x63drpYrb9CFRySkmUQZhGLNBmBi3ikpI\nFxKVHHLgwIGqogyiZuTIkb7GraIS0oVEJWccOXKk4o9/OdRSMVsvqrjNFhKVnHHPPd4/ChmlMRuE\naVSCjNvkI1HJET/84Q8915lEGUSNY9yeeOKJnn1k3CYfiUpO2LFjB2+88YbnetMog6hpaWlh4cKF\nnutl3CYfiUoOSJoxG4SM23QjUck4fX19iTNmgzCtuD1y5EiMoxKmSFQyTlRRBlFjUnF7zz33yLhN\nIMk7mkRoBBmzV111VV1RBlGjjNt0IlHJKEHG7NKlSxNhzAahitv0IVHJICY//hVmlEHU6MfJ0oVE\nJWMERRncdNNNiTNmg1DFbbqQqGSMoCiDSZMmJdKYDUJRCenB5Oh6HOgBdrjaTgE2AruA54GTXOsW\nA7uBncAVrvaL7G3sBu6vfcjCC78og2KxmHhjNghFJaQDE1F5AphZ1nYnlqicC2y2XwNMAq6zn2cC\nDwPOdfYjwM3ARPtRvk1RBx0dHbFkzDYaVdwmHxNR+W/gnbK2q4G19vJa4Fp7+RpgPdAHdAGdwGRg\nPDAG2Gb3W+d6j6iTAwcO8NOf/tRzfdIqZuslKONWxm1jqfXmeizWLRH281h7+XRgr6vfXuCMCu3d\ndruok6AogzQas0HIuE02Ydxgl+xHaLS3tx9dbm1tpbW1NczNZwq/KINPfOITia2YrRcnKmH69Ols\n3ry5Yp/ly5dz1113ZUpQw6Srq4uurq7Qt1urqPQA44B9WLc2++32buAsV78zsa5Quu1ld3u318bb\n2tpqHFa+CKqYnTVrVkOjDKKmWCwydepUtm/fzvvvv1+xz+rVq/nSl74U88jSQfkH9pYtW0LZbq0f\nYc8Ac+zlOcDTrvbPAi3A2ViG7DYs8Xkfy18pADe53iNqoKOjIxVRBlEzfPjwQON227ZtnutF+JiI\nynrgN8A/AH8B/g1YDlyONaU8zX4N8DrwY/v5OWAeQ7dG84DVWFPKncCGUP6CHJI3YzaIoKiE5557\nTsZtjJiIyvVYRmsL1q3NE8DbwAysKeUrgHdd/e8GzgHOA37hav8d8BF73a31DjyvpDHKIGpk3CaL\n7Dl4GcevYjbLxmwQbuPWC1XcxkP+jr4U42fMFgoFZs2aleqK2XopFotMmzZNFbcNRqKSEoKM2WXL\nluXCmA3CpOJWxm20SFRSQG9vr68xu3DhwlwZs0EEVdzKuI0WiUrCMYkyGDt2bK6M2SBMM25l3EaD\nRCXhBBmzaY0yiBpFJTQOHY0Jxs+YbW5uZvbs2bk2ZoPQj5M1BolKQgkyZpcsWZKJKIOoaWlpYcGC\nBZ7rZdyGj0QlgQRVzMqYrQ4Zt/EiUUkYJlEGMmarQxW38SJRSRClUim3UQZRo4rb+NDRmSCCMmZn\nz56d6SiDqHEqbmXcRotEJSEEGbNZyZhtNCYVt1u3bo1xRNlDopIAFGUQL0FRCRs2bJBxWwcSlQaT\nx4zZRuMYtzfffLNnHxm3tSNRaSAyZhtHoVBg4sSJMm4jQEdrA8l7xmyjkXEbDRKVBtHR0cGuXbs8\n1+clY7bRKCohfCQqDUBRBslCGbfhIlGJGUUZJA9V3IaLRCVGSqWSogwSilNxq6iE+tHRGyN+xqxT\nMasog8bhRCV84AMf8OyzatWqGEeUTiQqMRFkzKpiNhkERSX89a9/lXEbgEQlBhRlkC5k3NaHRCVi\nFGWQPkwrbv/+97/HOKr0IFGJEFXMpheTitsVK1bIuK2AjuYICTJmVTGbbIrFIlOnTvWtuJVxezwS\nlYgIMmaXLFmiitkUMHz48EDjVlEJxyJRiQCTKAO/n+YUySIo41ZRCcciUQkZRRlkD+fHyRSVYIZE\nJURkzGaXpqYmJk6cyIwZMzz7qOLWQkd3iCjKINs4xq0qbv2RqISEogzygSpug5GohIAqZvOFKm79\nkajUiSpm84epcZvXilsTUTkL+CXwe+A14Fa7/RRgI7ALeB44yfWexcBuYCdwhav9ImCHve7+egae\nBEyMWUUZZBPHuPWLSshrxa3J0d4HLAAuAC4G5gPnA3diicq5wGb7NcAk4Dr7eSbwMOB8TD8C3AxM\ntB8zw/gjGsX69es91ynKIPuYRCV873vfi3FEycBEVPYBr9jL/wf8ATgDuBpYa7evBa61l68B1mOJ\nURfQCUwGxgNjAMfFWud6T+ro6Ohg9+7dnusVZZAPgozbffv25a7ittrr8lbgn4GtwFigx27vsV8D\nnA7sdb1nL5YIlbd32+2pQ8ascKMfJzuWakRlNPAUcBvwv2XrSvYj8wQZs5///OdlzOYMRSUci+kN\n/zAsQfk+8LTd1gOMw7o9Gg/st9u7scxdhzOxrlC67WV3e3elnbW3tx9dbm1tpbW11XCY0RJkzF5y\nySUyZnOKE5UwY8YMNm3aVLHPihUrWLZsWWKOj66uLrq6ukLfromoFIDHgNeB/3C1PwPMAVbYz0+7\n2n8A/DvW7c1ELB+lBLyP5a9sA24CHqi0w7a2tur+ipjwM2abm5tVMZtznIrb7du3895771Xss2rV\nKubOnRvzyCpT/oG9ZcuWULZrIpmXAjcCU4GX7cdMYDlwOdaU8jT7NVji82P7+TlgHkO3RvOA1VhT\nyp3AhjD+iDgIMmaXLFkiY1bIuMXsSuUFvMXH69tVd9uPcn4HfMRgn4lCxqyoBse4/e53v1tx/YYN\nG5gwYQLjxo2LeWTxkIybuwQjY1ZUS96NW4mKDzJmRa24jVsvVqxYweDgYIyjigedDT784Ac/8FzX\n3NzMlVdeKWNWeJLXqASJigevvvoqnZ2dnutlzAoT8mjcSlQq0Nvby9NPP+253jFm5aMIE0aOHMnX\nv/51z/VZq7iVqJRx5MgRVq5c6bleUQaiWvIWlSBRcWFizCpjVtRCnoxbnR0u/IzZQqGgillRFybG\nbRaiEiQqNh0dHb7G7LJly2TMiroJMm57enpSb9xKVFDFrIiXrBu3uRcVVcyKuMm6cZtrUVHFrGgU\njnE7ffp0zz5pNW5zfbYEVczKmBVRUiwWmTZtWuaM29yKSpAxq4pZEQdZNG5zKSoyZkWSyJpxmztR\nOXz4sIxZkSiyZtzmSlRKpRLLly/3XH/ppZfKmBUNwTFuL7/8cs8+aTFuc3X2+Bmzw4YNY+bMmTJm\nRcMoFou0tbVx0kknefZJg3GbG1EJijJYvHixjFnRcFpaWrj99ts916fBuM2FqCjKQKSJtBu3mReV\noCgDGbMiaRQKBT74wQ+m1rjNtKgEVczKmBVJpampKbVRCZk+m/yMWUAZsyLROFEJaTNuMysqQRWz\nd911Fy0tLTGOSIjqSaNxm0lRMa2YlY8i0kDajNvMiUpQxawyZkXaSFvFbaZExaRiVhmzIo2kKeM2\nU2fXk08+6blOP/4l0o4TleBn3D766KMxjqgymRGVV199lT179niuX7JkiYxZkXqam5t9jdv9+/c3\n3LjNhKgEVcwuWrRIxqzIBIVCIfHGbepF5fDhw4EVs6eddpoERWQG04rbQ4cOxTiqIVItKooyEHnF\nqbj1i0q49957G2LcpvpsU5SByDMmUQmNMG5TKyqKMhDC+vBMmnGbSlFRlIEQFkk0blMnKibGrCpm\nRZ5ImnEbJCojgK3AK8DrgJMjcAqwEdgFPA+4b+oWA7uBncAVrvaLgB32uvtrGayJMXv++efLmBW5\nI0nGbdDZ93dgKvBPwD/ay1OAO7FE5Vxgs/0aYBJwnf08E3gYcC4ZHgFuBibaj5nVDtavYrZYLHLl\nlVfS3Nxc7WaFyARJMW5NPtIP2s8tQBF4B7gaWGu3rwWutZevAdYDfUAX0AlMBsYDY4Btdr91rvcY\n0dHR4Vsxu3TpUlXMitwzbNgwbrvtNs/1cRi3JqLShHX70wP8Evg9MNZ+jf081l4+Hdjreu9e4IwK\n7d12uxG9vb2KMhDCgEKhwKhRoxpq3JqIyiDW7c+ZwCexboHclOxHJMiYFaI6nKiEL37xi559ooxK\nqMaAeA94Fstw7QHGAfuwbm322326gbNc7zkT6wql2152t3d77ai9vR2wjNlf/epXngO65JJLVDEr\nRAUKhQLnnHMOM2bMYNOmTRX7rFixgk9+8pOhfyAHnY2nMjSzcwJwOfAy8Awwx26fAzhFI88An8Xy\nX87GMmS3YYnP+1j+SgG4yfWe42hra6OtrY3ubk/dobm5mVmzZqliVggPnKiEk08+2bPPzp07j55v\nYREkKuOB/8LyVLYCP8Oa7VmOJTC7gGn2a7CmnX9sPz8HzGPo1mgesBprSrkT2OC3Y0UZCFE/jTBu\nk2hElObPn+/royxatEg+ihCGlEolent7uffeez37zJ0715lurvukSuJZ6Wv6zpkzRz6KEFVSKpXY\ntWsXq1evDupatyak6sxUxawQteEYt34Vt2GRmrNTUQZC1IdTcetn3IZBakRFUQZC1E9LS4uvcRsG\nqRAVVcwKER5BUQn1knhRUcWsEOFiUnFbD4kUlVKpRKlU4r333uPCCy+UMStEyBQKBc4777yj51qp\nFN43bRJ9tp544onMmzev0cMQInMUCgUWLFgQybYTLSoAo0aNYvLkyY0ehhCZoq2tLbI7gMSLCsBl\nl13G+PHjGz0MITLBhz/8YS666KLItp8KUQG44YYbOOGEExo9DCFSzahRo/jMZz4T6T5SIyoA8+fP\nV/GbEDXS3NzMV77ylcj3k0hROXz4sOe6L3/5yzGORIhsUCgUuPXWWz3X9/b2hravRIrKAw884LlO\nxq0Q1RNkzD722GOh7SuRonLw4EHf6IPLLruMcePGxTgiIdJLkDH7wAMPhBotmUhRKZVKHDhwwPe3\nkm+88UZGjBgR46iESB+jR4/2NWbXrVvHO++8k4/it8HBQfbs2cNvfvMbzz5f/epXZdwK4UFzczO3\n3HKL5/rNmzfz5ptvhiookGBRKRQKDA4O0t7e7ntpNnfu3BhHJUQ6MDFmX3zxRQYGBmIPvm4YzvcR\n+vr6ePDBBz37jRw5UsatEGWYGLMDAwMA+bpSAesPNjFuVXErhIWpMRu2mDgkVlTcmBi3qrgVItiY\nXbt2bUVjNsxboMSKivNHO3+siXE7f/58xSSI3FIsFgON2T//+c8Vr1ByMfvj4P5jZdwKUZlCoeAb\nE/m3v/2NF198kcHBwcjHknhRcWNi3KriVuSRIGP2iSeeYGBgIDIfxU2qRAWGjNuHHnrIs4+MW5En\nJkyY4GvMPvjggxw6dCgWQYEUigpYwvL222+zfv16zz4ybkUeGD16NJ/+9Kc9169bt4633347NkGB\nlIoKWP5KZ2cnv/71rz37KCpBZBmTitmurq5YBQVSLCoAAwMDtLe3+0YlyLgVWcTUmI1bUCDlogLQ\n39/vG5WgiluRRdra2nxrS+I0ZstJvaiYVtwqKkFkBZOK2WqN2VwUv1WDU3HrZ9wqKkFkgUZEGVRL\nJkQFzIxbRSWINBOlMZuritpqGBgYYMuWLaq4FZkjKMqg1orZKH5OOFOiAigqQWSSoIrZxx9//GiU\ngSmFQoFSqZSfPJVaUVSCyBpRRRk4/XOTp1IPpVKJ3t5eRSWI1FNrlEEQUdz2OJiKShF4GfiZ/foU\nYCOwC3geOMnVdzGwG9gJXOFqvwjYYa+7v/Yhm1EqldizZ09gxa2iEkRSCTJmN23a5BllEESUs0Om\nZ9RtwOuAM5I7sUTlXGCz/RpgEnCd/TwTeBhwJPER4GZgov2YWefYA3GM20OHDnn28ftPE6JRmBiz\nL730Uk0+StSYiMqZwCxgNUMCcTWw1l5eC1xrL18DrAf6gC6gE5gMjAfGANvsfutc74mUvr4+3280\ny7gVSSQKYxaivUJxMBGV+4A7APdc1Vigx17usV8DnA7sdfXbC5xRob3bbo8cRSWItBEUZRB1xmy9\nNAesvwrYj+WntHn0KTF0WxQK7e3tR5dbW1tpbW2ta3vuqITrr7++Yp8bbriBlStX+t4qCRE1QVEG\ntRqzlejq6qKrq6vu7ZQTJCqXYN3qzAJGACcC38e6OhkH7MO6tdlv9+8GznK9/0ysK5Rue9nd3u21\n07a2NtPxG+NU3L7wwgtMmTKlYp/58+dz33331XRZKUS9mFTM1mrMVqL8A3vLli2hbDfo9mcJlkic\nDXwW+C/gJuAZYI7dZw7wtL38jN2vxX7PRCwfZR/wPpa/UrC34bwnNkwrbuMws4RwYxplkIYPvGrn\nUx2JXA5cjjWlPM1+DdYM0Y/t5+eAea73zMMye3djGbgbah51HZhU3H784x+PcURCmEcZpIFqRGUL\n1q0QwNvADKwp5SuAd1397gbOAc4DfuFq/x3wEXud91xZDKjiViQJE2M2zozZesll5Zdpxe3IkSNj\nHJXII3Eas3GRS1EBs4rbefPmKSpBRIaJMfvmm2+mSlAgx6ICQxm3Mm5F3AQZs/v374/tx7/CJtei\nAsq4FY0hS8ZsObkXlVKpxKFDh3yN2ylTpsi4FaFRTcVsGq+Scy8qoKgEER/VGLNOiFLakKjYmEYl\nyLgVtWJaMQscJyhpumKRqLgwiUqQcStqoVpjNq2CAhKV41BUgoiCqVOnGv/4V7mgpO0WSKJShklU\ngoxbUQ0TJkzgox/9qOf6SlEGjgCV58gWCoXEX7lIVCrgRCWo4lbUy5gxY2qqmC1/XS4mSRYWiYoH\ng4OD7NmzhxdeeMGzjypuhR/Dhg3z/Z0p04xZr1ugpAqLRMUHGbeiVoIyZvfv32+UMescW5WOsaR6\nLRKVAPr7+wONW0UliHKCjNk1a9b4CoqXp5KGWSGJSgCmP052+umnxzgqkWRMjNmgKINKnkql9UkU\nFomKASYVt5/73Odk3IqajVk/Kk0xJ9m0lagY4lTcyrgVXjQ3Nwcas7VEGbiN2kq3QUm7apGoVIFJ\nxu0tt9ySmP9cER8mFbMvvfRSXeaqc4VSaXo5ScecRKVK+vv7fTNuTzjhBFXc5hDTKINqRKXcrHUv\nl1feJgmJSpWYGLequM0XUfz4VyXvxMF9xRJUJNcIJCo1UCqVOHDggCpuhVGUwbvvvltzmb2XEFXy\nUZIy3SxRqRGn4lYZt/klKMrAXTFbfmURdNLXcmtTvs1GCUuqRaXRl3mDg4OBFbcybrOJqTHrzpit\n9OVA9/aqPU7kqYRMEr4SXiqVAqMSZNxmk2qjDNz+SCUDtlYqTTG78Svzj4rUikqjBcXBNCpBFbfZ\noZYoAzdukSn3RGo5+YPe52f6RkHqRCWJlYQmxq0qbrNBPRWzlSpjy0/2qD4s47yyT52olJte5Zd3\njRIaVdxmH5OK2aAoA79clKiP3UqFc1GQClHxm6t331M22mdxKm4PHjzo2cfvoBTJJSjKoKenxyjK\nIGjat5YT3n3VU8176tmnH4kXlUpXIl5fA0+Cz9Lf3+9bGKeM23QydepUmpq8T5egKINqqGUWqFYx\niuKqJZGi4vWHVhNY08jboIMHD/qW8isqIV0EGbP3339/1RWzboKiJCsR5kxS2B/GiRSVSs54Oe5/\n1Epi08irFpOMWxm36SDImF2zZs0xFbO1UMtv/PhNJTf6ij2RolLrVUaj/zHdyLhNP0HG7MaNG40y\nZivhVQjn9kaqOQ+8PlxN3xsmiRQV00Ss8rLnev5ho0AZt+mlUChw++23e67v6elh69atx1TMmuJ3\nlVHJQ6zlu0K1hkCFQSJFxeS2p9LrJH4dXD9Olk7a2tp8169Zs4b+/v6691MpJrLWArhKJRa1jKFe\nEikqUP8XqpIkLtVU3CZp3HnFNMogDMo/CGuZpAiDMLdrKipdQAfwMrDNbjsF2AjsAp4HTnL1Xwzs\nBnYCV7jaLwJ22Ovu99pZPWqdRJyK2yeffNKzj2Pc6laosZhEGVSbMWtCpRmgSjOaft/xScqxYyoq\nJaAN+GfA+T2KO7FE5Vxgs/0aYBJwnf08E3gYcP7aR4CbgYn2Y6bnDkP46nct76t2G11dXUb9SqUS\nf/zjHxNh3O7ZsyfyfYRJXOMNijLYuHGjccas6XHhxuuKJejLgl7baBTV3P6U/wVXA2vt5bXAtfby\nNcB6oA/rCqcTmAyMB8YwdKWzzvWe0KhUihyUM1HPl62qOXgGBgZob2/3rbiNIypBonI8QVEG1Rqz\ntYgKVBaWSsdypWM4zi8N+lHNlcom4LfAl+y2sUCPvdxjvwY4Hdjreu9e4IwK7d12u/dOq1Rdr/xO\nky9t1eKa18LAwICiEhKIyY9/hWHMBmFSGV5N2FMjMBWVS7Fufa4E5gOXla0v2Y/QMZ1eBv8iIpOv\nh8dBqVTi0KFDikpIEKYVs3Fg+qHmFfBUzwdxI/kWsAjLhB1nt423X4Plrdzp6r8B6/ZnHPAHV/v1\nwH9W2P4rDImUHnroEd+jnZgYieWFAIwCfo01o3Mv8A27/U5gub08CUsYWoCzgT0M+TFbsQSmAPwc\nH6NWCJFdzsYSiVeA17Cmi8GaUt5E5SnlJVgG7U7gX13tzpRyJ/BApKMWQgghhAiTmVhXNrsZuq1q\nBI9jzWbtcLVFVugXEmcBvwR+j3U16aQJJXXcI7BuhV8BXgfuSfh43RSxikB/Zr9O+pi7iLFwNUkU\nsW6JWoFhWAfb+Q0ay2VYM11uUbkX+Lq9/A2O94+GYY29kyH/aBtDhYJR+0fjgH+yl0cDb2D9+yV5\n3E7uQzPwEjAl4eN1WAg8CTxjv076mP+EJSJukj7mUPgE1iyRQ/kMUty0cqyo7GSoDmccQzNdizn2\nqmoDcDHWbJh7puuzVJ7pioqngRmkY9wjge3ABSR/vGdi+YhTGbpSSfqY/wR8sKwt0jEn5QuFZwB/\ncb12CuaSQuSFfiHSinWltZVkj7sJ61Oxh6FbtySPF+A+4A7AXVab9DGXiLlwtbnuIYdDqdEDqAJn\nTj+JjAaeAm4D/rdsXdLGPYh1y/YB4BdYn/5ukjbeq4D9WN5Em0efpI0ZrMLVvwIfwvJRdpatD33M\nSblS6cYyGx3O4lhlbDQ9HFvot99eLh/3mVjj7raX3e3dEY9xGJagfB/r9gfSMe73gGexjMAkj/cS\nrO+7/Qnru23TsP6tkzxmsAQF4G/AT7F8kaSPORSasYrkWrGK5hpp1MLxnkrSC/0KWF/QvK+sPanj\nPpWhGYcTgF8B0xM83nL+hSFPJcljzn3h6pVYsxadDBXYNYL1wFvAESyf599IfqHfFKzbiVewLs9f\nxvpPT+q4PwL8jz3eDiyfggSPt5x/YWj2J8ljVuGqEEIIIYQQQgghhBBCCCGEEEIIIYQQQghhwv8D\nDT5CJjWErhoAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x142a24050>"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit -r 3 -n 1 model1d_f.model1d(h,dt,vtr,dtr,w)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 loops, best of 3: 99.9 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit -r 3 -n 1 model1d_f.model1d_vec(h,dt,vtr,dtr,w)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 loops, best of 3: 102 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 23
    }
   ],
   "metadata": {}
  }
 ]
}