rtm-hw2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import segyio\n",
    "import matplotlib.pylab as plt\n",
    "from examples.seismic import plot_shotrecord\n",
    "%matplotlib inline\n",
    "from examples.seismic import plot_velocity, plot_perturbation\n",
    "from scipy import ndimage\n",
    "from examples.seismic import Model\n",
    "names = locals()\n",
    "\n",
    "shape = (201, 101)   # Number of grid point (nx, nz)\n",
    "spacing = (5, 5)  # Grid spacing in m.\n",
    "origin = (0., 0.)  # What is the location of the top left corner.\n",
    "\n",
    "v = np.ones(shape, dtype = np.float32)*2.0\n",
    "'''\n",
    "for i in range(shape[0]):\n",
    "    v[i, :20] = 1.5\n",
    "    v[i, 20:80] = 1.7\n",
    "    v[i, 80:] = 1.9\n",
    "'''\n",
    "\n",
    "model = Model(vp=v, origin=origin, shape=shape, spacing=spacing,mytime=0.3,\n",
    "                      space_order=4)\n",
    "\n",
    "nshots = 51  # Need good covergae in shots, one every two grid points\n",
    "nreceivers = 204  # One recevier every grid point\n",
    "fwi_iterations=5 #number of outer FWI iterations\n",
    "t0 = 0.\n",
    "tn = 600.  # Simulation last 3.5 second (600 ms)\n",
    "f0 = 0.015  # Source peak frequency is 25Hz (0.015 kHz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAADPCAYAAABRPx4bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAd0klEQVR4nO3dfbxmc73/8dd7xuReaFRuG4XklKKheugcCR06B904RRGiSXFKOvWT/NSDyk3F8SulIXclKkpTTYkSqci4HWNSc6hM9EOEfpgxM+/fH2tdulyufe21b9be19r7/Xw8rodrrfVda33W3vPYH+u7vuvzlW0iIiL6yZTxDiAiIqJTklNERPSdJKeIiOg7SU4REdF3kpwiIqLvJDlFRETfqTU5SdpN0h2SFkk6qsv2lSV9s9x+naQZdcYTERHdSdpY0pWSFkpaIOmDXdq8U9Kt5edXkl5eVzy1JSdJU4HTgd2BrYB9JW3V0exg4CHbmwGnAifVFU9ERPS0DPiw7ZcArwYO6/I3+y5gR9tbA8cDs+sKps47p+2BRbbvtL0UuAjYq6PNXsB55feLgZ0lqcaYIiKiC9v32r6x/P4osBDYsKPNr2w/VC5eC2xUVzx1JqcNgbvblhfTcaHtbWwvAx4GnlNjTBERMYjyEcs2wHU9mh0M/KiuGFaq68BAtzugzlpJVdogaRYwC2Da6tNeOX3L5K+ImPjuveEvD9heD2AzyY+N5FiwAHiibdVs28/olpO0BnAJcITtR7odS9JOFMnptSMIqac6k9NiYOO25Y2AewZos1jSSsCzgQc7D1T+AGcDbDBzfR8876BaAo6I6Cef0gl/bH1/HDhsBMc6Bp6wPbNXG0nTKBLTBba/M0CbrYGzgN1t/3UEIfVUZ3K6Hthc0qbAn4F9gHd0tJkDHAD8Gtgb+JlTiTYi4hmmAKvWePzyef9XgYW2TxmgzSbAd4D9bf+uxnDqS062l0k6HLgMmAqcbXuBpOOAebbnUPwgviZpEcUd0z51xRMR0WQCptV7ih2A/YH5km4u1x0NbAJg+wzgWIpxAV8qx64tG+xubLjqvHPC9lxgbse6Y9u+PwH8R50xRERMBHXfOdm+hu7jANrbHAIcUmMYT6k1OUVExOgQk+sP9mS61oiIxqr7zqnfJDlFRDTAFGC18Q5iDCU5RUQ0xGT6gz2ZrjUiorHSrRcREX0nySkiIvrOGLzn1FeSnCIiGiB3ThER0XfynlNERPSdDCWPiIi+I2DVkfzFXjZakYyNJKeIiAaQYKUkp4iI6CdTpsCqK4/gAE8M3qSfJDlFRDSABNNGkpwaJskpIqIJJtlwvUl0qRERDSZglfEOYuwkOUVENIEo5hSfJJKcIiKaYAq5c4qIiD4jIAMiIiKir2RARERE9J3cOUVERN/JM6eIiOhLk2i03pQ6Dy5pN0l3SFok6agu2/9F0o2Slknau85YIiIarXXnNNxPw9SWnCRNBU4Hdge2AvaVtFVHsz8BBwLfqCuOiIgJofWe03A/DVNnt972wCLbdwJIugjYC7i91cD2H8ptK2qMIyKi+VIhYtRsCNzdtrwYeFWN54uImLimkNF6o0Rd1nlYB5JmAbMA1tpkrZHEFBHRTHnPadQsBjZuW94IuGc4B7I9G5gNsMHM9YeV4CIiGi3vOY2a64HNJW0K/BnYB3hHjeeLiJi4Jtl7TrWN1rO9DDgcuAxYCHzL9gJJx0naE0DSdpIWA/8BfEXSgrriiYhovIzWGx225wJzO9Yd2/b9eoruvoiI6KXmOydJGwPnA88HVgCzbZ/W0WZL4BxgW+Djtj9XVzyT6PFaRESD1T8gYhnwYds3SloTuEHS5bZvb2vzIPAB4E21RkKSU0REM9Q8IML2vcC95fdHJS2keCWo/d3U+4D7JP1bfZEUkpwiIppgDEfrSZoBbANcNzZnfKYkp4iIJhh5t950SfPalmeXr+k8/TTSGsAlwBG2HxnRGUcgySkioglGXiHiAdszezWQNI0iMV1g+zsjOtsIJTlFRDRBzbX1JAn4KrDQ9imjeNzVgSdsLx/KfklOERFNUP9ovR2A/YH5km4u1x0NbAJg+wxJzwfmAWsBKyQdAWzV3v0naQpF0YV3AtsBS4CVJd1P8WrRbNu/HyyYJKeIiCaoufCr7WvoXhO1vc1fGPzd1CuBK4CPAbfZXgEgaV1gJ+BESd+1/fVeB0lyiohoimb8xd7F9pOdK20/SPE865Ly2VZPzbjUiIjJriFTZrQSk6QXAYttL5H0OmBr4Hzbf+uWvDrVOk17RESMktZ7TsP9jL1LgOWSNqMYaLEpQ5j1PHdOERFN0Lz5nFbYXibpzcB/2/6CpJuq7tysS42ImKyaN2XGk5L2BQ4A9ijXDfqsqSXJKSKiCZo32eBBwKHAp23fVc7t13OEXrskp4iIJmhIt56k2cCPgCtsf6C13vZdwIlVj9PzUiVtRPEy1T8DGwCPA7cBPwR+1Bq/HhERNWvOndPZwG7AkZKWAj8Bfmz7lqEcZMDkJOkcinLpPwBOAu6j6PHcojzxxyUdZfvq4cUfERGVNeTOyfa1wLXAJyU9B3gD8GFJWwM3UiSqbw12nF6X+nnbt3VZfxvwHUnPoixrERERNRO4GXdOT7H9V+DC8oOkV1Lc3AxqwOQ0QGJq374UWFQ9zIiIGC5PgaUNGq0naW3gXcAM2nJN+3OoXga9SZT078DxwAvK9iqO77WGEW9ERAyDBcumjqRuwpgPEZhL0b03fzgnr9KD+d/AW4D5tj3UE0RExMhZYunKI+nXe3zUYqloFdtHDnfnKsnpborKsklMERHjZAVTWMKzRnCEMU9OX5P0HopBdUtaK8sCsIOqkpw+CsyVdFXHCQadjErSbsBpwFTgLNsndmw/EjgEWAbcD7zb9h+rBB4RMdksb8JwvX9YCnwW+DjQurkx8MIqO1e50k8Df6cYRl45bUuaCpwO7AosBq6XNMf27W3NbgJm2n5M0vuAk4G3Vz1HRMRksYIpLB3RndOYOxLYzPYDw9m5SnJa1/YbhnHs7YFFtu8EkHQRsBfwVHKyfWVb+2uB/YZxnoiISWE5U8c7hKFYADw23J2rJKcrJL3B9k+GeOwNKZ5XtSwGXtWj/cEUJS8iIqJD8cypUS86LQdulnQlT38kNDpDyYHDgI9KWgI8SfWh5N2m++06qELSfsBMYMcBts8CZgGstUlGsEfE5GPUtG69S8vPsAyanGyv2blOUs955kuLgY3bljcC7ulyrF0oHpjtaHtJ5/YyhtnAbIANZq6fUYMRMekYWNasbr3bbN/QvkLSHgM17jToG12SjutYnkK1sufXA5tL2rQsdbQPMKfjWNsAXwH2tH1f1aAjIiYbM4WlrDzszzg4U9LLWgvl3E7HVN25yuvGm0j6WHnwlSlu034/2E62lwGHA5cBC4Fv2V4g6ThJe5bNPgusAXxb0s2S5gxwuIiISW85U4f9GQd7A+dJekn5vtP7KYrAVlLlmdNBwAVlgtqJYqqMU6sc3PZcihIW7euObfu+S9VAIyIms5G/hDu2bN8paR+KG5q7gTfYrvwmcK8pM7ZtWzyNovvtl8BVkra1feMwY46IiCEqBkT0/2g9SfN5+uC3dSkKMVwnCdtbVzlOzykzOpYfArYq1xt4ffVwIyJiJIya8p7Tv4/GQXpNmbHTaJwgIiJGzqgp7zn91fbfezWQtMZgbQYcECFpv15DxiW9SNJrB48zIiJGqvWe03A/Y+h7kj4v6V8krd5aKemFkg6WdBkVJhzs1a33HIq3e28AbqAozLoKsBnFy7IPAEeN5AoiIqIao0a852R7Z0lvBN4L7CBpHYri3ncAPwQOsP2XwY7Tq1vvNElfpHi2tAOwNUXN9YXA/rb/NPLLiIiIKpoyIAK6j9Qeqp5DyW0vBy4vPxERMU4aNCBiVDRqcpCIiMmqGBBR37MjSRsD5wPPp5hWfbbt0zraiOLVojdSVBw/sK7XipKcIiIaYEVZvqhGy4AP275R0prADZIu75iDb3dg8/LzKuDL9J5tYtiqlC+KiIg+UGf5Itv3tu6CbD9KMb5gw45mewHnu3AtsLak9bsdT9LnJP3TcK910Dunsp7eW4EZ7e1tHzfQPhERMbrGcj4nSTOAbYDrOjZ1m6dvQ+DeLof5LTBb0krAOcCFth+uGkOVbr3vAQ9TDCfvOqVFRETUaxTmc5ouaV7b8uxyOqKnkbQGcAlwhO1HOjd3Da3bSvss4CxJL6ao0XqrpF8CZ3bMgt5VleS0ke1BX5iKiIh6jXC03gO2Z/ZqIGkaRWK6wPZ3ujSpNE9f2/GmAluWnweAW4AjJb3X9j69YqmSnH4l6WW251doGxERNai7W68cifdVYKHtUwZoNgc4XNJFFAMhHrbdrUsPSacAewA/Az5j+zflppMk3TFYPL2qkrcqy64EHCTpTopuvdY07ZUqy0ZExOio+T2nHYD9gfmSbi7XHQ1sAmD7DIoXa98ILKIYSn5Qj+PdBhxj+7Eu27YfLJhed06jUlk2IiJGru75nGxfQ/dnSu1tDBxW8ZDvtH12+wpJP7W9c5WBEb3KF/2xPNjXbO/fcYKvUWTYiIgYA00pXyRpFWA1igEY6/CPhLcWsEHV41R55vS0cerlA65XVj1BRESMXIPKF70XOIIiEbVXj3gEOL3qQXo9c/oYRX/jqpIe4R/ZbynwjOGHERFRn7rLF42WsuTRaZL+0/YXhnucXt16JwAnSDrB9seGe4KIiBi5BnXrvd72z4A/S3pL5/YBhqg/Q5VuvaPLE7yWYvTeL2xfOqRoIyJiRBrUrbcjxfDxPbpsMzBqyel0igkGLyyXD5W0q+2qIzYiImKEmjJNu+1PlP/tNcx8UFWS047AS8shhEg6D8gLuRERY6hBd04ASPoMcLLtv5XL61BUPT+myv5VqpLfQfkSVmlj4NaKwe0m6Q5JiyQ9Y0p3SYdKmi/pZknXSNqqynEjIiab1oCI4X7Gwe6txARg+yGKF3grqXLn9BxgoaRW6YntgF9LmlOecM9uO5VDzk8HdqWox3S9pDkdc4N8o3zrGEl7AqcAqeMXEdFhDOZzGm1TJa1sewmApFWh+gVUSU7HDjOw7YFFtu8sA7uIYi6Qp5JTR8Xb1Rmgum1ERNRevmi0fR34qaRzKP62vxs4r+rOgyYn21dJegGwue0ryuy3UjkZVS/d5v14xoyJkg4DjgSeBby+auAREZPJioa859Ri+2RJtwK7lKuOt31Z1f2rTDb4HmAWsC7wIooS6WcAOw+2a7d4n7HCPh04XdI7gGOAA7rEMKuMgbU2WWuwkCMiJhw3r1sP4CZgGsXf/puGsmOVARGHUVSrfQTA9u+B51bYb0jzfgAXAW/qtsH2bNszbc9cfb3VKpw6ImLiqXOa9tEm6W3Ab4C9gbcB10nau+r+VZ45LbG9tJjqA8opd6s8G7oe2FzSpsCfgX2Ad3QEv3mZ7AD+Dfg9ERHxDGM5Tfso+Tiwne37ACStB1wBXFxl5yrJ6SpJrRp7uwLvB74/2E62l0k6HLgMmAqcbXuBpOOAebZbk1btAjwJPESXLr2IiCg0bEDElFZiKv2Var11QLXkdBRwMMWLt++lmGzqrCoHtz23bN++7ti27x+sGmhExGS2giksXdGcARHAjyVdxj+qC72djnzQS5XReiskXQpcavv+4cUYEREj4RViyRPN6daz/RFJb6UYsyBgtu3vVt2/15QZAj4BHF4eWJKWA1+wfdzIwo6IiKGwxfJljerWw/YlwCXD2bfXndMRFBlvO9t3AUh6IfBlSR+yfepwThgREUPnFWLJ4/3frSfpUboPmhPFTO+V3gfqlZzeBexq+4HWCtt3StoP+AmQ5BQRMVYsVizp/24922uOxnF6Jadp7Ymp7cT3S5o2GiePiIiKLGhYt56k11JUFzpH0nRgzVZP3GB6Jaelw9wWERGjzcAT3Qrv9CdJnwBmAi8GzqEoUfd1isdFg+qVnF4u6ZEu6wWsMsQ4IyJiJAwsG+8ghuTNwDbAjQC275FUuctvwORku1n3jxERE9kK4PHxDmJIltq2pNZEtasPZecqL+FGRMR4M7BkvIMYkm9J+gqwdllA/N3AmVV3TnKKiGiChnXr2f5cWfLuEYrnTsfavrzq/klOERFNsAJ4YryDGJykL1LMcv6rMhlVTkjtkpwiIprANCI5Ucwu8XlJ6wPfBC60ffNQD1K5QmxERIwjU8zfMNzPWIVpn2b7NcCOwIPAOZIWSjpW0hZVj5PkFBHRBCsoBkQM9zPGbP/R9km2t6GYy+/NwMKq+yc5RUQ0QWtAxHA/g5B0tqT7JN02wPZ1JH1X0q2SfiPppYMcb5qkPSRdAPwI+B3w1iqXCnnmFBHRDPU/czoX+CJw/gDbjwZutv1mSVsCpwM7dzYqR+jtSzG7+W+Ai4BZtv/fUIJJcoqIaIKaR+vZvlrSjB5NtgJOKNv+VtIMSc+z/X872h0NfAP4L9sPDjeeJKeIiCZoDYgYP7cAbwGukbQ98AJgI+Bpycn2TqNxsiSniIgmGHmFiOmS5rUtz7Y9ewj7nwicJulmYD5wEzW+FpzkFBHRBCPv1nvA9szh7mz7EeAgeGqm9LvKTy2SnCIimmCcu/UkrQ08ZnspcAhwdZmwapHkFBHRBDUXfpV0IfA6iu6/xcAngGkAts8AXgKcL2k5cDtwcH3R1JycJO0GnAZMBc6yfeIA7fYGvg1sZ3tetzYREZNazYVfbe87yPZfA5vXF8HT1ZacJE2lGAe/K7AYuF7SHNu3d7RbE/gAcF1dsURENF5DCr+OljorRGwPLLJ9Z9lHeRGwV5d2xwMnM6l+7BERQ9R6CXe4n4apMzltCNzdtry4XPcUSdsAG9v+QY1xREQ0X0MKv46WOp85qcs6P7VRmgKcChw46IGkWcAsgLU2WWuUwouIaJBW4ddJos7ktBjYuG15I+CetuU1gZcCPy+GzPN8YI6kPTsHRZQvis0G2GDm+iYiYrJZATw+3kGMnTqT0/XA5pI2Bf4M7ENRNh0A2w8D01vLkn5OUYspo/UiIrpZPt4BjJ3akpPtZZIOBy6jGEp+tu0Fko4D5tmeU9e5IyImnEk2Wq/W95xszwXmdqw7doC2r6szloiIRhv/wq9jKhUiIiKaoOYKEf0mySkiogkyICIiIvqOyYCIiIjoM/VP095XkpwiIpogAyIiIqLvpEJERET0nQyIiIiIvpQBERER0XcmUWXROqfMiIiIGJbcOUVENMLkGq6X5BQR0QiTa0REklNERCMYWDbeQYyZJKeIiEZYATw23kGMmSSniIhGMOnWi4iIPpNuvYiI6DsZEBEREX0nySkiIvpO3nOKiIi+kzuniIjoSxkQERERfSXvOUVERN+ZXN16tVYll7SbpDskLZJ0VJftB0q6X9LN5eeQOuOJiGiu1ntOw/30JulsSfdJum2A7c+W9H1Jt0haIOmgkV/TwGq7c5I0FTgd2BVYDFwvaY7t2zuaftP24XXFERExMdR+53Qu8EXg/AG2HwbcbnsPSesBd0i6wPbSOoKps1tve2CR7TsBJF0E7AV0JqeIiBhUvcnJ9tWSZvRqAqwpScAawIPUOEKjzm69DYG725YXl+s6vVXSrZIulrRxjfFERDRY6z2n4X5G7IvAS4B7gPnAB22vGI0Dd1PnnZO6rOucZPj7wIW2l0g6FDgPeP0zDiTNAmaVi0s+pRO69olOAtOBB8Y7iHGSa5+cJvO1A7z4H1/vuQyOmT6CY60iaV7b8mzbs4ew/78CN1P8jX4RcLmkX9h+ZAQxDajO5LQYaL8T2ogi4z7F9l/bFs8ETup2oPIHOBtA0jzbM0c31GbItefaJ5vJfO1QXH/ru+3dxjMW4CDgRNsGFkm6C9gS+E0dJ6uzW+96YHNJm0p6FrAPMKe9gaT12xb3BBbWGE9ERAzfn4CdASQ9j+Ku7s66TlbbnZPtZZIOBy4DpgJn214g6Thgnu05wAck7UnxUO1B4MC64omIiIFJuhB4HTBd0mLgE8A0ANtnAMcD50qaT/HY5n/Zrq3LVcUdWnNImjXEftIJI9eea59sJvO1w+S+/sYlp4iImPhqrRARERExHH2bnCqUPlpZ0jfL7dcN8vJYo1S49iMl3V6+H/ZTSS8YjzjrMNi1t7XbW5IlTZiRXFWuXdLbyt/9AknfGOsY61Lh3/wmkq6UdFP57/6N4xFnHSqUDZKk/1P+bG6VtO1YxzgubPfdh2IAxf8ALwSeBdwCbNXR5v3AGeX3fSjKII177GN07TsBq5Xf3zeZrr1styZwNXAtMHO84x7D3/vmwE3AOuXyc8c77jG89tnA+8rvWwF/GO+4R/H6/wXYFrhtgO1vBH5EMQjh1cB14x3zWHz69c7pqdJHLuo2tUoftduL4qVdgIuBncuyGk036LXbvtJ2q3b+tRTvkE0EVX7vUIwaOhl4YiyDq1mVa38PcLrthwBs3zfGMdalyrUbWKv8/mw63plsMttXU4xWHshewPkuXAus3fEazoTUr8mpSumjp9rYXgY8DDxnTKKrV9WyTy0HU/xf1UQw6LVL2gbY2PYPxjKwMVDl974FsIWkX0q6VtJ4v5Q5Wqpc+yeB/cohznOB/xyb0PrCUP8mTAj9Op9TldJHVdo0UeXrkrQfMBPYsdaIxk7Pa5c0BTiVifk+XJXf+0oUXXuvo7hb/oWkl9r+W82x1a3Kte8LnGv785JeA3ytvPbaarv1kYn6t66nfr1zGrT0UXsbSStR3Or3ujVuiirXjqRdgI8De9peMkax1W2wa18TeCnwc0l/oOh/nzNBBkVU/Tf/PdtP2r4LuIMiWTVdlWs/GPgWgO1fA6tQ1N2bDCr9TZho+jU5DVr6qFw+oPy+N/Azl08PG65K2adtgK9QJKaJ8twBBrl22w/bnm57hu0ZFM/b9rQ9r/vhGqXKv/lLKQbDIGk6RTdfbeVjxlCVa28vnfMSiuR0/5hGOX7mAO8qR+29GnjY9r3jHVTd+rJbz9VKH32V4tZ+EcUd0z7jF/HoqXjtn6WYT+Xb5RiQP9nec9yCHiUVr31CqnjtlwFvkHQ7sBz4iJ9ePLmRKl77h4EzJX2IokvrwAnyP6NVygbNpRixtwh4jKIA64SXChEREdF3+rVbLyIiJrEkp4iI6DtJThER0XeSnCIiou8kOUVERN9Jcoq+ImljSXdJWrdcXqdcrqXyuqRDJb2r/H6gpA3atp0laatROs+bJB1bfj9X0t7DPM56kn48GjFF9LMkp+grtu8GvgycWK46EZht+481ne8M2+eXiwcCG7RtO8T27aN0qo8CXxrpQWzfD9wraYeRhxTRv5Kcoh+dCrxa0hHAa4HPdzaQNEPSbyWdV85xc7Gk1cptO5fz/swv58pZuVx/Yts8WJ8r131S0n+VdzIzgQsk3SxpVUk/b5VGkrRvebzbJJ3UFsffJX1a0i1lMdbndYl1C2CJ7Qe6bDu+vJOaIukPkj4j6deS5knaVtJlkv5H0qFtu10KvHP4P96I/pfkFH3H9pPARyiS1BHlNArdvJjirmpr4BHg/ZJWAc4F3m77ZRRVUN5XdhO+Gfinsv2nOs55MTAPeKftV9h+vLWt7Oo7CXg98ApgO0lvKjevDlxr++UUc0y9p0ucOwA3dq6UdDLwXOCgtgKmd9t+DfCL8jr2pqgheFzbrvOAfx7gZxIxISQ5Rb/aHbiXotDrQO62/cvy+9cp7rJeDNxl+3fl+vMoJnN7hGL+p7MkvYWiDExV2wE/t31/OT3LBeUxAZYCrek7bgBmdNl/fZ5ZB+5/A2vbfm9HGZ5Wiab5FJPKPVp25T0hae1y2320dT9GTERJTtF3JL0C2JXijuFDPSZW66y9ZbpPL9Ca82t74BLgTcBQBhX0msTyybbkspzu9SofpyhU2u564JWtgR9tWhXmV7R9by23jr1KecyICSvJKfqKikq2X6bozvsTRZHbzw3QfJNybh8o5vu5BvgtMEPSZuX6/YGrJK0BPNv2XOAIiu65To9STMvR6TpgR0nTJU0tz3XVEC5rIbBZx7ofUwz2+KGkbufsZQvgtiHuE9EoSU7Rb95DUWX98nL5S8CWkrpNqLgQOEDSrcC6wJdtP0FRtfnbkuZT3HGcQZF0flC2vQr4UJfjnQuc0RoQ0VpZTk/wMeBK4BbgRtvfG8I1XQ1sUybep9j+NnAmxZxUq3bds7udgB8OoX1E46QqeTSSpBnAD2z3eibVNySdBnzf9hWjcKyrgb1sPzTyyCL6U+6cIsbGZ4DVRnoQSesBpyQxxUSXO6eIiOg7uXOKiIi+k+QUERF9J8kpIiL6TpJTRET0nSSniIjoO0lOERHRd/4/t1AIO2UCqn4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the true and initial model and the perturbation between them\n",
    "plot_velocity(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo = np.load('/home1/wangrongqian/hw/project_data/common.npy', allow_pickle=True).item()\n",
    "for i in range(nshots):\n",
    "    names['data' + str(i)] = np.load('/home1/wangrongqian/hw/project_data/project2_data_shot_'+str(i)+'.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['time_values', 'source', 'src_coordinates', 'rec_coordinates'])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "geo.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2001, 204)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data1.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXSc9X3v8fdXM9JoteVFXuIFY2MWB4gBXQINbSBAr8lGm5I2tKW0h8anJyFrb3vISZvmpj09bW7TJNxyc+OQNKRJQ0i6QFI3kFCSm9sLBAPBgB2wbMA2GEt4l2akWfS9f8wz8lgejWZ/RtLndY6O5nnm0fP8HkvWR7/l+f3M3REREZlKS9gFEBGR5qagEBGRohQUIiJSlIJCRESKUlCIiEhR0bALUGuLFy/2NWvWhF0MEZEZ5fHHH3/N3fsKvTfrgmLNmjVs27Yt7GKIiMwoZvbSVO+p6UlERIpSUIiISFEKChERKUpBISIiRSkoRESkqFCDwsy+YmaDZvbMFO+bmd1uZgNmtt3MLm50GUVE5rqwaxRfBTYVef86YH3wsRn4QgPKJCIieUINCnf/P8DhIodcD3zNsx4Bes1seWNKF47joym+8n9fYGQsHXZRRESA8GsU01kB7Mvb3h/sO4WZbTazbWa2bWhoqGGFq4e/feB5PvW9Hdz+4K6wiyIiAjR/UFiBfaettOTuW9y93937+/oKPoE+Yzz484PB58GQSyIiktXsQbEfWJW3vRJ4JaSy1N2J0RT7DieIRVsYGBwmkcyEXSQRkaYPivuA3wlGP10GHHP3A2EXql4GBocB+K+vXwbAi4dGwiyOiAgQ/vDYbwIPA+eY2X4zu8XM/sDM/iA4ZCuwBxgAvgS8L6SiNsSeoWwwXLthKQAvvKagEJHwhTp7rLvfOM37Dry/QcUJ3avHRwG4bO0iAPYdjodZHBERoPmbnuaUg8dHmdcepa8nRmdbhMETY2EXSUREQdFMXj02yrL57QAs6YlxMKhhiIiESUHRRA6eGGPpvFxQtKtGISJNQUHRRA4eG50Iir55MYYUFCLSBBQUTcLdGRoeY0lPDMg2PQ2q6UlEmoCCokmcGEuTGXcWdLYB0NcTYySZIZ7UnE8iEi4FRZM4OpICoLezFWAiMI4lUqGVSUQEFBRN42giCUBvEBC9HdnAOBpXUIhIuBQUTeJIEAgLghrF/E4FhYg0BwVFkzgan1yjyDU9JUMrk4gIKCiaRq7mkOuj6FWNQkSahIKiSUwERdA3MT/XR6HObBEJmYKiSRyJJ+mJRYlGst+SzrYIrRHTqCcRCZ2CokkcS6To7Wqd2DYz5ne0qelJREKnoGgSR+PJiQ7snN7OVnVmi0joFBRN4vhomnkdpy4P0tvRqhqFiIROQdEkhkfTdMdODYr5Ha3qoxCR0CkomsTwWJruWOsp+7rbo4yMaa4nEQlX2GtmbzKz58xswMxuK/D+ajN7yMyeNLPtZvbWMMrZCMNjaXraT61RdMeiDCsoRCRkoQWFmUWAO4DrgA3AjWa2YdJhfwLc4+4XAe8B/ldjS9kY7s7wWJquWOSU/d2xKCdGFRQiEq4waxSXAgPuvsfdk8DdwPWTjnFgXvB6PvBKA8vXMKOpcTLjfnrTUyzKWHqcVGY8pJKJiIQbFCuAfXnb+4N9+T4J/LaZ7Qe2Ah9oTNEa68RYtsO6e3LTU7CtfgoRCVOYQWEF9vmk7RuBr7r7SuCtwD+Y2WllNrPNZrbNzLYNDQ3Voaj1NTKWAaBn0qinrmBbzU8iEqYwg2I/sCpveyWnNy3dAtwD4O4PA+3A4skncvct7t7v7v19fX11Km79DAdB0DUpKHLBoQ5tEQlTmEHxGLDezM40szayndX3TTpmL3A1gJmdRzYoZl6VYRoTTU+xwk1PCgoRCVNoQeHuaeBW4H5gJ9nRTc+a2afM7J3BYX8IvNfMngK+Cfyuu09unprxcjWKycNju1SjEJEmEJ3+kPpx961kO6nz930i7/UO4E2NLlej5YJgco1ioulJfRQiEiI9md0EcqOaJvdRqOlJRJqBgqIJnBgr3PTUrRqFiDQBBUUTGB5NE20xYtFTvx1dbapRiEj4FBRNYHgsTXd7FLNTHy1paTG62iIKChEJlYKiCQyPpidqD5N1t0fV9CQioVJQNIFCM8fmdGkGWREJmYKiCWTXoigcFD0KChEJmYKiCWSnGC8cFB1tERLJTINLJCJykoKiCcSTmdPWosjpbIsyklSNQkTCo6BoAolkho7WwjWKTtUoRCRkCoomEE+m6Wgr/K3obIsQV1CISIgUFE0gnszQOcXwWDU9iUjYFBQhy4w7Y+lxOlqn6qNQ05OIhEtBEbLRVDYEOtumDor0uJNMa91sEQmHgiJkuf6HqYMiGhyn5icRCYeCImS5ZqWOKfsosgGiDm0RCYuCImTxVLamMFUfRcdEUKhGISLhUFCEbLqmp66JpifVKEQkHAqKkI1OND1N3ZkNCgoRCU+oQWFmm8zsOTMbMLPbpjjm181sh5k9a2b/2Ogy1tu0ndkxdWaLSLgK96A2gJlFgDuAa4H9wGNmdp+778g7Zj3wMeBN7n7EzJaEU9r6iZcwPBZUoxCR8IRZo7gUGHD3Pe6eBO4Grp90zHuBO9z9CIC7Dza4jHWXCGoK7VN1ZrcqKEQkXGEGxQpgX972/mBfvrOBs83sP83sETPbVOhEZrbZzLaZ2bahoaE6Fbc+TjY9Tb1wEUBca1KISEjCDAorsM8nbUeB9cCVwI3AnWbWe9oXuW9x93537+/r66t5QespUWrTU0o1ChEJR5hBsR9Ylbe9EnilwDH3unvK3V8AniMbHLNGIpnBDGLRwt+KWLSFFkPzPYlIaMIMiseA9WZ2ppm1Ae8B7pt0zL8CVwGY2WKyTVF7GlrKOosnM3S2RjArVMECM8vOIDumoBCRcIQWFO6eBm4F7gd2Ave4+7Nm9ikze2dw2P3AITPbATwE/JG7HwqnxPURT2amfIYip6MtQiKlPgoRCUdow2MB3H0rsHXSvk/kvXbgo8HHrJRIpqcNii4tXiQiIdKT2SFLpDJ0TrEMak6Hmp5EJEQKipCV0vTUpaYnEQmRgiJkiWRmyqGxOR1tEdUoRCQ0CoqQxUsICi2HKiJhUlCELJHKTDl9R05XW3Ri3QoRkUZTUISs1KanuJqeRCQkCoqQxZPpKed5yunU8FgRCZGCImSJVCkP3EVJpDKMj0+eCktEpP4UFCFKZcZJZZzOafoock1TY+nxRhRLROQUCooQ5WaOna5GcXLxInVoi0jjKShClJhmvewcLV4kImFSUIRouvWyc3Kd3QmtSSEiIVBQhCjXlNQx7VxPLcHxCgoRaTwFRYgSJdYockGiPgoRCYOCIkTldmZrGg8RCYOCIkS5pqSOEofHqulJRMIwbVCYWaeZ/amZfSnYXm9mb69/0Wa/kpueVKMQkRCVUqP4e2AMuDzY3g/8Rd1KNIecHPU03RQeGvUkIuEpJSjWufungRSAuycAq8XFzWyTmT1nZgNmdluR424wMzez/lpct1lMjHoq+YE7BYWINF4pQZE0sw7AAcxsHdkaRlXMLALcAVwHbABuNLMNBY7rAT4IPFrtNZvNaKq0pqdYtAWz7PraIiKNVkpQ/BnwfWCVmX0DeBD44xpc+1JgwN33uHsSuBu4vsBxfw58GhitwTWbSjyZIdpitEaKfxvMjM5WzSArIuGYNijc/QfAu4DfBb4J9Lv7j2pw7RXAvrzt/cG+CWZ2EbDK3b9X7ERmttnMtpnZtqGhoRoUrTFKWS87p6MtQlx9FCISgil7Uc3s4km7DgSfV5vZand/osprF+rnmJhH28xagM+SDaii3H0LsAWgv79/xszFXcqiRTkdWg5VREJSbLjNZ4LP7UA/8BTZX+4Xku0vuKLKa+8HVuVtrwReydvuAc4HfmRmAMuA+8zsne6+rcprN4V4KjPtiKecztaonswWkVBM2fTk7le5+1XAS8DF7t7v7pcAFwEDNbj2Y8B6MzvTzNqA9wD35V3/mLsvdvc17r4GeASYNSEB2RrFdA/b5XRolTsRCUkpndnnuvvTuQ13fwbYWO2F3T0N3ArcD+wE7nH3Z83sU2b2zmrPPxMkUumS+yg61fQkIiEppd1jp5ndCXydbB/Cb5P9xV41d98KbJ207xNTHHtlLa7ZTOLJDN2xEpue2iIcjafqXCIRkdOVUqP4PeBZ4EPAh4EdwT6pUnlNT9GJ5y5ERBpp2j9n3X2U7Oijz9a/OHNLvIxRT3qOQkTCMm1QmNkL5A1bzXH3tXUp0RySSGXoKHHUU7YzW6OeRKTxSvktlT+/UjvwbmBhfYozt5T9HIWankQkBKU8mX0o7+Nld/8c8JYGlG1Wc3fiyXTJfRSdrRFSGSeVGa9zyURETlVK01P+E9otZGsYPXUr0Rwxlh5n3KefOTanI28G2fkdWm9KRBqnlKanz+S9TgMvAL9en+LMHaUuWpQzsSZFMsP8jta6lUtEZLJSguIWd9+Tv8PMzqxTeeaMRIlTjOecXJNCHdoi0liltGF8p8R9UoaJ9bLLGPWU/3UiIo1SbPbYc4HXA/PN7F15b80jO/pJqpBreiq5Mzu3brZGPolIgxX7c/Yc4O1AL/COvP0ngPfWs1BzQa4JqdymJ833JCKNNmVQuPu9wL1mdrm7P9zAMs0JuUWISh711Jr9VqnpSUQarVjT0x+7+6eB3zSzGye/7+4frGvJZrnRskc95Zqe1JktIo1VrOkpN0PsrFn/oZnkagadrerMFpHmVqzp6bvB57saV5y5I9f01N5W2sNzHeqjEJGQFGt6+i4FJgPMcfc5sbhQvSQmOrNLXQpVNQoRCUex31J/07BSzEHxMofHRiMttEVaFBQi0nDFmp5+nHsdrGl9LtkaxnPunmxA2Wa1RDJDLNpCpMVK/pqOtshETUREpFGmbSA3s7cBu4Hbgb8DBszsulpc3Mw2mdlzZjZgZrcVeP+jZrbDzLab2YNmdkYtrtsMylm0KKezTYsXiUjjldKT+hngKne/0t3fDFxFDVa7M7MIcAdwHbABuNHMNkw67Emg390vJDttyKervW6zyAZFaf0TOVqTQkTCUEpQDLr7QN72HmCwBte+FBhw9z1BU9bdwPX5B7j7Q+4eDzYfAVbW4LpNIZFKl/ywXU5nW0SjnkSk4Ur5k/ZZM9sK3EO2j+LdwGO5+Z/c/Z8rvPYKYF/e9n7gjUWOvwX49wqv1XQqanpqjarpSUQarpSgaAcOAm8OtofILoX6DrLBUWlQFOrFLTgc18x+m+yCSW+e4v3NwGaA1atXV1icxoonMyWPeMppb4twLJGqU4lERAqbNijc/ffqdO39wKq87ZXAK5MPMrNrgI8Db3b3sUIncvctwBaA/v7+KZ/9aCaJZIbF3W1lfU1na4RXjyXqVCIRkcJKWQr1TOADwJr842vwwN1jwPrg/C8D7wF+c9K1LwK+CGxy91r0izSNeDJNZ1tnWV+jUU8iEoZSmp7+Ffgy8F1gvFYXdve0md0K3A9EgK+4+7Nm9ilgm7vfB/wPoBv4tpkB7J0tT4QnkpmyO7M71JktIiEoJShG3f32elzc3bcCWyft+0Te62vqcd1mkEiV30ehGoWIhKGUoPi8mf0Z8AAw0Ufg7k/UrVRzQCWjnjraoiRSGcbHnZYynugWEalGKUFxAXAT8BZONj15sC0VyIw7Y+nxip6jABhNl/+wnohIpUr5bfOrwFrN71Q7uaerK5nCA7L9GwoKEWmUUp7MforsutlSI7n1sjvKncJDU42LSAhK+U21FPi5mT3GyT4Kd/fri3yNFJGYWN2u3BpF9tul+Z5EpJFKCYo/y3ttwBXAaWtoS+niZa6XndMRrIanGoWINNK0TU/BuhTHgLcBXwWuBv53fYs1u00sWlRuUATra8e1JoWINFCxpVDPJvu09I3AIeBbgLn7VQ0q26w10fRUZh9Fp9bNFpEQFPtN9XPgJ8A7ctOMm9lHGlKqWS4+sV52ZaOe1PQkIo1UrOnp14BXgYfM7EtmdjWFZ3yVMuU6oyuZwgNUoxCRxpoyKNz9X9z9N8iulf0j4CPAUjP7gpn9coPKNytV2pmda6pSH4WINFIpndkj7v4Nd3872anAfwactr61lG6iM7uCuZ4AEqmazc0oIjKtUh64m+Duh939i+6u6TuqkJh44K68oIhFWzA7+fUiIo1QVlBIbcSTGSItRlukvH9+M6OzVTPIikhjKShCEE9m6GyNEKyxUZaOtghxPZktIg2koAhBJYsW5WjxIhFpNAVFCOKp8teiyOlsjWrUk4g0lIIiBIlkuuyZY3M6tMqdiDSYgiIEiWpqFGp6EpEGCzUozGyTmT1nZgNmdtqzGWYWM7NvBe8/amZrGl/K2qtkGdQcrZstIo0WWlCYWQS4A7gO2ADcaGYbJh12C3DE3c8CPgv8dWNLWR+JZKbsh+1ycutmi4g0SpjraV4KDLj7HgAzuxu4HtiRd8z1wCeD198B/s7MzN291oUZTWV479e2saK3gxW9HZy/Yj6XnrmQrljt/4mqqlG01q/paXgszcO7D7HzwHH2Ho5zNJ4inkwzksyQGR/HHdxhPPjndwen5t8KEanQucvmcfuNF9X8vGEGxQpgX972fuCNUx3j7mkzOwYsAl7LP8jMNgObAVavXl1RYU6MpjmeSLHzwAleG84u5BeLtnDDJSv5yLVns7g7VtF5C4knM1V2Ztd21NPg8VE++8Pn+acnXiaZzk4PsnRejIVdMbpjEeZ3tBJtMYzsQ39mBK/ByG6LSPhWLOioy3nDDIpCv14m/3layjG4+xZgC0B/f39Ff+L29cS499YrgOyke0+8dJTvPvUK3962n/uffZUv3nQJl5yxsJJTn2ZkLE1XNZ3ZNWx6enj3Id73jccZSWa44ZKVvOPC17FxVW/Fz3mIyOwTZmf2fmBV3vZK4JWpjjGzKDAfOFzvgnW2Rbli/WL++oYL+e4HrmBeeys3ffmnPPPysarPnRl3EqlMxU1aHa0RUhknlal+YsAn9x7h9776UxZ1x9j6wV/kL3/1Ai5ft0ghISKnCDMoHgPWm9mZZtZGdjW9+yYdcx9wc/D6BuA/6tE/Ucw5y3q4e/NlzGtv5f3/+ATDY9U1++SajborDYoaLV50LJHi1n98ksXdMe7efBlnLemu6nwiMnuFFhTungZuBe4HdgL3uPuzZvYpM3tncNiXgUVmNgB8lJCmN18yr53Pv2cjLx2K84UfDVR1rpGx7C/4SmsUuTUpqu3Q/p8P7uKVYwluv/Gimva/iMjsE2YfBe6+Fdg6ad8n8l6PAu9udLkKeePaRfzKxtfxpZ+8wM2Xr2HJvPaKzpOrkXTFKu+jgOoWL3r5aIK7Hn6Rd1+ykotXL6j4PCIyN+jJ7DJ8+JqzSWXG+fojL1V8jpFcUFQx6gmqa3r6h4dfIjPufOiasys+h4jMHQqKMqxZ3MXV5y7l64/unRhGWq6JoKi46Sm3yl1lQTGayvCtx/byyxuWsaK3PkPpRGR2UVCU6TffuIrDI0l+smuooq8fCWoClXZmTwRFhTWKH+48yJF4ipsuP6OirxeRuUdBUaYrzuqjt7OV7z41eSRvaUaq7KPoaM0GTKVNT99/5lUWd7dx2dpFFX29iMw9CooytUVb2PT6Zfxgx8GKnmXIdWZXXaNIld+ZPZrK8NDPB7l2wzIiLXqcWkRKo6CowJXnLGEkmeGJl46U/bXV9lFU05n9/3a/xkgyw6bzl1V0bRGZmxQUFfiFsxYRaTF+suu16Q+eJBcUlc4em6uJjFTw4N9/DhyiLdrCG8+szVQkIjI3KCgqMK+9lYtW9VbUoT2SzNDVFqGlwqafzrYIZjA8Wn5QPLz7EJesXkB7hSElInOTgqJCV6xfzPaXj3F8NFXW142MpauautzM6I5FOVFmjeJoPMnOV49z+Tp1YotIeRQUFbrkjAW4w1P7jpb1dcNj6Yo7snN6YlFOlFmjePSFw7ijoBCRsikoKrRxVS9m8MRL5QVFtTUKgJ721rKbnp7ce5TWiHHhyvlVXVtE5h4FRYV62ls5e0kPT+wtb+TTyFim4mcocrrbo2XPYrt9/1HOWz6PWFT9EyJSHgVFFS5a3cuTe48wPl76zOcjyXTF8zzllNtHMT7uPL3/GBesUG1CRMqnoKjCxlW9HB9Ns/dwvOSvqUXTU3d7lBNldKK/cGiEE2Np3rCyt6rrisjcpKCownnL5wGw88Dxkr9meKzy1e1y5rVHy+qjeHp/dmW+C1epRiEi5VNQVOGcZT20WHlBMTKWprvaPopYeX0U2/cfo721hbP6tIqdiJRPQVGF9tYIa/u62VFiUKQz4yRSGXraW6u6bneslXgyQ6bEvpFdgyc4e2kP0Yi+3SJSPv3mqNJ5y+ex88CJko7NPfswr736Pgoo/ensgcFh1qk2ISIVCiUozGyhmf3AzHYFn09bj9PMNprZw2b2rJltN7PfCKOs09mwfB4vH01wNJ6c9tjcU9zV1ih6gj6OE2PTd2gPj6U5cGyUs5YoKESkMmHVKG4DHnT39cCDwfZkceB33P31wCbgc2bWdMN2zlveA1BSreJ4IqhRdFQZFLkaRQn9FC8MjQCwrq+rqmuKyNwVVlBcD9wVvL4L+JXJB7j78+6+K3j9CjAI9DWshCU6e2k2KAaGhqc9NlejaGTT08BQNsBUoxCRSoUVFEvd/QBA8HlJsYPN7FKgDdg9xfubzWybmW0bGqpsidJKLZ/fTmdbhN2D0wdF7tmHamsUubmiSpnvaffgCJEWY/VC1ShEpDLV/WlbhJn9ECi0Qs7HyzzPcuAfgJvdveCScu6+BdgC0N/fX/pj0jVgZqzr62Z3KTWKGjc9lfJ09sDgMGcs7KQtqnELIlKZugWFu18z1XtmdtDMlrv7gSAIBqc4bh7wb8CfuPsjdSpq1c5a0s0jew5Ne9zJzuzqJwWE0pqedg8Ns07NTiJShbD+zLwPuDl4fTNw7+QDzKwN+Bfga+7+7QaWrWxnLenmwLHRaTuXjydSmEF3DeZ6AhieZtRTOjPOi4dGNDRWRKoSVlD8FXCtme0Crg22MbN+M7szOObXgV8CftfMfhZ8bAynuMXlRhTtmab56fhomp5YtOLV7XI62yJEWoxjieJBsfdwnFTGNeJJRKpSt6anYtz9EHB1gf3bgN8PXn8d+HqDi1aR3IiigcFhLiwy8d7x0VTV/ROQ7Rfp7WjlaLx4UOwOhsZqxJOIVEM9nDWwemEXkRabtkP7eCLNvCoftsuZ39nK0WlqFLnyqI9CRKqhoKiBtmgLZyzqZGCaIbLHR1NVd2TnLOhsm/Zp8IHBYZb0xGoWTiIyNykoauSsvu7pgyJRm6YnoMSmJ83xJCLVU1DUyLol3UHnccFHPYDsA3K1qlHM7yweFO6enQxwiTqyRaQ6CooaWbu4i1TG2VdktbvDI0kWdrbV5HrTNT0NDY9xYjStNShEpGoKihrJdRjvCUYaTZZIZkikMizoqk1Q9Ha0MpLMkEwXrsHsHhw5pVwiIpVSUNTIusXZX8hTjXw6Evz1v7BWQRGcZ6pnKXKTFGporIhUS0FRI/M7W1ncHZsyKA6P1Dgogk7xqZqfdg8O09kWYdm89ppcT0TmLgVFDa3t65p4yG2ymgdFZxAUU9QociOezKp7ClxEREFRQ+v6uqecxiPX9LSghp3ZwJQjn3YPDqvZSURqQkFRQ+v6ujgST03UHvLVukYxP2h6OlKg6WlkLM0rx0Y1x5OI1ISCooZyI4wK9VMcGUlidvIXfLUWd8cAeG147LT39miOJxGpIQVFDeWeWSi02t2hkSQLOtuIVDlzbE5HW4SeWJTB46cHxcQcT3qGQkRqQEFRQ6/r7aAt2sKe107v0D40nGRBZ23nXOqbF2PoxOlBMTA4TKTFOGORmp5EpHoKihqKtBhrF3cVrFEcPDHKsvm1Haq6pCfG4InR0/bvGjzBGYu0/KmI1IZ+k9TYVOtnDx4fY2mNn2no62lnsECNYtfgMOvVPyEiNaKgqLF1fV3sO5JgLJ2Z2Dc+7hw8PlrzoFjSE2Pw+BjuPrFvLJ3hpUNxzl7aU9NricjcpaCosbV93WTGnb2HTk4OeDieJD3uNX9KeklPjEQqc8pa3S+8NkJm3FmvoBCRGgklKMxsoZn9wMx2BZ8XFDl2npm9bGZ/18gyVio30ii/+eng8Ww/wtJ5sZpea0lwvoN5I5+eP5i9rpqeRKRWwqpR3AY86O7rgQeD7an8OfDjhpSqBtYGD7nlL2L0ytFsUCyf31HTa61c0AnAviMnay8DB09kO9X1sJ2I1EhYQXE9cFfw+i7gVwodZGaXAEuBBxpUrqp1xaKsXtjJjgPHJ/a9dCg7XHZNjYernrEwGxT5zVw7Xz3BmkWdxKKRml5LROausIJiqbsfAAg+L5l8gJm1AJ8B/qjBZavaBSvns33/sYntFw+N0NvZyvxaP0fRE6O9tYW9eYslbd9/lAtX9tb0OiIyt9VmXc4CzOyHwLICb328xFO8D9jq7vummwHVzDYDmwFWr15dTjHr4sIV8/m37Qc4NDzGou4YL74Wr8vDb2bG6oWdvBTUKA4eH+Xg8TEuWDG/5tcSkbmrbkHh7tdM9Z6ZHTSz5e5+wMyWA4MFDrsc+EUzex/QDbSZ2bC7n9af4e5bgC0A/f39Pvn9Rsv9Rf/0y8e48pwl7Bka5tIzF9blWmct6eaZl7PNXLlazIUrFRQiUjthNT3dB9wcvL4ZuHfyAe7+W+6+2t3XAP8N+FqhkGhG56+YR6TF2PbiEV4bHuOVY6O8/nX1+eV9wYpe9h6Oc2QkySN7DtEWbeF81ShEpIbCCoq/Aq41s13AtcE2ZtZvZneGVKaa6Wlv5ZLVC/jR84P8bO9RINtvUQ+52sP2l4/xk11DvPHMhbS3qiNbRGonlKBw90PufrW7rw8+Hw72b3P33y9w/Ffd/dbGl7RyV527hGdePs6Wn+yhOxZl46r6dDC/YVUvsWgLdzw0wPMHh7nynNPGBYiIVEVPZtfJu/tXEou28NMXDrPp/GV1+yu/OxZl0/nL+OkLh+lsi/Cui1bU5ToiMnfVrd5ZmpEAAAcaSURBVDN7rlvcHeOLN13Cj58f4tarzqrrtf707RtY0NnG1ectYUGNVtATEcmx/AnlZoP+/n7ftm1b2MUQEZlRzOxxd+8v9J6ankREpCgFhYiIFKWgEBGRohQUIiJSlIJCRESKUlCIiEhRCgoRESlKQSEiIkXNugfuzGwIeKmKUywGXqtRccI0W+4DdC/Narbcy2y5D6juXs5w975Cb8y6oKiWmW2b6unEmWS23AfoXprVbLmX2XIfUL97UdOTiIgUpaAQEZGiFBSn2xJ2AWpkttwH6F6a1Wy5l9lyH1Cne1EfhYiIFKUahYiIFKWgEBGRohQUATPbZGbPmdmAmd0WdnmmY2ZfMbNBM3smb99CM/uBme0KPi8I9puZ3R7c23Yzuzi8kp/KzFaZ2UNmttPMnjWzDwX7Z+K9tJvZT83sqeBe/nuw/0wzezS4l2+ZWVuwPxZsDwTvrwmz/IWYWcTMnjSz7wXbM/JezOxFM3vazH5mZtuCfTPxZ6zXzL5jZj8P/s9c3oj7UFCQ/c8A3AFcB2wAbjSzDeGWalpfBTZN2ncb8KC7rwceDLYhe1/rg4/NwBcaVMZSpIE/dPfzgMuA9wf/9jPxXsaAt7j7G4CNwCYzuwz4a+Czwb0cAW4Jjr8FOOLuZwGfDY5rNh8CduZtz+R7ucrdN+Y9ZzATf8Y+D3zf3c8F3kD2e1P/+3D3Of8BXA7cn7f9MeBjYZerhHKvAZ7J234OWB68Xg48F7z+InBjoeOa7QO4F7h2pt8L0Ak8AbyR7JOy0ck/a8D9wOXB62hwnIVd9rx7WBn84nkL8D3AZvC9vAgsnrRvRv2MAfOAFyb/uzbiPlSjyFoB7Mvb3h/sm2mWuvsBgODzkmD/jLi/oLniIuBRZui9BE01PwMGgR8Au4Gj7p4ODskv78S9BO8fAxY1tsRFfQ74Y2A82F7EzL0XBx4ws8fNbHOwb6b9jK0FhoC/D5oD7zSzLhpwHwqKLCuwbzaNG276+zOzbuCfgA+7+/FihxbY1zT34u4Zd99I9q/xS4HzCh0WfG7aezGztwOD7v54/u4Chzb9vQTe5O4Xk22Oeb+Z/VKRY5v1XqLAxcAX3P0iYISTzUyF1Ow+FBRZ+4FVedsrgVdCKks1DprZcoDg82Cwv6nvz8xayYbEN9z9n4PdM/Jectz9KPAjsv0uvWYWDd7KL+/EvQTvzwcON7akU3oT8E4zexG4m2zz0+eYmfeCu78SfB4E/oVsiM+0n7H9wH53fzTY/g7Z4Kj7fSgosh4D1gcjOtqA9wD3hVymStwH3By8vplse39u/+8EoyAuA47lqqphMzMDvgzsdPe/zXtrJt5Ln5n1Bq87gGvIdjY+BNwQHDb5XnL3eAPwHx40JofN3T/m7ivdfQ3Z/w//4e6/xQy8FzPrMrOe3Gvgl4FnmGE/Y+7+KrDPzM4Jdl0N7KAR9xF2B02zfABvBZ4n26b88bDLU0J5vwkcAFJk/3K4hWyb8IPAruDzwuBYIzuqazfwNNAfdvnz7uMKstXh7cDPgo+3ztB7uRB4MriXZ4BPBPvXAj8FBoBvA7Fgf3uwPRC8vzbse5jivq4EvjdT7yUo81PBx7O5/98z9GdsI7At+Bn7V2BBI+5DU3iIiEhRanoSEZGiFBQiIlKUgkJERIpSUIiISFEKChERKSo6/SEic4OZ5YYZAiwDMmSnTACIu/sv1OGaFwHvd/ffr/I8twIj7v73tSmZyEkaHitSgJl9Ehh297+p83W+DfyFuz9V5Xk6gf/07NQOIjWlpieREpjZcPD5SjP7sZndY2bPm9lfmdlvWXYdiqfNbF1wXJ+Z/ZOZPRZ8vKnAOXuAC3MhYWafNLO7zOyBYP2Ed5nZp4Pzfj+Y6oTgmjuCNQb+BsDd48CLZnZpo/5NZO5QUIiU7w1k12m4ALgJONvdLwXuBD4QHPN5sus2/Bfg14L3Jusn+wR3vnXA24Drga8DD7n7BUACeJuZLQR+FXi9u18I/EXe124DfrH62xM5lfooRMr3mAdz5pjZbuCBYP/TwFXB62uADdmprACYZ2Y97n4i7zzLOdkHkvPv7p4ys6eBCPD9vHOvIbsuxChwp5n9W7CdMwicW+W9iZxGQSFSvrG81+N52+Oc/D/VQnYhn0SR8yTIzpF02rndfdzMUn6yE3Gc7IJB6aB56Wqyk/XdSnZmV4JzFbueSEXU9CRSHw+Q/SUOgJltLHDMTuCsck4arNsx3923Ah8mO0lcztmc3pQlUjUFhUh9fBDoDzqcdwB/MPkAd/85MD83BXaJeoDvmdl24MfAR/LeexPwwyrKLFKQhseKhMjMPgKccPdCnd3lnOci4KPuflNtSiZykmoUIuH6Aqf2eVRqMfCnNTiPyGlUoxARkaJUoxARkaIUFCIiUpSCQkREilJQiIhIUQoKEREp6v8DzvUsslFVq0cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from examples.seismic import AcquisitionGeometry\n",
    "\n",
    "# First, position source centrally in all dimensions, then set depth\n",
    "src_coordinates = geo['src_coordinates']\n",
    "\n",
    "# Define acquisition geometry: receivers\n",
    "\n",
    "# Initialize receivers for synthetic and imaging data\n",
    "rec_coordinates = geo['rec_coordinates']\n",
    "\n",
    "# Geometry\n",
    "\n",
    "geometry = AcquisitionGeometry(model, rec_coordinates, src_coordinates, t0, tn, f0=f0, src_type='Ricker')\n",
    "# We can plot the time signature to see the wavelet\n",
    "geometry.src.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "src_coordinates[:, 1] = 500 - src_coordinates[:, 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "rec_coordinates[:, 1] = 500 - rec_coordinates[:, 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define rtm operator for imaging\n",
    "from devito import TimeFunction, Operator, Eq, solve\n",
    "from examples.seismic import PointSource\n",
    "\n",
    "def ImagingOperator(model, image):\n",
    "    objective=0\n",
    "    # Define the wavefield with the size of the model and the time dimension\n",
    "    v = TimeFunction(name='v', grid=model.grid, time_order=2, space_order=4)\n",
    "\n",
    "    u = TimeFunction(name='u', grid=model.grid, time_order=2, space_order=4,\n",
    "                     save=geometry.nt)\n",
    "    \n",
    "    # Define the wave equation, but with a negated damping term\n",
    "    eqn = model.m * v.dt2 - v.laplace + model.damp * v.dt.T\n",
    "\n",
    "    # Use `solve` to rearrange the equation into a stencil expression\n",
    "    stencil = Eq(v.backward, solve(eqn, v.backward))\n",
    "    \n",
    "    # Define residual injection at the location of the forward receivers\n",
    "    dt = model.critical_dt\n",
    "    residual = PointSource(name='residual', grid=model.grid,\n",
    "                           time_range=geometry.time_axis,\n",
    "                           coordinates=geometry.rec_positions)    \n",
    "    res_term = residual.inject(field=v.backward, expr=residual * dt**2 / model.m)\n",
    "    \n",
    "    # Correlate u and v for the current time step and add it to the image\n",
    "    image_update = Eq(image, image - u * v)\n",
    "    return Operator([stencil] + res_term + [image_update],\n",
    "                    subs=model.spacing_map)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "from examples.seismic.acoustic import AcousticWaveSolver\n",
    "solver = AcousticWaveSolver(model, geometry, space_order=4)\n",
    "from examples.seismic import plot_shotrecord"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAADPCAYAAABRPx4bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAfyklEQVR4nO3de7gdZX328e+dEAE5CIgHzkEIKlUUDagXtoiIL9gCHmgLKgpF4wFqEWuriBjQCljF8iqKW+SkCCpYjRqlWhE8gQkQCBApKahE8AWEAhZyvt8/ZlZYrKy99uzD7L1m7/tzXetizcwzz/wmybV+zMwzv0e2iYiI6CfTJjqAiIiITklOERHRd5KcIiKi7yQ5RURE30lyioiIvpPkFBERfafW5CTpQEm3SVoq6QNdtm8o6Wvl9mslzawznoiI6E7SDpKulLRE0i2S/qFLmzdJuqn8/ELSC+qKp7bkJGk6cDZwELA7cISk3TuaHQM8aHtX4NPAGXXFExERPa0G3mf7ucBLgWO7/GbfCexrew/go8BAXcHUeeW0N7DU9h22VwKXAod2tDkUuLD8fhmwvyTVGFNERHRh+x7b15ffHwGWANt1tPmF7QfLxWuA7euKp87ktB1wV9vyMjpOtL2N7dXAQ8BTa4wpIiKGUD5i2RO4tkezY4Dv1xXDBnV1DHS7AuqslVSlDZLmAHMAZmwy48VbPyf5KyImv3uu+8P9tp8GsKvkR0fTF9wCLG9bNWB7vdtykjYFLgeOt/1wt74k7UeRnF4+ipB6qjM5LQN2aFveHrh7kDbLJG0APAV4oLOj8g9wAGDb2dv4mIVH1xJwREQ/+ZhO+23r+2PAsaPo6yRYbnt2rzaSZlAkpottf3OQNnsA5wIH2f7jKELqqc7ktACYJWln4PfA4cAbO9rMA94K/BI4DPixU4k2ImI904CNa+y/fN7/JWCJ7TMHabMj8E3gSNv/VWM49SUn26slHQdcAUwHzrN9i6RTgYW251H8QXxZ0lKKK6bD64onIqLJBMyo9xD7AEcCiyUtKtedCOwIYPsc4GSKcQGfK8eurR7qamyk6rxywvZ8YH7HupPbvi8H/rrOGCIiJoO6r5xs/4zu4wDa27wNeFuNYaxTa3KKiIixIabWD/ZUOteIiMaq+8qp3yQ5RUQ0wDTgyRMdxDhKcoqIaIip9IM9lc41IqKxclsvIiL6TpJTRET0nXF4z6mvJDlFRDRArpwiIqLv5D2niIjoOxlKHhERfUfAxqP5xV49VpGMjySniIgGkGCDJKeIiOgn06bBxhuOooPlQzfpJ0lOERENIMGM0SSnhklyiohogik2XG8KnWpERIMJ2Giigxg/SU4REU0gijnFp4gkp4iIJphGrpwiIqLPCMiAiIiI6CsZEBEREX0nV04REdF38swpIiL60hQarTetzs4lHSjpNklLJX2gy/a/kHS9pNWSDqszloiIRmtdOY300zC1JSdJ04GzgYOA3YEjJO3e0ex3wFHAV+uKIyJiUmi95zTST8PUeVtvb2Cp7TsAJF0KHArc2mpg+zfltrU1xhER0XypEDFmtgPualteBrykxuNFRExe08hovTGiLus8oo6kOcAcgM133Hw0MUVENFPecxozy4Ad2pa3B+4eSUe2B4ABgG1nbzOiBBcR0Wh5z2nMLABmSdoZ+D1wOPDGGo8XETF5TbH3nGobrWd7NXAccAWwBPi67VsknSrpEABJe0laBvw18AVJt9QVT0RE42W03tiwPR+Y37Hu5LbvCyhu90VERC81XzlJ2gG4CHgmsBYYsH1WR5vnAOcDLwI+ZPuTdcUzhR6vRUQ0WP0DIlYD77N9vaTNgOsk/dD2rW1tHgDeA7y21khIcoqIaIaaB0TYvge4p/z+iKQlFK8Etb+bei9wr6S/rC+SQpJTREQTjONoPUkzgT2Ba8fniOtLcoqIaILR39bbWtLCtuWB8jWdJx5G2hS4HDje9sOjOuIoJDlFRDTB6CtE3G97dq8GkmZQJKaLbX9zVEcbpSSniIgmqLm2niQBXwKW2D5zDPvdBFhue81w9ktyiohogvpH6+0DHAkslrSoXHcisCOA7XMkPRNYCGwOrJV0PLB7++0/SdMoii68CdgLWAFsKOk+ileLBmzfPlQwSU4REU1Qc+FX2z+je03U9jZ/YOh3U68EfgR8ELjZ9loASVsB+wGnS/p321/p1UmSU0REUzTjF/tVtld1rrT9AMXzrMvLZ1s9NeNUIyKmuoZMmdFKTJJ2AZbZXiHpFcAewEW2/6db8upU6zTtERExRlrvOY30M/4uB9ZI2pVioMXODGPW81w5RUQ0QfPmc1pre7Wk1wH/Zvszkm6ounOzTjUiYqpq3pQZqyQdAbwVOLhcN+SzppYkp4iIJmjeZINHA+8E/sX2neXcfj1H6LVLcoqIaIKG3NaTNAB8H/iR7fe01tu+Ezi9aj89T1XS9hQvU/05sC3wGHAz8D3g+63x6xERUbPmXDmdBxwInCBpJfAfwA9s3zicTgZNTpLOpyiX/l3gDOBeijueu5UH/pCkD9i+emTxR0REZQ25crJ9DXANMFfSU4FXA++TtAdwPUWi+vpQ/fQ61U/ZvrnL+puBb0p6EmVZi4iIqJnAzbhyWsf2H4FLyg+SXkxxcTOkQZPTIImpfftKYGn1MCMiYqQ8DVY2aLSepC2AtwAzacs17c+hehnyIlHSXwEfBXYq26vo35uPIN6IiBgBC1ZPH03dhHEfIjCf4vbe4pEcvModzH8DXg8stu3hHiAiIkbPEis3HM19vcfGLJaKNrJ9wkh3rpKc7qKoLJvEFBExQdYyjRU8aRQ9jHty+rKkt1MMqlvRWlkWgB1SleT0T8B8SVd1HGDIyagkHQicBUwHzrV9esf2E4C3AauB+4C/s/3bKoFHREw1a5owXO9xK4F/BT4EtC5uDDyrys5VzvRfgD9RDCOvnLYlTQfOBg4AlgELJM2zfWtbsxuA2bYflfQu4BPA31Y9RkTEVLGWaawc1ZXTuDsB2NX2/SPZuUpy2sr2q0fQ997AUtt3AEi6FDgUWJecbF/Z1v4a4M0jOE5ExJSwhukTHcJw3AI8OtKdqySnH0l6te3/GGbf21E8r2pZBrykR/tjKEpeREREh+KZU6NedFoDLJJ0JU98JDQ2Q8mBY4F/krQCWEX1oeTdpvvtOqhC0puB2cC+g2yfA8wB2HzHjGCPiKnHqGm39b5VfkZkyORke7POdZJ6zjNfWgbs0La8PXB3l75eRfHAbF/bKzq3lzEMAAMA287eJqMGI2LKMbC6Wbf1brZ9XfsKSQcP1rjTkG90STq1Y3ka1cqeLwBmSdq5LHV0ODCvo689gS8Ah9i+t2rQERFTjZnGSjYc8WcCfFHS81sL5dxOJ1XducrrxjtK+mDZ+YYUl2m3D7WT7dXAccAVwBLg67ZvkXSqpEPKZv8KbAp8Q9IiSfMG6S4iYspbw/QRfybAYcCFkp5bvu/0booisJVUeeZ0NHBxmaD2o5gq49NVOrc9n6KERfu6k9u+v6pqoBERU9noX8IdX7bvkHQ4xQXNXcCrbVd+E7jXlBkvals8i+L228+BqyS9yPb1I4w5IiKGqRgQ0f+j9SQt5omD37aiKMRwrSRs71Gln55TZnQsPwjsXq438Mrq4UZExGgYNeU9p78ai056TZmx31gcICIiRs+oKe85/dH2n3o1kLTpUG0GHRAh6c29hoxL2kXSy4eOMyIiRqv1ntNIP+Po25I+JekvJG3SWinpWZKOkXQFFSYc7HVb76kUb/deB1xHUZh1I2BXipdl7wc+MJoziIiIaowa8Z6T7f0lvQZ4B7CPpC0pinvfBnwPeKvtPwzVT6/bemdJ+izFs6V9gD0oaq4vAY60/bvRn0ZERFTRlAER0H2k9nD1HEpuew3ww/ITERETpEEDIsZEoyYHiYiYqooBEfU9O5K0A3AR8EyKadUHbJ/V0UYUrxa9hqLi+FF1vVaU5BQR0QBry/JFNVoNvM/29ZI2A66T9MOOOfgOAmaVn5cAn6f3bBMjVqV8UURE9IE6yxfZvqd1FWT7EYrxBdt1NDsUuMiFa4AtJG3TrT9Jn5T0ZyM91yGvnMp6em8AZra3t33qYPtERMTYGs/5nCTNBPYEru3Y1G2evu2Ae7p082tgQNIGwPnAJbYfqhpDldt63wYeohhO3nVKi4iIqNcYzOe0taSFbcsD5XRETyBpU+By4HjbD3du7hpat5X2ucC5kp5NUaP1Jkk/B77YMQt6V1WS0/a2h3xhKiIi6jXK0Xr3257dq4GkGRSJ6WLb3+zSpNI8fW39TQeeU37uB24ETpD0DtuH94qlSnL6haTn215coW1ERNSg7tt65Ui8LwFLbJ85SLN5wHGSLqUYCPGQ7W639JB0JnAw8GPg47Z/VW46Q9JtQ8XTqyp5q7LsBsDRku6guK3Xmqa9UmXZiIgYGzW/57QPcCSwWNKict2JwI4Ats+heLH2NcBSiqHkR/fo72bgJNuPdtm291DB9LpyGpPKshERMXp1z+dk+2d0f6bU3sbAsRW7fJPt89pXSPpP2/tXGRjRq3zRb8vOvmz7yI4DfJkiw0ZExDhoSvkiSRsBT6YYgLEljye8zYFtq/ZT5ZnTE8aplw+4Xlz1ABERMXoNKl/0DuB4ikTUXj3iYeDsqp30eub0QYr7jRtLepjHs99KYL3hhxERUZ+6yxeNlbLk0VmS/t72Z0baT6/beqcBp0k6zfYHR3qAiIgYvQbd1nul7R8Dv5f0+s7tgwxRX0+V23onlgd4OcXovZ/a/tawoo2IiFFp0G29fSmGjx/cZZuBMUtOZ1NMMHhJufxOSQfYrjpiIyIiRqkp07Tb/kj5317DzIdUJTntCzyvHEKIpAuBvJAbETGOGnTlBICkjwOfsP0/5fKWFFXPT6qyf5Wq5LdRvoRV2gG4qWJwB0q6TdJSSetN6S7pnZIWS1ok6WeSdq/Sb0TEVNMaEDHSzwQ4qJWYAGw/SPECbyVVrpyeCiyR1Co9sRfwS0nzygMe0m2ncsj52cABFPWYFkia1zE3yFfLt46RdAhwJpA6fhERHcZhPqexNl3ShrZXAEjaGKqfQJXkdPIIA9sbWGr7jjKwSynmAlmXnDoq3m7CINVtIyKi9vJFY+0rwH9KOp/it/3vgAur7jxkcrJ9laSdgFm2f1Rmvw3Kyah66Tbvx3ozJko6FjgBeBLwyqqBR0RMJWsb8p5Ti+1PSLoJeFW56qO2r6i6f5XJBt8OzAG2AnahKJF+DrD/ULt2i3e9FfbZwNmS3gicBLy1SwxzyhjYfMfNhwo5ImLScfNu6wHcAMyg+O2/YTg7VhkQcSxFtdqHAWzfDjy9wn7DmvcDuBR4bbcNtgdsz7Y9e5OnPbnCoSMiJp86p2kfa5L+BvgVcBjwN8C1kg6run+VZ04rbK8spvqAcsrdKs+GFgCzJO0M/B44HHhjR/CzymQH8JfA7URExHrGc5r2MfIhYC/b9wJIehrwI+CyKjtXSU5XSWrV2DsAeDfwnaF2sr1a0nHAFcB04Dzbt0g6FVhouzVp1auAVcCDdLmlFxERhYYNiJjWSkylP1Ltbh1QLTl9ADiG4sXbd1BMNnVulc5tzy/bt687ue37P1QNNCJiKlvLNFaubc6ACOAHkq7g8epCf0tHPuilymi9tZK+BXzL9n0jizEiIkbDa8WK5c25rWf7/ZLeQDFmQcCA7X+vun+vKTMEfAQ4ruxYktYAn7F96ujCjoiI4bDFmtWNuq2H7cuBy0eyb68rp+MpMt5etu8EkPQs4POS3mv70yM5YEREDJ/XihWP9f9tPUmP0H3QnChmeq/0PlCv5PQW4ADb97dW2L5D0puB/wCSnCIixovF2hX9f1vP9mZj0U+v5DSjPTG1Hfg+STPG4uAREVGRBQ27rSfp5RTVhc6XtDWwWetO3FB6JaeVI9wWERFjzcDyboV3+pOkjwCzgWcD51OUqPsKxeOiIfVKTi+Q9HCX9QI2GmacERExGgZWT3QQw/I6YE/gegDbd0uqfMtv0ORku1nXjxERk9la4LGJDmJYVtq2pNZEtZsMZ+cqL+FGRMREM7BiooMYlq9L+gKwRVlA/O+AL1bdOckpIqIJGnZbz/Yny5J3D1M8dzrZ9g+r7p/kFBHRBGuB5RMdxNAkfZZilvNflMmockJql+QUEdEEphHJiWJ2iU9J2gb4GnCJ7UXD7aRyhdiIiJhAppi/YaSf8QrTPsv2y4B9gQeA8yUtkXSypN2q9pPkFBHRBGspBkSM9DPObP/W9hm296SYy+91wJKq+yc5RUQ0QWtAxEg/Q5B0nqR7Jd08yPYtJf27pJsk/UrS84bob4akgyVdDHwf+C/gDVVOFfLMKSKiGep/5nQB8FngokG2nwgssv06Sc8Bzgb272xUjtA7gmJ2818BlwJzbP/vcIJJcoqIaIKaR+vZvlrSzB5NdgdOK9v+WtJMSc+w/f862p0IfBX4R9sPjDSeJKeIiCZoDYiYODcCrwd+JmlvYCdge+AJycn2fmNxsCSniIgmGH2FiK0lLWxbHrA9MIz9TwfOkrQIWAzcQI2vBSc5RUQ0wehv691ve/ZId7b9MHA0rJsp/c7yU4skp4iIJpjg23qStgAetb0SeBtwdZmwapHkFBHRBDUXfpV0CfAKitt/y4CPADMAbJ8DPBe4SNIa4FbgmPqiqTk5SToQOAuYDpxr+/RB2h0GfAPYy/bCbm0iIqa0mgu/2j5iiO2/BGbVF8ET1ZacJE2nGAd/ALAMWCBpnu1bO9ptBrwHuLauWCIiGq8hhV/HSp0VIvYGltq+o7xHeSlwaJd2HwU+wZT6Y4+IGKbWS7gj/TRMnclpO+CutuVl5bp1JO0J7GD7uzXGERHRfA0p/DpW6nzmpC7rvG6jNA34NHDUkB1Jc4A5AJvvuPkYhRcR0SCtwq9TRJ3JaRmwQ9vy9sDdbcubAc8DflIMmeeZwDxJh3QOiihfFBsA2Hb2NiYiYqpZCzw20UGMnzqT0wJglqSdgd8Dh1OUTQfA9kPA1q1lST+hqMWU0XoREd2smegAxk9tycn2aknHAVdQDCU/z/Ytkk4FFtqeV9exIyImnSk2Wq/W95xszwfmd6w7eZC2r6gzloiIRpv4wq/jKhUiIiKaoOYKEf0mySkiogkyICIiIvqOyYCIiIjoM/VP095XkpwiIpogAyIiIqLvpEJERET0nQyIiIiIvpQBERER0XemUGXROqfMiIiIGJFcOUVENMLUGq6X5BQR0QhTa0REklNERCMYWD3RQYybJKeIiEZYCzw60UGMm+lz586d6BiGZeCUU+a+7tlb8cAuW8I0oTVr2W3+Up7/tSU86ZEV69YDg24b7vp+7atp8fZrX02Lt1/7alq8TehrwSW3csLcuWcCnHLK5+bCGyiunkby+cw9c+fOHRjtb/C4sd2oz4vByzed4aX77+STV7zfS/ffycs3neE1enz9h1f/sz+8+p+7bhtsn6b11bR4+7WvpsXbr301Ld6m9PVicOu3D/7McOsoPiyc6N/v4XwaeeX07pVr2eS+R5m2ag3Pu+w2NvzTKgRsUK7/wwuezlOXPsjLzlq43rbB9mlaX02Lt1/7alq8/dpX0+JtSl8DwJy5c08BOOWUz86FQxj5ldNAo66cGvvMacb/rmKnn/+eGf+7ar312yy6F9ldtw22T9P6alq8/dpX0+Lt176aFm/T+ipMrdF6jX0Jd9UmM/jtPtuxapMZ662/54VP5+49n9F122D7NK2vpsXbr301Ld5+7atp8Tatr0LrPaeRfpqlkclpxaYzWPaSbbnqpH1Y9pJtWbHpDNbq8fW3H7QLtx+0S9dtg+3TtL6aFm+/9tW0ePu1r6bF25S+nqh15TTST7OoeNDWHLttt5k//IUDuf2gXfD0aWjNWmZ9/7/ZZtG93PPCp69bDwy6bbjr+7WvpsXbr301Ld5+7atp8TahrxMPvuzuu+ztAKRZhjNH8et5yHW2Z4+ig3HVuOS07extfMzCoyc6jIiI2n1Mp61LKNIuho+PorfDG5WcGjsgIiJiasmAiDEj6UBJt0laKukDXbYfJek+SYvKz9vqjCciorla5YtG+ulN0nmS7pV08yDbnyLpO5JulHSLpFpvYdV25SRpOnA2cACwDFggaZ7tWzuafs32cXXFERExOdR+5XQB8FngokG2HwvcavtgSU8DbpN0se2VdQRT5229vYGltu8AkHQpcCjQmZwiImJI9SYn21dLmtmrCbCZJAGbAg9QYyXaOm/rbQfc1ba8rFzX6Q2SbpJ0maQdaownIqLBJvw9p88CzwXuBhYD/2B77Vh03E2dV07qsq5zaOB3gEtsr5D0TuBC4JXrdSTNAeaUiys+ptO63hOdArYG7p/oICZIzn1qmsrnDvDsx7/efQWctPUo+tpI0sK25QHbwyln9H+ARRS/0bsAP5T0U9sPjyKmQdWZnJYB7VdC21Nk3HVs/7Ft8YvAGd06Kv8ABwAkLWzScMixlHPPuU81U/ncoTj/1nfbB05kLMDRwOku3j9aKulO4DnAr+o4WJ239RYAsyTtLOlJwOHAvPYGkrZpWzwEWFJjPBERMXK/A/YHkPQMiqu6O+o6WG1XTrZXSzoOuAKYDpxn+xZJp1KUbp8HvEdSq8zuA8BRdcUTERGDk3QJ8Apga0nLgI8AMwBsnwN8FLhA0mKKxzb/bLu2W66NqxAhac4w75NOGjn3nPtUM5XPHab2+TcuOUVExOTXyKrkERExufVtcqpQ+mhDSV8rt187xMtjjVLh3E+QdGv5fth/StppIuKsw1Dn3tbuMEmWNGlGclU5d0l/U/7d3yLpq+MdY10q/JvfUdKVkm4o/92/ZiLirEOFskGS9H/LP5ubJL1ovGOcEBM9T3y3D8UAiv8GngU8CbgR2L2jzbuBc8rvh1OUQZrw2Mfp3PcDnlx+f9dUOvey3WbA1cA1wOyJjnsc/95nATcAW5bLT5/ouMfx3AeAd5Xfdwd+M9Fxj+H5/wXwIuDmQba/Bvg+xSCElwLXTnTM4/Hp1yundaWPXNRtapU+ancoxUu7AJcB+5dlNZpuyHO3faXtR8vFayjeIZsMqvy9QzFq6BPA8vEMrmZVzv3twNm2HwSwfe84x1iXKuduYPPy+1PoeGeyyWxfTTFaeTCHAhe5cA2wRcdrOJNSvyanKqWP1rWxvRp4CHjquERXr6pln1qOofi/qslgyHOXtCewg+3vjmdg46DK3/tuwG6Sfi7pGkkT/VLmWKly7nOBN5dDnOcDfz8+ofWF4f4mTAr9Op9TldJHVdo0UeXzkvRmYDawb60RjZ+e5y5pGvBpJuf7cFX+3jeguLX3Coqr5Z9Kep7t/6k5trpVOfcjgAtsf0rSy4Avl+deW223PjJZf+t66tcrpyFLH7W3kbQBxaV+r0vjpqhy7kh6FfAh4BDbK8YptroNde6bAc8DfiLpNxT33+dNkkERVf/Nf9v2Ktt3ArdRJKumq3LuxwBfB7D9S2Ajirp7U0Gl34TJpl+T05Clj8rlt5bfDwN+7PLpYcNVKfu0J/AFisQ0WZ47wBDnbvsh21vbnml7JsXztkNsL+zeXaNU+Tf/LYrBMEjamuI2X23lY8ZRlXNvL53zXIrkdN+4Rjlx5gFvKUftvRR4yPY9Ex1U3frytp6rlT76EsWl/VKKK6bDJy7isVPx3P+VYj6Vb5RjQH5n+5AJC3qMVDz3SaniuV8BvFrSrcAa4P1+YvHkRqp47u8DvijpvRS3tI6aJP8zWqVs0HyKEXtLgUcpCrBOeqkQERERfadfb+tFRMQUluQUERF9J8kpIiL6TpJTRET0nSSniIjoO0lO0Vck7SDpTklblctblsu1VF6X9E5Jbym/HyVp27Zt50rafYyO81pJJ5ffL5B02Aj7eZqkH4xFTBH9LMkp+ortu4DPA6eXq04HBmz/tqbjnWP7onLxKGDbtm1vs33rGB3qn4DPjbYT2/cB90jaZ/QhRfSvJKfoR58GXirpeODlwKc6G0iaKenXki4s57i5TNKTy237l/P+LC7nytmwXH962zxYnyzXzZX0j+WVzGzgYkmLJG0s6Set0kiSjij7u1nSGW1x/EnSv0i6sSzG+owuse4GrLB9f5dtHy2vpKZJ+o2kj0v6paSFkl4k6QpJ/y3pnW27fQt408j/eCP6X5JT9B3bq4D3UySp48tpFLp5NsVV1R7Aw8C7JW0EXAD8re3nU1RBeVd5m/B1wJ+V7T/WcczLgIXAm2y/0PZjrW3lrb4zgFcCLwT2kvTacvMmwDW2X0Axx9Tbu8S5D3B950pJnwCeDhzdVsD0LtsvA35ansdhFDUET23bdSHw54P8mURMCklO0a8OAu6hKPQ6mLts/7z8/hWKq6xnA3fa/q9y/YUUk7k9TDH/07mSXk9RBqaqvYCf2L6vnJ7l4rJPgJVAa/qO64CZXfbfhvXrwH0Y2ML2OzrK8LRKNC2mmFTukfJW3nJJW5Tb7qXt9mPEZJTkFH1H0guBAyiuGN7bY2K1ztpbpvv0Aq05v/YGLgdeCwxnUEGvSSxXtSWXNXSvV/kYRaHSdguAF7cGfrRpVZhf2/a9tdzqe6Oyz4hJK8kp+oqKSrafp7id9zuKIrefHKT5juXcPlDM9/Mz4NfATEm7luuPBK6StCnwFNvzgeMpbs91eoRiWo5O1wL7Stpa0vTyWFcN47SWALt2rPsBxWCP70nqdsxedgNuHuY+EY2S5BT95u0UVdZ/WC5/DniOpG4TKi4B3irpJmAr4PO2l1NUbf6GpMUUVxznUCSd75ZtrwLe26W/C4BzWgMiWivL6Qk+CFwJ3Ahcb/vbwzinq4E9y8S7ju1vAF+kmJNq4657drcf8L1htI9onFQlj0aSNBP4ru1ez6T6hqSzgO/Y/tEY9HU1cKjtB0cfWUR/ypVTxPj4OPDk0XYi6WnAmUlMMdnlyikiIvpOrpwiIqLvJDlFRETfSXKKiIi+k+QUERF9J8kpIiL6TpJTRET0nf8PbxkbPWWQWfAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Prepare the varying source locations\n",
    "source_locations = src_coordinates\n",
    "\n",
    "plot_velocity(model, source=source_locations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 1 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.25 s\n",
      "Operator `Kernel` run in 0.19 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 2 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 3 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.19 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 4 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 5 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.19 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 6 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.16 s\n",
      "Operator `Kernel` run in 0.21 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 7 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.19 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 8 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 9 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.22 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 10 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.16 s\n",
      "Operator `Kernel` run in 0.21 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 11 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.16 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 12 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 13 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.16 s\n",
      "Operator `Kernel` run in 0.21 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 14 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.21 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 15 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.16 s\n",
      "Operator `Kernel` run in 0.21 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 16 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.16 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 17 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 18 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.21 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 19 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.21 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 20 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.21 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 21 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.21 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 22 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.21 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 23 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 24 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 25 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 26 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 27 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.16 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 28 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 29 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.21 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 30 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 31 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 32 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.16 s\n",
      "Operator `Kernel` run in 0.21 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 33 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.16 s\n",
      "Operator `Kernel` run in 0.21 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 34 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 35 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 36 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 37 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 38 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 39 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 40 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 41 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 42 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.19 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 43 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.16 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 44 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 45 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 46 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.19 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 47 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.19 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 48 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 49 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.19 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 50 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.19 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Imaging source 51 out of 51\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Forward` run in 0.15 s\n",
      "Operator `Kernel` run in 0.20 s\n"
     ]
    }
   ],
   "source": [
    "# Run imaging loop over shots\n",
    "from devito import Function\n",
    "\n",
    "# Create image symbol and instantiate the previously defined imaging operator\n",
    "image = Function(name='image', grid=model.grid)\n",
    "op_imaging = ImagingOperator(model, image)\n",
    "solver = AcousticWaveSolver(model, geometry, space_order=4)\n",
    "objective=0\n",
    "for i in range(nshots):\n",
    "    print('Imaging source %d out of %d' % (i+1, nshots))\n",
    "    \n",
    "    # Update source location\n",
    "    geometry.src_positions[0, :] = source_locations[i, :]\n",
    "\n",
    "    # Generate synthetic data from true model\n",
    "    #true_d, _, _ = solver.forward(vp=model.vp)\n",
    "    \n",
    "    # Compute smooth data and full forward wavefield u0\n",
    "    smooth_d, u0, _ = solver.forward(vp=model.vp, save=True)\n",
    "    dat = smooth_d.data\n",
    "    \n",
    "    # Compute gradient from the data residual  \n",
    "    v = TimeFunction(name='v', grid=model.grid, time_order=2, space_order=4)\n",
    "    residual = dat - names['data' + str(i)]\n",
    "    residual = np.float32(residual)\n",
    "    op_imaging(u=u0, v=v, vp=model.vp, dt=1., \n",
    "               residual=residual)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAADJCAYAAADGvkL8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAPbElEQVR4nO3df6xkd1nH8ffHlqIUtC0spLbVbk1FiYm2u8EqQgxFpRVpVTAlRDfapDEBBdFIkcTwp/UHoImBrBStptJigbQxUSG1aPyDyra0tGUpXUopS5d2QX4FjVh5/GPONdPbuffOnTl3zrnffb+Sm5k58+u533PmM89858xMqgpJUpu+begCJEk7x5CXpIYZ8pLUMENekhpmyEtSwwx5SWrYUiGf5KVJ7k9yJMnVfRUlSepHFt1PPslJwKeAnwKOAh8FXlVVn+ivPEnSMpbp5J8PHKmqB6vqm8ANwGX9lCVJ6sPJS1z3LOBzU6ePAj+62RWS+PFaSdq+L1bVnkWuuEzIZ8ayJ4V4kquAq5a4H0k60X120SsuE/JHgXOmTp8NPLL+QlV1EDgIdvKStGrLzMl/FDg/yd4kpwBXALf0U5YkqQ8Ld/JV9XiS1wL/BJwEvLuq7uutMknS0hbehXKhO3O6RpIWcUdV7V/kin7iVZIaZshLUsMMeekEtNE0bVUtfF6ft6f+GPKS1LBl9pOXtEslsz7LuPHyRc9b9PbUHzt5SWqYIS9JDTPkJalhhrwkNayZN15PhF2xfKNKLWj1sTrWx6edvCQ1bJSdfKvP9MvaTeMy1q6mVbtp22jVIutgFY+TwUPejbNNrldpa7MeJ30Hv9M1ktSwQTp5uzxJmm06H/vo6u3kJalhKw35ffv22cVL0pz6+KZOO3lJapghL0kNM+QlqWGGvCQ1zJCXpIYZ8pLUMENekhpmyEtSw7YM+STnJLktyeEk9yV5Xbf8jCQfSvJAd3j6zpcrSdqOeTr5x4HfrqofBC4CXpPkecDVwK1VdT5wa3dakjQiW4Z8VR2rqju7418HDgNnAZcB13UXuw64fKeKlCQtZltz8knOBS4AbgeeU1XHYPJEADx7g+tcleRQkkPHjx9frlpJ0rbMHfJJng68D3h9VX1t3utV1cGq2l9V+/fs2bNIjZKkBc0V8kmewiTgr6+q93eLH01yZnf+mcBjO1OiJGlR8+xdE+Ba4HBVvXXqrFuAA93xA8DN/ZcnSVrGPL8M9QLgl4F7ktzVLfs94A+A9ya5EngYeOXOlChJWtSWIV9V/wZs9BtUF/dbjiSpT37iVZIaZshLUsMMeUlqmCEvSQ0z5CWpYYa8JDXMkJekhs3zYaiVm3zIdqKqBqxEkp5sOqP61nfmDR7yswZr7Z9MstRg+gQhaTt2MrwXqaGPDHO6RpIaNkgnv9WzZV/PpvPejh2/dOIYQ7c+rz5qtZOXpIZllV1skrnubHpOfizs9qXdawxZsmSu3VFV+xe54uBvvM4yhhWynnv8SLvH2DNklZyukaSGjTLkq2rU3fKyu3ZKOvEMlWujDHlJUj+ck5fUlLHmh3PyU8Y+XTP2+qQT2Vgfm07XSJJ6t9KQ37dv31zPZEO9sbn2TLvVn6RxG+Njdqhcs5OXpIYN8sbrWJ5ZJZ0YZmXOqrrqIe8bRrZ3zfp/3CcDSTtlVfmy2depL3Mb85p7uibJSUk+luTvu9N7k9ye5IEkNyY5ZeEqJEk7Yjtz8q8DDk+dvgZ4W1WdD3wZuHLZYsb4ZokkLWPonTfmCvkkZwM/C7yrOx3gxcBN3UWuAy7fiQIlSYubt5N/O/C7wLe6088EvlJVj3enjwJn9VybJDVjtLtQJnkZ8FhV3TG9eMZFZ74GSXJVkkNJDh0/fnzBMiVpdxvzJ15fALw8yUPADUymad4OnJZkbe+cs4FHZl25qg5W1f6q2r9nz54eSpYkzWvLkK+qN1XV2VV1LnAF8M9V9WrgNuAV3cUOADfvWJWSpIUs84nXNwJvSHKEyRz9tf2UJEnqy7Y+DFVVHwY+3B1/EHh+/yVJkvrid9dIUsMMeUlqmCEvSQ0z5CWpYYa8JDXMkJekhhnyktQwQ16SGmbIS1LDRhnyQ30lpyTtlKFybVS/8brGX4WS1Jqhcm2UnbwkqR+jDHmnayS1ZrS/DCVJ2r2ck5ekFXBOXpLUu1F18mOZh/eVhKQ1q8ilncycQUJ+q0Fb+4eHCv2duF+fOKT+jKUh7MtG/08fueF0jSQ1LKvsMJPYzq7YTq/f1jqqE9ky24rbwY67o6r2L3JFO3lJatio3nhV/+ywNC+3lTbZyUtSwwx5SWrYXCGf5LQkNyX5ZJLDSX4syRlJPpTkge7w9J0uVpK0PfN28n8K/GNV/QDww8Bh4Grg1qo6H7i1Oy1JGpEtd6FM8p3A3cB5NXXhJPcDP1lVx5KcCXy4qp67xW25C6Ukbd+O7kJ5HnAc+MskH0vyriSnAs+pqmMA3eGzFylAkrRz5gn5k4ELgXdU1QXAN9jG1EySq5IcSnJowRolSQuaJ+SPAker6vbu9E1MQv/RbpqG7vCxWVeuqoNVtX/RlxqSpMVtGfJV9QXgc0nW5tsvBj4B3AIc6JYdAG7ekQolSQub9xOvvwFcn+QU4EHgV5k8Qbw3yZXAw8Ard6ZESdKiRvkFZUN/1bAk9W3JXPMLyiRJTzbKLyizg5fUmqFyzU5ekho2SCc/pp/C81WDpM2MIa+WyamVdvL79u2ba8CSrCx8q2rTP0ntW0UGrDLXpjldI0kNG+Ubr2PqoGfV4hSPtHsNlS9D3a+dvCQ1bJSd/FqnPKaOftp0XXb10u4wdJ4MlWujDPmhV8Z2+OlcafzGkClO10iSejfKTn7s0zWSdocxZchQuWYnL0kNG2UnP6Zn3604Fz9uQ21LbhfjMKZZgaFqGGXIj2nFzOIDeH5jXYc7bUz/t9vrsDYb/1VsJ07XSFLDRtPJz3q2swMZvzF1rJpto3V0Ij2+pv/XvrfZZcZx1nX7rs9OXpIaNkgnfyJ1EC2xa2/LifpBvrH/v32/6hjlb7xqHAz1E8vYw+8E52+8SpKezJDXTHbxJx7XeZsMeUlq2Gh2odQ42M1JbZmrk0/yW0nuS3Jvkvck+fYke5PcnuSBJDcmOWWni5W0c3zjtU1bhnySs4DfBPZX1Q8BJwFXANcAb6uq84EvA1fuZKGSpO2bd07+ZOA7kpwMPA04BrwYuKk7/zrg8v7L06qt/aK8XV37pte167tdW4Z8VX0e+GPgYSbh/lXgDuArVfV4d7GjwFk7VaQkaTHzTNecDlwG7AW+GzgVuGTGRWe+Y5fkqiSHkhxaplCt3vpOb6M/jYvrTdPm2bvmJcBnquo4QJL3Az8OnJbk5K6bPxt4ZNaVq+ogcLC7rrtuNMjAkMZrnjn5h4GLkjwtk0fzxcAngNuAV3SXOQDcvDMlSpIWNc+c/O1M3mC9E7inu85B4I3AG5IcAZ4JXLuDdUqSFuAXlEnS+PkFZZKkJzPkJalhhrwkNcyQl6SGGfKS1DBDXpIaZshLUsMMeUlqmCEvSQ0z5CWpYYa8JDXMkJekhhnyktQwQ16SVqCqWOW3/q4x5CWpYfP8/J8kaQ5rnfqsn8Qc6mcyV9rJ79u37/9fsqz/AzZcvuh5Q93eLMuc5+15e97eYrc3hjxY5q8PTtdIUsNGM12z3WfqRc9bxe3txpq9PW+vxdtbdS19W7ufZaZ67OQlqWGGvCQ1zJCXpIYZ8pLUMENekhpmyEtSwwx5SWqYIS9JDTPkJalhWdUntwCSHAe+AXxxZXe6mGcx/hrBOvu0G2oE6+zTbqgRJnWeWlV7FrnySkMeIMmhqtq/0jvdpt1QI1hnn3ZDjWCdfdoNNcLydTpdI0kNM+QlqWFDhPzBAe5zu3ZDjWCdfdoNNYJ19mk31AhL1rnyOXlJ0uo4XSNJDVtZyCd5aZL7kxxJcvWq7ncrSc5JcluSw0nuS/K6bvlbknw+yV3d36UD1/lQknu6Wg51y85I8qEkD3SHpw9c43OnxuuuJF9L8voxjGWSdyd5LMm9U8tmjl8m/qzbVj+e5MIBa/yjJJ/s6vhAktO65ecm+a+pMX3nKmrcpM4N13GSN3VjeX+Snxm4zhunanwoyV3d8kHGc5P86W/b7Ou3CLf4ncKTgE8D5wGnAHcDz1vFfc9R25nAhd3xZwCfAp4HvAX4naHrm6rzIeBZ65b9IXB1d/xq4Jqh61y3zr8AfO8YxhJ4EXAhcO9W4wdcCvwDEOAi4PYBa/xp4OTu+DVTNZ47fbkRjOXMddw9lu4Gngrs7XLgpKHqXHf+nwC/P+R4bpI/vW2bq+rknw8cqaoHq+qbwA3AZSu6701V1bGqurM7/nXgMHDWsFXN7TLguu74dcDlA9ay3sXAp6vqs0MXAlBV/wr8x7rFG43fZcBf18RHgNOSnDlEjVX1wap6vDv5EeDsna5jKxuM5UYuA26oqv+uqs8AR5jkwY7brM5Mfk/vl4D3rKKWjWySP71tm6sK+bOAz02dPsoIgzTJucAFwO3dotd2L4nePfRUCFDAB5PckeSqbtlzquoYTDYW4NmDVfdkV/DEB9CYxnLNRuM31u3115h0cWv2JvlYkn9J8sKhipoyax2PdSxfCDxaVQ9MLRt0PNflT2/b5qpCftav0I5qt54kTwfeB7y+qr4GvAP4PuBHgGNMXtoN6QVVdSFwCfCaJC8auJ4NJTkFeDnwd92isY3lVka3vSZ5M/A4cH236BjwPVV1AfAG4G+TfOdQ9bHxOh7dWHZexRObkEHHc0b+bHjRGcs2Hc9VhfxR4Jyp02cDj6zovreU5ClMBvj6qno/QFU9WlX/W1XfAv6CFb3E3EhVPdIdPgZ8oKvn0bWXat3hY8NV+ASXAHdW1aMwvrGcstH4jWp7TXIAeBnw6uomZrvpjy91x+9gMtf9/UPVuMk6HtVYAiQ5GfgF4Ma1ZUOO56z8ocdtc1Uh/1Hg/CR7uy7vCuCWFd33prq5uWuBw1X11qnl0/NcPw/cu/66q5Lk1CTPWDvO5M24e5mM4YHuYgeAm4ep8Eme0CWNaSzX2Wj8bgF+pduT4SLgq2svnVctyUuBNwIvr6r/nFq+J8lJ3fHzgPOBB4eosatho3V8C3BFkqcm2cukzn9fdX3rvAT4ZFUdXVsw1HhulD/0uW2u8F3kS5m8c/xp4M2rut856voJJi93Pg7c1f1dCvwNcE+3/BbgzAFrPI/JHgp3A/etjR/wTOBW4IHu8IwRjOfTgC8B3zW1bPCxZPKkcwz4Hybd0JUbjR+Tl8R/3m2r9wD7B6zxCJM52LVt853dZX+x2xbuBu4Efm7gsdxwHQNv7sbyfuCSIevslv8V8OvrLjvIeG6SP71tm37iVZIa5ideJalhhrwkNcyQl6SGGfKS1DBDXpIaZshLUsMMeUlqmCEvSQ37PzY4ejGslTw+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from examples.seismic import plot_image\n",
    "\n",
    "# Plot the inverted image\n",
    "data0 = image.data[20:-20,20:-20]\n",
    "plt.imshow(np.diff(data0, axis=1).T,vmin=-1,vmax=10,cmap='gray')\n",
    "plt.savefig('rtm4sr')\n",
    "#plot_image(image.data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define bounding box constraints on the solution.\n",
    "def apply_box_constraint(vp):\n",
    "    # Maximum possible 'realistic' velocity is 3.5 km/sec\n",
    "    # Minimum possible 'realistic' velocity is 2 km/sec\n",
    "    return np.clip(vp, 2.0, 3.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#NBVAL_SKIP\n",
    "\n",
    "# Run FWI with gradient descent\n",
    "history = np.zeros((fwi_iterations, 1))\n",
    "for i in range(0, fwi_iterations):\n",
    "    # Compute the functional value and gradient for the current\n",
    "    # model estimate\n",
    "    phi, direction = fwi_gradient(model.vp)\n",
    "    \n",
    "    # Store the history of the functional values\n",
    "    history[i] = phi\n",
    "    \n",
    "    # Artificial Step length for gradient descent\n",
    "    # In practice this would be replaced by a Linesearch (Wolfe, ...)\n",
    "    # that would guarantee functional decrease Phi(m-alpha g) <= epsilon Phi(m)\n",
    "    # where epsilon is a minimum decrease constant\n",
    "    alpha = .05 / np.abs(direction).max()\n",
    "    \n",
    "    # Update the model estimate and enforce minimum/maximum values\n",
    "    model.vp = apply_box_constraint(model.vp.data - alpha * direction)\n",
    "    \n",
    "    # Log the progress made\n",
    "    print('Objective value is %f at iteration %d' % (phi, i+1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}