feat: added parabolic3d method

docs/source/api/others.rst

@@ -118,6 +118,7 @@ Synthetics
     seismicevents.parabolic2d
     seismicevents.hyperbolic2d
     seismicevents.linear3d
+    seismicevents.parabolic3d
     seismicevents.hyperbolic3d
 
 .. currentmodule:: pylops.waveeqprocessing

examples/plot_seismicevents.py

@@ -106,11 +106,21 @@
 ############################################
 # Let's finally repeat the same exercise in 3d
 phi = [20, 0, -10]
+py = [0, 0, 0]
+pyy = [3e-5, 1e-5, 5e-6]
 
 mlin, mlinwav = pylops.utils.seismicevents.linear3d(
     x, y, t, v, t0, theta, phi, amp, wav
 )
 
+mpar, mparwav = pylops.utils.seismicevents.parabolic3d(
+    x, y, t, t0, px, py, pxx, pyy, amp, wav
+)
+
+mhyp, mhypwav = pylops.utils.seismicevents.hyperbolic3d(
+    x, y, t, t0, vrms, vrms, amp, wav
+)
+
 fig, axs = plt.subplots(1, 2, figsize=(7, 5), sharey=True)
 fig.suptitle("Linear events in 3d", fontsize=12, fontweight="bold", y=0.95)
 axs[0].imshow(
@@ -135,9 +145,30 @@
 )
 axs[1].set_xlabel(r"$y(m)$")
 
-mhyp, mhypwav = pylops.utils.seismicevents.hyperbolic3d(
-    x, y, t, t0, vrms, vrms, amp, wav
+fig, axs = plt.subplots(1, 2, figsize=(7, 5), sharey=True)
+fig.suptitle("Parabolic events in 3d", fontsize=12, fontweight="bold", y=0.95)
+axs[0].imshow(
+    mparwav[par["ny"] // 2].T,
+    aspect="auto",
+    interpolation="nearest",
+    vmin=-2,
+    vmax=2,
+    cmap="gray",
+    extent=(x.min(), x.max(), t.max(), t.min()),
+)
+axs[0].set_xlabel(r"$x(m)$")
+axs[0].set_ylabel(r"$t(s)$")
+axs[1].imshow(
+    mparwav[:, par["nx"] // 2].T,
+    aspect="auto",
+    interpolation="nearest",
+    vmin=-2,
+    vmax=2,
+    cmap="gray",
+    extent=(y.min(), y.max(), t.max(), t.min()),
 )
+axs[1].set_xlabel(r"$y(m)$")
+
 
 fig, axs = plt.subplots(1, 2, figsize=(7, 5), sharey=True)
 fig.suptitle("Hyperbolic events in 3d", fontsize=12, fontweight="bold", y=0.95)

pylops/signalprocessing/fourierradon2d.py

@@ -64,6 +64,11 @@ class FourierRadon2D(LinearOperator):
         Operator contains a matrix that can be solved explicitly (``True``) or
         not (``False``)
 
+    Raises
+    ------
+    NotImplementedError
+        If ``engine`` is neither ``numpy``, ``numba``, nor ``cuda``.
+
     Notes
     -----
     The FourierRadon2D operator applies the Radon transform in the frequency domain.
@@ -118,7 +123,7 @@ def __init__(
     ) -> None:
         # engine
         if engine not in ["numpy", "numba", "cuda"]:
-            raise KeyError("engine must be numpy or numba or cuda")
+            raise NotImplementedError("engine must be numpy or numba or cuda")
         if engine == "numba" and jit_message is not None:
             engine = "numpy"
 

pylops/signalprocessing/fourierradon3d.py

@@ -37,7 +37,7 @@ class FourierRadon3D(LinearOperator):
     taxis : :obj:`np.ndarray`
         Time axis
     hxaxis : :obj:`np.ndarray`
-        Fast patial axis
+        Fast spatial axis
     hyaxis : :obj:`np.ndarray`
         Slow spatial axis
     pyaxis : :obj:`np.ndarray`
@@ -51,7 +51,8 @@ class FourierRadon3D(LinearOperator):
         the application of the Radon matrix (if ``None``, use entire axis)
     kind : :obj:`tuple`, optional
         Curves to be used for stacking/spreading along the y- and x- axes
-        (``linear``, ``parabolic``)
+        (``("linear", "linear")``, ``("linear", "parabolic")``,
+         ``("parabolic", "linear")``, or  ``("parabolic", "parabolic")``)
     engine : :obj:`str`, optional
         Engine used for computation (``numpy`` or ``numba`` or ``cuda``)
     num_threads_per_blocks : :obj:`tuple`, optional
@@ -69,6 +70,13 @@ class FourierRadon3D(LinearOperator):
         Operator contains a matrix that can be solved explicitly (``True``) or
         not (``False``)
 
+    Raises
+    ------
+    NotImplementedError
+        If ``engine`` is neither ``numpy``, ``numba``, nor ``cuda``.
+    ValueError
+        If ``kind`` is not a tuple of two elements.
+
     Notes
     -----
     The FourierRadon3D operator applies the Radon transform in the frequency domain.
@@ -117,18 +125,22 @@ def __init__(
         pyaxis: NDArray,
         nfft: int,
         flims: Optional[Tuple[int, int]] = None,
-        kind: Optional[str] = ("linear", "linear"),
+        kind: Optional[tuple] = ("linear", "linear"),
         engine: Optional[str] = "numpy",
         num_threads_per_blocks: Tuple[int, int] = (32, 32),
         dtype: Optional[DTypeLike] = "float64",
         name: Optional[str] = "C",
     ) -> None:
         # engine
         if engine not in ["numpy", "numba", "cuda"]:
-            raise KeyError("engine must be numpy or numba or cuda")
+            raise NotImplementedError("engine must be numpy or numba or cuda")
         if engine == "numba" and jit_message is not None:
             engine = "numpy"
 
+        # kind
+        if not isinstance(kind, (tuple, list)) and len(kind) != 2:
+            raise ValueError("kind must be a tuple of two elements")
+
         # dimensions and super
         dims = len(pyaxis), len(pxaxis), len(taxis)
         dimsd = len(hyaxis), len(hxaxis), len(taxis)

pylops/utils/seismicevents.py

@@ -4,6 +4,7 @@
     "parabolic2d",
     "hyperbolic2d",
     "linear3d",
+    "parabolic3d",
     "hyperbolic3d",
 ]
 
@@ -160,7 +161,7 @@ def parabolic2d(
     r"""Parabolic 2D events
 
     Create 2d parabolic events given intercept time,
-    slowness, curvature, and amplitude of each event
+    slowness, curvature, and amplitude of each event.
 
     Parameters
     ----------
@@ -330,7 +331,7 @@ def linear3d(
     v : :obj:`float`
         propagation velocity
     t0 : :obj:`tuple` or :obj:`float`
-        intercept time at :math:`x=0` of each linear event
+        intercept time at :math:`x=0` and :math:`y=0` of each linear event
     theta : :obj:`tuple` or :obj:`float`
         angle in x direction (in degrees) of each linear event
     phi : :obj:`tuple` or :obj:`float`
@@ -397,6 +398,104 @@ def linear3d(
     return d, dwav
 
 
+def parabolic3d(
+    x: npt.NDArray,
+    y: npt.NDArray,
+    t: npt.NDArray,
+    t0: Union[float, Tuple[float]],
+    px: Union[float, Tuple[float]],
+    py: Union[float, Tuple[float]],
+    pxx: Union[float, Tuple[float]],
+    pyy: Union[float, Tuple[float]],
+    amp: Union[float, Tuple[float]],
+    wav: npt.NDArray,
+) -> Tuple[npt.NDArray, npt.NDArray]:
+    r"""Parabolic 3D events
+
+    Create 3d parabolic events given intercept time,
+    slowness, curvature, and amplitude of each event.
+
+    Parameters
+    ----------
+    x : :obj:`numpy.ndarray`
+        space axis in x direction
+    y : :obj:`numpy.ndarray`
+        space axis in y direction
+    t : :obj:`numpy.ndarray`
+        time axis
+    t0 : :obj:`tuple` or :obj:`float`
+        intercept time at :math:`x=0` and :math:`y=0` of each parabolic event
+    px : :obj:`tuple` or :obj:`float`
+        slowness of each parabolic event in x direction
+    py : :obj:`tuple` or :obj:`float`
+        slowness of each parabolic event in y direction
+    pxx : :obj:`tuple` or :obj:`float`
+        curvature of each parabolic event
+    amp : :obj:`tuple` or :obj:`float`
+        amplitude of each linear event
+    wav : :obj:`numpy.ndarray`
+        wavelet to be applied to data
+
+    Returns
+    -------
+    d : :obj:`numpy.ndarray`
+        data without wavelet of size
+        :math:`[n_y \times n_x \times n_t]`
+    dwav : :obj:`numpy.ndarray`
+        data with wavelet of size
+        :math:`[n_y \times n_x \times n_t]`
+
+    Notes
+    -----
+    Each event is created using the following relation:
+
+    .. math::
+        t_i(x, y) = t_{0,i} + p_{x,i} x + p_{y,i} x + p_{xx,i} x^2 + p_{yy,i} y^2
+
+    """
+    if isinstance(t0, (float, int)):
+        t0 = (t0,)
+    if isinstance(px, (float, int)):
+        px = (px,)
+    if isinstance(py, (float, int)):
+        py = (py,)
+    if isinstance(pxx, (float, int)):
+        pxx = (pxx,)
+    if isinstance(pyy, (float, int)):
+        pyy = (pyy,)
+
+    # identify dimensions
+    dt = t[1] - t[0]
+    wcenter = int(len(wav) / 2)
+    nx = np.size(x)
+    ny = np.size(y)
+    nt = np.size(t) + wcenter
+    nevents = np.size(t0)
+
+    # create events
+    d = np.zeros((ny, nx, nt))
+    for ievent in range(nevents):
+        for iy in range(ny):
+            tevent = (
+                t0[ievent]
+                + px[ievent] * x
+                + py[ievent] * y[iy]
+                + pxx[ievent] * x**2
+                + pyy[ievent] * y[iy] ** 2
+            )
+            tevent = (tevent - t[0]) / dt
+            itevent = tevent.astype(int)
+            dtevent = tevent - itevent
+            for ix in range(nx):
+                if itevent[ix] < nt - 1 and itevent[ix] >= 0:
+                    d[iy, ix, itevent[ix]] += amp[ievent] * (1 - dtevent[ix])
+                    d[iy, ix, itevent[ix] + 1] += amp[ievent] * dtevent[ix]
+
+    # filter events with certain wavelet
+    d, dwav = _filterdata(d, nt, wav, wcenter)
+    return d, dwav
+
+
 def hyperbolic3d(
     x: npt.NDArray,
     y: npt.NDArray,