{
  "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"
      ]
    },
    {
      "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\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.0, 40.0\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 = 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 = 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 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 = pylops.avo.avo.AVOLinearModelling(\n    theta, vsvp=vsvp, nt0=nt0, linearization=\"akirich\", dtype=np.float64\n)\n\n# depth-variant vsvp\nPPop_variant = pylops.avo.avo.AVOLinearModelling(\n    theta, vsvp=vsvp_z, linearization=\"akirich\", dtype=np.float64\n)"
      ]
    },
    {
      "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\ndPP_variant = PPop_variant * m"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To visualize these responses, we will plot their anomaly - how much they\ndeveiate from their mean\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mean_dPP_const = dPP_const.mean()\ndPP_const -= mean_dPP_const\nmean_dPP_variant = dPP_variant.mean()\ndPP_variant -= mean_dPP_variant\n\nfig, axs = plt.subplots(1, 2, figsize=(10, 5), sharey=True)\nim = axs[0].imshow(\n    dPP_const,\n    cmap=\"RdBu_r\",\n    extent=(theta[0], theta[-1], t0[-1], t0[0]),\n    vmin=-dPP_const.max(),\n    vmax=dPP_const.max(),\n)\ncax = make_axes_locatable(axs[0]).append_axes(\"right\", size=\"5%\", pad=\"2%\")\ncb = fig.colorbar(im, cax=cax)\ncb.set_label(f\"Deviation from mean = {mean_dPP_const:.2f}\")\naxs[0].set(xlabel=r\"$\\theta$ [\u00b0]\", ylabel=r\"$t$ [s]\", title=\"Data with constant VP/VS\")\naxs[0].axis(\"tight\")\n\nim = axs[1].imshow(\n    dPP_variant,\n    cmap=\"RdBu_r\",\n    extent=(theta[0], theta[-1], t0[-1], t0[0]),\n    vmin=-dPP_variant.max(),\n    vmax=dPP_variant.max(),\n)\ncax = make_axes_locatable(axs[1]).append_axes(\"right\", size=\"5%\", pad=\"2%\")\ncb = fig.colorbar(im, cax=cax)\ncb.set_label(f\"Deviation from mean = {mean_dPP_variant:.2f}\")\naxs[1].set(xlabel=r\"$\\theta$ [\u00b0]\", 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 = pylops.avo.avo.AVOLinearModelling(\n    theta, vsvp=vsvp, nt0=nt0, linearization=\"ps\", dtype=np.float64\n)"
      ]
    },
    {
      "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\nmean_dPS = dPS.mean()\ndPS -= mean_dPS\n\nfig, axs = plt.subplots(1, 2, figsize=(10, 5), sharey=True)\nim = axs[0].imshow(\n    dPP_const,\n    cmap=\"RdBu_r\",\n    extent=(theta[0], theta[-1], t0[-1], t0[0]),\n    vmin=-dPP_const.max(),\n    vmax=dPP_const.max(),\n)\ncax = make_axes_locatable(axs[0]).append_axes(\"right\", size=\"5%\", pad=\"2%\")\ncb = fig.colorbar(im, cax=cax)\ncb.set_label(f\"Deviation from mean = {mean_dPP_const:.2f}\")\naxs[0].set(xlabel=r\"$\\theta$ [\u00b0]\", ylabel=r\"$t$ [s]\", title=\"PP Data\")\naxs[0].axis(\"tight\")\n\nim = axs[1].imshow(\n    dPS,\n    cmap=\"RdBu_r\",\n    extent=(theta[0], theta[-1], t0[-1], t0[0]),\n    vmin=-dPS.max(),\n    vmax=dPS.max(),\n)\ncax = make_axes_locatable(axs[1]).append_axes(\"right\", size=\"5%\", pad=\"2%\")\ncb = fig.colorbar(im, cax=cax)\ncb.set_label(f\"Deviation from mean = {mean_dPS:.2f}\")\naxs[1].set(xlabel=r\"$\\theta$ [\u00b0]\", 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.9.15"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}