{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# AVO modelling\nThis example shows how to create pre-stack angle gathers using\nthe :py:class:`pylops.avo.avo.AVOLinearModelling` 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\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, np.random.normal(0, 80, nt0))\nvs = 600 + vp/2 + filtfilt(np.ones(5)/5., 1, np.random.normal(0, 20, nt0))\nrho = 1000 + vp + filtfilt(np.ones(5)/5., 1, np.random.normal(0, 30, nt0))\nvp[201:] += 500\nvs[201:] += 200\nrho[201:] += 100\n\n# Wavelet\nntwav = 41\nwavoff = 10\nwav, twav, wavc = ricker(t0[:ntwav//2+1], 20)\nwav_phase = np.hstack((wav[wavoff:], np.zeros(wavoff)))\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=(9, 7), sharey=True)\naxs[0].plot(m[:, 0], t0, 'k', lw=6)\naxs[0].set_title('Vp')\naxs[0].set_ylabel(r'$t(s)$')\naxs[0].invert_yaxis()\naxs[0].grid()\naxs[1].plot(m[:, 1], t0, 'k', lw=6)\naxs[1].set_title('Vs')\naxs[1].invert_yaxis()\naxs[1].grid()\naxs[2].plot(m[:, 2], t0, 'k', lw=6, label='true')\naxs[2].set_title('Rho')\naxs[2].invert_yaxis()\naxs[2].grid()\naxs[2].legend()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We create now the operators to model the AVO responses for a set of\nelastic profiles\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# constant vsvp\nPPop_const = \\\n    pylops.avo.avo.AVOLinearModelling(theta, vsvp=vsvp,\n                                      nt0=nt0, linearization='akirich',\n                                      dtype=np.float64)\n\n# depth-variant vsvp\nPPop_variant = \\\n    pylops.avo.avo.AVOLinearModelling(theta, vsvp=vsvp_z,\n                                      linearization='akirich',\n                                      dtype=np.float64)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can then apply those operators to the elastic model and\ncreate some synthetic reflection responses\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dPP_const = PPop_const * m.ravel()\ndPP_const = dPP_const.reshape(nt0, ntheta)\n\ndPP_variant = PPop_variant * m.ravel()\ndPP_variant = dPP_variant.reshape(nt0, ntheta)\n\nfig, axs = plt.subplots(1, 2, figsize=(10, 5), sharey=True)\naxs[0].imshow(dPP_const, cmap='gray',\n              extent=(theta[0], theta[-1], t0[-1], t0[0]),\n              vmin=dPP_const.min(), vmax=dPP_const.max())\naxs[0].set_title('Data with constant VP/VS')\naxs[0].axis('tight')\naxs[1].imshow(dPP_variant, cmap='gray',\n              extent=(theta[0], theta[-1], t0[-1], t0[0]),\n              vmin=dPP_variant.min(), vmax=dPP_variant.max())\naxs[1].set_title('Data with variable VP/VS')\naxs[1].axis('tight')\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally we can also model the PS response by simply changing the\n``linearization`` choice as follows\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "PSop = \\\n    pylops.avo.avo.AVOLinearModelling(theta, vsvp=vsvp,\n                                      nt0=nt0, linearization='ps',\n                                      dtype=np.float64)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can then apply those operators to the elastic model and\ncreate some synthetic reflection responses\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dPS = PSop * m.ravel()\ndPS = dPS.reshape(nt0, ntheta)\n\nfig, axs = plt.subplots(1, 2, figsize=(10, 5), sharey=True)\naxs[0].imshow(dPP_const, cmap='gray',\n              extent=(theta[0], theta[-1], t0[-1], t0[0]),\n              vmin=dPP_const.min(), vmax=dPP_const.max())\naxs[0].set_title('PP Data')\naxs[0].axis('tight')\naxs[1].imshow(dPS, cmap='gray',\n              extent=(theta[0], theta[-1], t0[-1], t0[0]),\n              vmin=dPS.min(), vmax=dPS.max())\naxs[1].set_title('PS Data')\naxs[1].axis('tight')\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
}