{
  "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\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport matplotlib.pyplot as plt\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 = {'ot': 0,    'dt': 0.004, 'nt': 51,\n       'ox': -250, 'dx': 10, 'nx': 51,\n       'oy': -250, 'dy': 10, 'ny': 51,\n       'f0': 40}\ntheta = [0, ]\nt0 = [0.1, ]\namp = [1., ]\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 = \\\n    pylops.utils.seismicevents.linear2d(x, t, 1500., 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.ravel()\ndadj_chirp = R2Op.H * dL_chirp\ndinv_chirp = R2Op.inverse(dL_chirp)\n\ndL_chirp = dL_chirp.reshape(par['nx'], par['nt'])\ndadj_chirp = dadj_chirp.reshape(par['nx'], par['nt'])\ndinv_chirp = dinv_chirp.reshape(par['nx'], par['nt'])\n\nfig, axs = plt.subplots(1, 4, figsize=(12, 4))\naxs[0].imshow(d.T, vmin=-1, vmax=1, cmap='seismic_r',\n              extent=(x[0], x[-1], t[-1], t[0]))\naxs[0].set_title('Input model')\naxs[0].axis('tight')\naxs[1].imshow(dL_chirp.T, cmap='seismic_r', vmin=-dL_chirp.max(),\n              vmax=dL_chirp.max(),\n              extent=(px[0], px[-1], t[-1], t[0]))\naxs[1].set_title('Radon Chirp')\naxs[1].axis('tight')\naxs[2].imshow(dadj_chirp.T, cmap='seismic_r', vmin=-dadj_chirp.max(),\n              vmax=dadj_chirp.max(),\n              extent=(px[0], px[-1], t[-1], t[0]))\naxs[2].set_title('Adj Radon Chirp')\naxs[2].axis('tight')\naxs[3].imshow(dinv_chirp.T, cmap='seismic_r', vmin=-d.max(), vmax=d.max(),\n              extent=(px[0], px[-1], t[-1], t[0]))\naxs[3].set_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 = {'ot': 0,    'dt': 0.004, 'nt': 51,\n       'ox': -400, 'dx': 10, 'nx': 81,\n       'oy': -600, 'dy': 10, 'ny': 61,\n       'f0': 20}\ntheta = [10, ]\nphi = [0, ]\nt0 = [0.1, ]\namp = [1., ]\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 = \\\n    pylops.utils.seismicevents.linear3d(x, y, t, 1500., t0,\n                                        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(t, y, x,\n                                            (pymax*dy/dt, pxmax*dx/dt),\n                                            dtype='float64')\ndL_chirp = R3Op * d.ravel()\ndadj_chirp = R3Op.H * dL_chirp\ndinv_chirp = R3Op.inverse(dL_chirp)\n\ndL_chirp = dL_chirp.reshape(par['ny'], par['nx'], par['nt'])\ndadj_chirp = dadj_chirp.reshape(par['ny'], par['nx'], par['nt'])\ndinv_chirp = dinv_chirp.reshape(par['ny'], par['nx'], par['nt'])\n\n\nfig, axs = plt.subplots(1, 4, figsize=(12, 4))\naxs[0].imshow(d[par['ny']//2].T, vmin=-1, vmax=1, cmap='seismic_r',\n              extent=(x[0], x[-1], t[-1], t[0]))\naxs[0].set_title('Input model')\naxs[0].axis('tight')\naxs[1].imshow(dL_chirp[par['ny']//2].T, cmap='seismic_r',\n              vmin=-dL_chirp.max(), vmax=dL_chirp.max(),\n              extent=(px[0], px[-1], t[-1], t[0]))\naxs[1].set_title('Radon Chirp')\naxs[1].axis('tight')\naxs[2].imshow(dadj_chirp[par['ny']//2].T, cmap='seismic_r',\n              vmin=-dadj_chirp.max(), vmax=dadj_chirp.max(),\n              extent=(px[0], px[-1], t[-1], t[0]))\naxs[2].set_title('Adj Radon Chirp')\naxs[2].axis('tight')\naxs[3].imshow(dinv_chirp[par['ny']//2].T, cmap='seismic_r',\n              vmin=-d.max(), vmax=d.max(),\n              extent=(px[0], px[-1], t[-1], t[0]))\naxs[3].set_title('Inv Radon Chirp')\naxs[3].axis('tight')\nplt.tight_layout()\n\nfig, axs = plt.subplots(1, 4, figsize=(12, 4))\naxs[0].imshow(d[:, par['nx']//2].T, vmin=-1, vmax=1, cmap='seismic_r',\n              extent=(x[0], x[-1], t[-1], t[0]))\naxs[0].set_title('Input model')\naxs[0].axis('tight')\naxs[1].imshow(dL_chirp[:, 2*par['nx']//3].T, cmap='seismic_r',\n              vmin=-dL_chirp.max(), vmax=dL_chirp.max(),\n              extent=(px[0], px[-1], t[-1], t[0]))\naxs[1].set_title('Radon Chirp')\naxs[1].axis('tight')\naxs[2].imshow(dadj_chirp[:, par['nx']//2].T, cmap='seismic_r',\n              vmin=-dadj_chirp.max(), vmax=dadj_chirp.max(),\n              extent=(px[0], px[-1], t[-1], t[0]))\naxs[2].set_title('Adj Radon Chirp')\naxs[2].axis('tight')\naxs[3].imshow(dinv_chirp[:, par['nx']//2].T, cmap='seismic_r',\n              vmin=-d.max(), vmax=d.max(),\n              extent=(px[0], px[-1], t[-1], t[0]))\naxs[3].set_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.6.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}