{
  "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"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\nfrom mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable\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 = (\n    2000\n    + 5 * np.arange(nt0)\n    + 2 * filtfilt(np.ones(5) / 5.0, 1, np.random.normal(0, 160, nt0))\n)\nvs = 600 + vp / 2 + 3 * filtfilt(np.ones(5) / 5.0, 1, np.random.normal(0, 100, nt0))\nrho = 1000 + vp + filtfilt(np.ones(5) / 5.0, 1, np.random.normal(0, 120, nt0))\nvp[201:] += 1500\nvs[201:] += 500\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 = vs / vp\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(vp, t0, \"k\", lw=3)\naxs[0].set(xlabel=\"[m/s]\", ylabel=r\"$t$ [s]\", ylim=[t0[0], t0[-1]], title=\"Vp\")\naxs[0].grid()\naxs[1].plot(vp / vs, t0, \"k\", lw=3)\naxs[1].set(title=\"Vp/Vs\")\naxs[1].grid()\naxs[2].plot(rho, t0, \"k\", lw=3)\naxs[2].set(xlabel=\"[kg/m\u00b3]\", 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 = pylops.avo.prestack.PrestackLinearModelling(\n    wav, theta, vsvp=vsvp, nt0=nt0, linearization=\"akirich\"\n)\n\n# depth-variant vsvp\nPPop_variant = pylops.avo.prestack.PrestackLinearModelling(\n    wav, theta, vsvp=vsvp_z, linearization=\"akirich\"\n)"
      ]
    },
    {
      "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\ndPP_variant = PPop_variant * m"
      ]
    },
    {
      "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, sharey=ax1)\nax3 = plt.subplot2grid((3, 2), (2, 0), sharex=ax1)\nax4 = plt.subplot2grid((3, 2), (2, 1), sharex=ax2)\nim = ax1.imshow(\n    dPP_const,\n    cmap=\"bwr\",\n    extent=(theta[0], theta[-1], t0[-1], t0[0]),\n    vmin=-0.2,\n    vmax=0.2,\n)\ncax = make_axes_locatable(ax1).append_axes(\"bottom\", size=\"5%\", pad=\"3%\")\ncb = fig.colorbar(im, cax=cax, orientation=\"horizontal\")\ncb.ax.xaxis.set_ticks_position(\"bottom\")\nax1.set(ylabel=r\"$t$ [s]\")\nax1.set_title(r\"Data with constant $VP/VS$\", fontsize=10)\nax1.tick_params(labelbottom=False)\nax1.axhline(t0[nt0 // 4], color=\"k\", linestyle=\"--\")\nax1.axhline(t0[nt0 // 2], color=\"k\", linestyle=\"--\")\nax1.axis(\"tight\")\n\nim = ax2.imshow(\n    dPP_variant,\n    cmap=\"bwr\",\n    extent=(theta[0], theta[-1], t0[-1], t0[0]),\n    vmin=-0.2,\n    vmax=0.2,\n)\ncax = make_axes_locatable(ax2).append_axes(\"bottom\", size=\"5%\", pad=\"3%\")\ncb = fig.colorbar(im, cax=cax, orientation=\"horizontal\")\ncb.ax.xaxis.set_ticks_position(\"bottom\")\nax2.set_title(r\"Data with depth-variant $VP/VS$\", fontsize=10)\nax2.tick_params(labelbottom=False, labelleft=False)\nax2.axhline(t0[nt0 // 4], color=\"k\", linestyle=\"--\")\nax2.axhline(t0[nt0 // 2], color=\"k\", linestyle=\"--\")\nax2.axis(\"tight\")\n\nax3.plot(theta, dPP_const[nt0 // 4], \"k\", lw=2)\nax3.plot(theta, dPP_variant[nt0 // 4], \"--r\", lw=2)\nax3.set(xlabel=r\"$\\theta$ [\u00b0]\")\nax3.set_title(\"AVO curve at t=%.2f s\" % t0[nt0 // 4], fontsize=10)\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(xlabel=r\"$\\theta$ [\u00b0]\")\nax4.set_title(\"AVO curve at t=%.2f s\" % t0[nt0 // 2], fontsize=10)\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.9.15"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}