{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Chirp Radon Transform\nThis example shows how to use the :py:class:`pylops.signalprocessing.ChirpRadon2D`\nand :py:class:`pylops.signalprocessing.ChirpRadon3D` operators to apply the\nlinear Radon Transform to 2-dimensional or 3-dimensional signals, respectively.\n\nWhen working with the linear Radon transform, this is a faster implementation\ncompared to in :py:class:`pylops.signalprocessing.Radon2D` and\n:py:class:`pylops.signalprocessing.Radon3D` and should be preferred.\nThis method provides also an analytical inverse.\n\nNote that the forward and adjoint definitions in these two pairs of operators\nare swapped.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\nimport pylops\n\nplt.close(\"all\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's start by creating a empty 2d matrix of size $n_x \\times n_t$\nwith a single linear event.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "par = {\n    \"ot\": 0,\n    \"dt\": 0.004,\n    \"nt\": 51,\n    \"ox\": -250,\n    \"dx\": 10,\n    \"nx\": 51,\n    \"oy\": -250,\n    \"dy\": 10,\n    \"ny\": 51,\n    \"f0\": 40,\n}\ntheta = [0.0]\nt0 = [0.1]\namp = [1.0]\n\n# Create axes\nt, t2, x, y = pylops.utils.seismicevents.makeaxis(par)\ndt, dx, dy = par[\"dt\"], par[\"dx\"], par[\"dy\"]\n\n# Create wavelet\nwav, _, wav_c = pylops.utils.wavelets.ricker(t[:41], f0=par[\"f0\"])\n\n# Generate data\n_, d = pylops.utils.seismicevents.linear2d(x, t, 1500.0, t0, theta, amp, wav)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can now define our operators and apply the forward, adjoint and inverse\nsteps.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "npx, pxmax = par[\"nx\"], 5e-4\npx = np.linspace(-pxmax, pxmax, npx)\n\nR2Op = pylops.signalprocessing.ChirpRadon2D(t, x, pxmax * dx / dt, dtype=\"float64\")\ndL_chirp = R2Op * d\ndadj_chirp = R2Op.H * dL_chirp\ndinv_chirp = R2Op.inverse(dL_chirp).reshape(R2Op.dimsd)\n\nfig, axs = plt.subplots(1, 4, figsize=(12, 4), sharey=True)\naxs[0].imshow(d.T, vmin=-1, vmax=1, cmap=\"bwr_r\", extent=(x[0], x[-1], t[-1], t[0]))\naxs[0].set(xlabel=r\"$x$ [m]\", ylabel=r\"$t$ [s]\", title=\"Input model\")\naxs[0].axis(\"tight\")\naxs[1].imshow(\n    dL_chirp.T,\n    cmap=\"bwr_r\",\n    vmin=-dL_chirp.max(),\n    vmax=dL_chirp.max(),\n    extent=(1e3 * px[0], 1e3 * px[-1], t[-1], t[0]),\n)\naxs[1].set(xlabel=r\"$p$ [s/km]\", title=\"Radon Chirp\")\naxs[1].axis(\"tight\")\naxs[2].imshow(\n    dadj_chirp.T,\n    cmap=\"bwr_r\",\n    vmin=-dadj_chirp.max(),\n    vmax=dadj_chirp.max(),\n    extent=(x[0], x[-1], t[-1], t[0]),\n)\naxs[2].set(xlabel=r\"$x$ [m]\", title=\"Adj Radon Chirp\")\naxs[2].axis(\"tight\")\naxs[3].imshow(\n    dinv_chirp.T,\n    cmap=\"bwr_r\",\n    vmin=-d.max(),\n    vmax=d.max(),\n    extent=(x[0], x[-1], t[-1], t[0]),\n)\naxs[3].set(xlabel=r\"$x$ [m]\", title=\"Inv Radon Chirp\")\naxs[3].axis(\"tight\")\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally we repeat the same exercise with 3d data.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "par = {\n    \"ot\": 0,\n    \"dt\": 0.004,\n    \"nt\": 51,\n    \"ox\": -400,\n    \"dx\": 10,\n    \"nx\": 81,\n    \"oy\": -600,\n    \"dy\": 10,\n    \"ny\": 61,\n    \"f0\": 20,\n}\ntheta = [10]\nphi = [0]\nt0 = [0.1]\namp = [1.0]\n\n# Create axes\nt, t2, x, y = pylops.utils.seismicevents.makeaxis(par)\ndt, dx, dy = par[\"dt\"], par[\"dx\"], par[\"dy\"]\n\n# Generate data\n_, d = pylops.utils.seismicevents.linear3d(x, y, t, 1500.0, t0, theta, phi, amp, wav)\n\nnpy, pymax = par[\"ny\"], 3e-4\nnpx, pxmax = par[\"nx\"], 5e-4\npy = np.linspace(-pymax, pymax, npy)\npx = np.linspace(-pxmax, pxmax, npx)\n\nR3Op = pylops.signalprocessing.ChirpRadon3D(\n    t, y, x, (pymax * dy / dt, pxmax * dx / dt), dtype=\"float64\"\n)\ndL_chirp = R3Op * d\ndadj_chirp = R3Op.H * dL_chirp\ndinv_chirp = R3Op.inverse(dL_chirp).reshape(R3Op.dimsd)\n\nfig, axs = plt.subplots(1, 4, figsize=(12, 4), sharey=True)\naxs[0].imshow(\n    d[par[\"ny\"] // 2].T,\n    vmin=-1,\n    vmax=1,\n    cmap=\"bwr_r\",\n    extent=(x[0], x[-1], t[-1], t[0]),\n)\naxs[0].set(xlabel=r\"$x$ [m]\", ylabel=r\"$t$ [s]\", title=\"Input model\")\naxs[0].axis(\"tight\")\naxs[1].imshow(\n    dL_chirp[par[\"ny\"] // 2].T,\n    cmap=\"bwr_r\",\n    vmin=-dL_chirp.max(),\n    vmax=dL_chirp.max(),\n    extent=(1e3 * px[0], 1e3 * px[-1], t[-1], t[0]),\n)\naxs[1].set(xlabel=r\"$p_x$ [s/km]\", title=\"Radon Chirp\")\naxs[1].axis(\"tight\")\naxs[2].imshow(\n    dadj_chirp[par[\"ny\"] // 2].T,\n    cmap=\"bwr_r\",\n    vmin=-dadj_chirp.max(),\n    vmax=dadj_chirp.max(),\n    extent=(x[0], x[-1], t[-1], t[0]),\n)\naxs[2].set(xlabel=r\"$x$ [m]\", title=\"Adj Radon Chirp\")\naxs[2].axis(\"tight\")\naxs[3].imshow(\n    dinv_chirp[par[\"ny\"] // 2].T,\n    cmap=\"bwr_r\",\n    vmin=-d.max(),\n    vmax=d.max(),\n    extent=(x[0], x[-1], t[-1], t[0]),\n)\naxs[3].set(xlabel=r\"$x$ [m]\", title=\"Inv Radon Chirp\")\naxs[3].axis(\"tight\")\nplt.tight_layout()\n\nfig, axs = plt.subplots(1, 4, figsize=(12, 4), sharey=True)\naxs[0].imshow(\n    d[:, par[\"nx\"] // 2].T,\n    vmin=-1,\n    vmax=1,\n    cmap=\"bwr_r\",\n    extent=(x[0], x[-1], t[-1], t[0]),\n)\naxs[0].set(xlabel=r\"$y$ [m]\", ylabel=r\"$t$ [s]\", title=\"Input model\")\naxs[0].axis(\"tight\")\naxs[1].imshow(\n    dL_chirp[:, 2 * par[\"nx\"] // 3].T,\n    cmap=\"bwr_r\",\n    vmin=-dL_chirp.max(),\n    vmax=dL_chirp.max(),\n    extent=(1e3 * py[0], 1e3 * py[-1], t[-1], t[0]),\n)\naxs[1].set(xlabel=r\"$p_y$ [s/km]\", title=\"Radon Chirp\")\naxs[1].axis(\"tight\")\naxs[2].imshow(\n    dadj_chirp[:, par[\"nx\"] // 2].T,\n    cmap=\"bwr_r\",\n    vmin=-dadj_chirp.max(),\n    vmax=dadj_chirp.max(),\n    extent=(x[0], x[-1], t[-1], t[0]),\n)\naxs[2].set(xlabel=r\"$y$ [m]\", title=\"Adj Radon Chirp\")\naxs[2].axis(\"tight\")\naxs[3].imshow(\n    dinv_chirp[:, par[\"nx\"] // 2].T,\n    cmap=\"bwr_r\",\n    vmin=-d.max(),\n    vmax=d.max(),\n    extent=(x[0], x[-1], t[-1], t[0]),\n)\naxs[3].set(xlabel=r\"$y$ [m]\", title=\"Inv Radon Chirp\")\naxs[3].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
}