{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Pre-stack modelling\nThis example shows how to create pre-stack angle gathers using\nthe :py:class:`pylops.avo.prestack.PrestackLinearModelling` operator.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport matplotlib.pyplot as plt\nfrom scipy.signal import filtfilt\n\nimport pylops\nfrom pylops.utils.wavelets import ricker\n\nplt.close('all')\nnp.random.seed(0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's start by creating the input elastic property profiles and wavelet\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nt0 = 501\ndt0 = 0.004\nntheta = 21\n\nt0 = np.arange(nt0)*dt0\nthetamin, thetamax = 0, 40\ntheta = np.linspace(thetamin, thetamax, ntheta)\n\n# Elastic property profiles\nvp = 1200 + np.arange(nt0) + filtfilt(np.ones(5)/5., 1,\n                                      np.random.normal(0, 160, nt0))\nvs = 600 + vp/2 + filtfilt(np.ones(5)/5., 1, np.random.normal(0, 100, nt0))\nrho = 1000 + vp + filtfilt(np.ones(5)/5., 1, np.random.normal(0, 120, nt0))\nvp[201:] += 500\nvs[201:] += 200\nrho[201:] += 100\n\n# Wavelet\nntwav = 81\nwav, twav, wavc = ricker(t0[:ntwav//2+1], 5)\n\n# vs/vp profile\nvsvp = 0.5\nvsvp_z = np.linspace(0.4, 0.6, nt0)\n\n# Model\nm = np.stack((np.log(vp), np.log(vs), np.log(rho)), axis=1)\n\nfig, axs = plt.subplots(1, 3, figsize=(13, 7), sharey=True)\naxs[0].plot(vp, t0, 'k')\naxs[0].set_title('Vp')\naxs[0].set_ylabel(r'$t(s)$')\naxs[0].invert_yaxis()\naxs[0].grid()\naxs[1].plot(vs, t0, 'k')\naxs[1].set_title('Vs')\naxs[1].invert_yaxis()\naxs[1].grid()\naxs[2].plot(rho, t0, 'k')\naxs[2].set_title('Rho')\naxs[2].invert_yaxis()\naxs[2].grid()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We create now the operators to model a synthetic pre-stack seismic gather\nwith a zero-phase using both a constant and a depth-variant ``vsvp`` profile\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# constant vsvp\nPPop_const = \\\n    pylops.avo.prestack.PrestackLinearModelling(wav, theta, vsvp=vsvp, nt0=nt0,\n                                                linearization='akirich')\n\n# depth-variant vsvp\nPPop_variant = \\\n    pylops.avo.prestack.PrestackLinearModelling(wav, theta, vsvp=vsvp_z,\n                                                linearization='akirich')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's apply those operators to the elastic model and create some\nsynthetic data\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dPP_const = PPop_const *m.flatten()\ndPP_const = dPP_const.reshape(nt0, ntheta)\n\ndPP_variant = PPop_variant *m.flatten()\ndPP_variant = dPP_variant.reshape(nt0, ntheta)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally we visualize the two datasets\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# sphinx_gallery_thumbnail_number = 2\nfig = plt.figure(figsize=(6, 7))\nax1 = plt.subplot2grid((3, 2), (0, 0), rowspan=2)\nax2 = plt.subplot2grid((3, 2), (0, 1), rowspan=2)\nax3 = plt.subplot2grid((3, 2), (2, 0))\nax4 = plt.subplot2grid((3, 2), (2, 1))\nax1.imshow(dPP_const, cmap='bwr',\n           extent=(theta[0], theta[-1], t0[-1], t0[0]),\n           vmin=-0.1, vmax=0.1)\nax1.set_xlabel(r'$\\Theta$')\nax1.set_ylabel(r'$t(s)$')\nax1.set_title(r'Data with constant $VP/VS$', fontsize=10)\nax1.axis('tight')\nax2.imshow(dPP_variant, cmap='bwr',\n           extent=(theta[0], theta[-1], t0[-1], t0[0]),\n           vmin=-0.1, vmax=0.1)\nax2.set_title(r'Data with depth-variant $VP/VS$', fontsize=10)\nax2.set_xlabel(r'$\\Theta$')\nax2.axis('tight')\nax3.plot(theta, dPP_const[nt0//4], 'k', lw=2)\nax3.plot(theta, dPP_variant[nt0//4], '--r', lw=2)\nax3.set_title('AVO curve at t=%.2f s' % t0[nt0//4], fontsize=10)\nax3.set_xlabel(r'$\\Theta$')\nax4.plot(theta, dPP_const[nt0//2], 'k', lw=2, label=r'constant $VP/VS$')\nax4.plot(theta, dPP_variant[nt0//2], '--r', lw=2, label=r'variable $VP/VS$')\nax4.set_title('AVO curve at t=%.2f s' % t0[nt0//2], fontsize=10)\nax4.set_xlabel(r'$\\Theta$')\nax4.legend()\nplt.tight_layout()"
      ]
    }
  ],
  "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.6.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}