{
 "metadata": {
  "name": "",
  "signature": "sha256:bce36d54aa16d071bdcf1f06f6466931e0627af4168223050848a96cd01ef6a1"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "from matplotlib import animation\n",
      "from fatiando import utils\n",
      "from fatiando.seismic import wavefd\n",
      "from fatiando.vis import mpl"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Set the parameters of the finite difference grid\n",
      "area = [0, 20000, 0, 5000]\n",
      "velocity = utils.fromimage('scipy2014.png', ranges=[3500, 4500])\n",
      "density = utils.fromimage('scipy2014.png', ranges=[2200, 2500])\n",
      "shape = velocity.shape\n",
      "mu = wavefd.lame_mu(velocity, density)\n",
      "\n",
      "# Make a wave source from a mexican hat wavelet\n",
      "sources = [wavefd.MexHatSource(1000, 0, area, shape, 100000000, 10, delay=0.1),\n",
      "           wavefd.MexHatSource(19000, 0, area, shape, 100000000, 10, delay=0.1)]\n",
      "\n",
      "# Get the iterator for the simulation\n",
      "dt = wavefd.maxdt(area, shape, velocity.max())\n",
      "duration = 7\n",
      "maxit = int(duration/dt)\n",
      "snapshot = int(0.1/dt) # Plot a snapshot of the simulation every 0.5 seconds\n",
      "simulation = wavefd.elastic_sh(mu, density, area, dt, maxit, sources, None,\n",
      "                               snapshot, padding=100, taper=0.007)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# This part makes an animation using matplotlibs animation API\n",
      "fig = mpl.figure(figsize=(6, 1.5))\n",
      "ax = mpl.subplot(1, 1, 1)\n",
      "wavefield_plt = mpl.imshow(np.zeros(shape), vmin=-10**(-5), vmax=10**(-5), cmap=mpl.cm.seismic)\n",
      "mpl.imshow(np.ma.masked_array(velocity, velocity < velocity.max() - 10), cmap=mpl.cm.gray_r, alpha=0.9)\n",
      "ax.axison = False\n",
      "mpl.subplots_adjust(top=1, bottom=0, right=1, left=0)\n",
      "# Update the plot everytime the simulation yields\n",
      "def animate(i):\n",
      "    # Grab the iteration number, displacment panel and seismograms\n",
      "    t, u, seismograms = simulation.next()\n",
      "    wavefield_plt.set_array(u) # Revert the z axis so that 0 is top\n",
      "    mpl.savefig('snapshots/frame%05d.png' % (i + 1), dpi=50)\n",
      "    return wavefield_plt\n",
      "anim = animation.FuncAnimation(fig, animate, frames=maxit/snapshot, interval=1)\n",
      "mpl.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "! convert -coalesce -layers Optimize -delay 10 -size 300x75 -loop 0 snapshots/*.png scipy2014.gif"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.display import Image\n",
      "Image(url='files/scipy2014.gif')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<img src=\"files/scipy2014.gif\"/>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "<IPython.core.display.Image at 0x7fb71edb1250>"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}