From 86be5de12b792587a1f631fac3e3a62f41ec36c0 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sat, 30 Sep 2023 22:35:30 +0300
Subject: [PATCH 01/16] feature: restyling of Kirchhoff operator

---
 pylops/waveeqprocessing/kirchhoff.py | 607 +++++++--------------------
 pytests/test_kirchhoff.py            |  18 +-
 2 files changed, 149 insertions(+), 476 deletions(-)

diff --git a/pylops/waveeqprocessing/kirchhoff.py b/pylops/waveeqprocessing/kirchhoff.py
index a743ec6d..3ea6b6f4 100644
--- a/pylops/waveeqprocessing/kirchhoff.py
+++ b/pylops/waveeqprocessing/kirchhoff.py
@@ -35,10 +35,12 @@
 
 
 class Kirchhoff(LinearOperator):
-    r"""Kirchhoff Demigration operator.
+    r"""Kirchhoff demigration operator.
 
     Kirchhoff-based demigration/migration operator. Uses a high-frequency
-    approximation of  Green's function propagators based on ``trav``.
+    approximation of the Green's function propagators based on traveltimes
+    and amplitudes that are either computed internally by solving the Eikonal equation,
+    or passed directly by the user (which can use any other propagation engine of choice).
 
     Parameters
     ----------
@@ -73,23 +75,21 @@ class Kirchhoff(LinearOperator):
     dynamic : :obj:`bool`, optional
         .. versionadded:: 2.0.0
 
-        Include dynamic weights in computations (``True``) or not (``False``).
-        Note that when ``mode=byot``, the user is required to provide such weights
-        in ``amp``.
+        Apply dynamic weights in computations (``True``) or not (``False``). This includes both the amplitude
+        terms of the Green's function and the reflectivity-related scaling term (see equations below).
     trav : :obj:`numpy.ndarray` or :obj:`tuple`, optional
         Traveltime table of size
         :math:`\lbrack (n_y) n_x n_z \times n_s n_r \rbrack` or pair of traveltime tables
         of size :math:`\lbrack (n_y) n_x n_z \times n_s \rbrack` and :math:`\lbrack (n_y) n_x n_z \times n_r \rbrack`
         (to be provided if ``mode='byot'``). Note that the latter approach is recommended as less memory demanding
-        than the former.
+        than the former. Moreover, only ``mode='dynamic'`` is only possible when traveltimes are provided in
+        the latter form.
     amp : :obj:`numpy.ndarray`, optional
         .. versionadded:: 2.0.0
 
-        Amplitude table of size
-        :math:`\lbrack (n_y) n_x n_z \times n_s n_r \rbrack` or pair of amplitude tables
-        of size :math:`\lbrack (n_y) n_x n_z \times n_s \rbrack` and :math:`\lbrack (n_y) n_x n_z \times n_r \rbrack`
-        (to be provided if ``mode='byot'``). Note that the latter approach is recommended as less memory demanding
-        than the former.
+        Pair of amplitude tables of size :math:`\lbrack (n_y) n_x n_z \times n_s \rbrack` and
+        :math:`\lbrack (n_y) n_x n_z \times n_r \rbrack` (to be provided if ``mode='byot'``). Note that this parameter
+        is only required when ``mode='dynamic'`` is chosen.
     aperture : :obj:`float` or :obj:`tuple`, optional
         .. versionadded:: 2.0.0
 
@@ -104,18 +104,9 @@ class Kirchhoff(LinearOperator):
         Maximum allowed angle (either source or receiver side) in degrees. If ``None``,
         angle aperture limitations are not introduced. See ``aperture`` for implementation
         details regarding scalar and tuple cases.
-
-    anglerefl : :obj:`np.ndarray`, optional
-        .. versionadded:: 2.0.0
-
-        Angle between the normal of the reflectors and the vertical axis in degrees
     snell : :obj:`float` or :obj:`tuple`, optional
-        .. versionadded:: 2.0.0
-
-        Threshold on Snell's law evaluation. If larger, the source-receiver-image
-        point is discarded. If ``None``, no check on the validity of the Snell's
-        law is performed.  See ``aperture`` for implementation details regarding
-        scalar and tuple cases.
+        Deprecated, will be removed in v3.0.0. Simply kept for back-compatibility with previous implementation,
+        but effectively not affecting the behaviour of the operator.
     engine : :obj:`str`, optional
         Engine used for computations (``numpy`` or ``numba``).
     dtype : :obj:`str`, optional
@@ -141,29 +132,34 @@ class Kirchhoff(LinearOperator):
     Notes
     -----
     The Kirchhoff demigration operator synthesizes seismic data given a
-    propagation velocity model :math:`v` and a reflectivity model :math:`m`.
-    In forward mode [1]_, [2]_:
+    propagation velocity model :math:`v(\mathbf{x})` and a
+    reflectivity model :math:`m(\mathbf{x})`. In forward mode [1]_, [2]_:
 
     .. math::
         d(\mathbf{x_r}, \mathbf{x_s}, t) =
-        \widetilde{w}(t) * \int_V G(\mathbf{x_r}, \mathbf{x}, t)
-        m(\mathbf{x}) G(\mathbf{x}, \mathbf{x_s}, t)\,\mathrm{d}\mathbf{x}
+        \widetilde{w}(t) * \int_V \frac{2 cos(\theta)} {v(\mathbf{x})}
+        G(\mathbf{x_r}, \mathbf{x}, t) G(\mathbf{x}, \mathbf{x_s}, t)
+        m(\mathbf{x})  \,\mathrm{d}\mathbf{x}
+
+    where :math:`G(\mathbf{x}, \mathbf{x_s}, t)` and :math:`G(\mathbf{x_r}, \mathbf{x}, t)`
+    are the Green's functions from source-to-subsurface-to-receiver and finally
+    :math:`\widetilde{w}(t)` is a filtered version of the wavelet :math:`w(t)`
+    as explained below (or the wavelet itself when ``wavfilter=False``).
+    Moreover, an angle scaling is included in the modelling operator,
+    where :math:`\theta=(\theta_s-\theta_r)/2` is half of the opening angle,
+    where :math:`\theta_s` and :math:`\theta_r` are the angles between the source-side
+    and receiver-side rays and the vertical at the image point.
 
-    where :math:`m(\mathbf{x})` represents the reflectivity
-    at every location in the subsurface, :math:`G(\mathbf{x}, \mathbf{x_s}, t)`
-    and :math:`G(\mathbf{x_r}, \mathbf{x}, t)` are the Green's functions
-    from source-to-subsurface-to-receiver and finally :math:`\widetilde{w}(t)` is
-    a filtered version of the wavelet :math:`w(t)` [3]_ (or the wavelet itself when
-    ``wavfilter=False``). In our implementation, the following high-frequency
-    approximation of the Green's functions is adopted:
+    In our implementation, the following high-frequency approximation of the
+    Green's functions is adopted:
 
     .. math::
         G(\mathbf{x_r}, \mathbf{x}, \omega) = a(\mathbf{x_r}, \mathbf{x})
             e^{j \omega t(\mathbf{x_r}, \mathbf{x})}
 
     where :math:`a(\mathbf{x_r}, \mathbf{x})` is the amplitude and
-    :math:`t(\mathbf{x_r}, \mathbf{x})` is the traveltime. When ``dynamic=False`` the
-    amplitude is disregarded leading to a kinematic-only Kirchhoff operator.
+    :math:`t(\mathbf{x_r}, \mathbf{x})` is the traveltime. When ``dynamic=False``, the
+    amplitude correction terms are disregarded leading to a kinematic-only Kirchhoff operator.
 
     .. math::
         d(\mathbf{x_r}, \mathbf{x_s}, t) =
@@ -171,36 +167,32 @@ class Kirchhoff(LinearOperator):
         t(\mathbf{x}, \mathbf{x_s}))} m(\mathbf{x}) \,\mathrm{d}\mathbf{x}
 
     On the  other hand, when ``dynamic=True``, the amplitude scaling is defined as
-    :math:`a(\mathbf{x}, \mathbf{y})=\frac{1}{\|\mathbf{x} - \mathbf{y}\|}`,
-    that is, the reciprocal of the distance between the two points,
-    approximating the geometrical spreading of the wavefront.
-    Moreover an angle scaling is included in the modelling operator
-    added as follows:
 
-    .. math::
-        d(\mathbf{x_r}, \mathbf{x_s}, t) =
-        \tilde{w}(t) * \int_V a(\mathbf{x}, \mathbf{x_s}) a(\mathbf{x}, \mathbf{x_r})
-        \frac{|cos \theta_s + cos \theta_r|} {v(\mathbf{x})} e^{j \omega (t(\mathbf{x_r}, \mathbf{x}) +
-         t(\mathbf{x}, \mathbf{x_s}))} m(\mathbf{x}) \,\mathrm{d}\mathbf{x}
+    * ``2D``: :math:`a(\mathbf{x}, \mathbf{y})=\frac{1}{\sqrt{dist(\mathbf{x}, \mathbf{y})}}`
+    * ``3D``: :math:`a(\mathbf{x}, \mathbf{y})=\frac{1}{dist(\mathbf{x}, \mathbf{y})}`
 
-    where :math:`\theta_s` and :math:`\theta_r` are the angles between the source-side
-    and receiver-side rays and the normal to the reflector  at the image point (or
-    the vertical axis at the image point when ``reflslope=None``), respectively.
+    approximating the geometrical spreading of the wavefront. For ``mode=analytic``,
+    :math:`dist(\mathbf{x}, \mathbf{y})=\|\mathbf{x} - \mathbf{y}\|`, whilst for
+    ``mode=eikonal``, this is instead computed internally by the eikonal solver.
+
+    The wavelet filtering is applied as follows:
+
+    * ``2D``: :math:`\tilde{W}(f)=\sqrt{j\omega} \cdot W(f)`
+    * ``3D``: :math:`\tilde{W}(f)=j\omega \cdot W(f)`
 
     Depending on the choice of ``mode`` the traveltime and amplitude of the Green's
     function will be also computed differently:
 
-    * ``mode=analytic`` or ``mode=eikonal``: traveltimes, geometrical spreading, and angles
-      are computed for every source-image point-receiver triplets and the
-      Green's functions are implemented from traveltime look-up tables, placing
-      scaled reflectivity values at corresponding source-to-receiver time in the
-      data.
-    * ``byot``: bring your own tables. Traveltime table are provided
+    * ``mode=analytic`` or ``mode=eikonal``: traveltimes, amplitudes, and angles
+      are computed for every source-image point-receiver triplets upfront and the
+      Green's functions are implemented from traveltime and amplitude look-up tables,
+      placing scaled reflectivity values at corresponding source-to-receiver time
+      in the data.
+    * ``mode=byot``: bring your own tables. Traveltime table are provided
       directly by user using ``trav`` input parameter. Similarly, in this case one
-      can provide their own amplitude scaling ``amp`` (which should include the angle
-      scaling too).
+      can also provide their own amplitude scaling ``amp``.
 
-    Three aperture limitations have been also implemented as defined by:
+    Two aperture limitations have been also implemented as defined by:
 
     * ``aperture``: the maximum allowed aperture is expressed as the ratio of
       offset over depth. This aperture limitation avoid including grazing angles
@@ -208,15 +200,10 @@ class Kirchhoff(LinearOperator):
       edges of the aperture;
     * ``angleaperture``: the maximum allowed angle aperture is expressed as the
       difference between the incident or emerging angle at every image point and
-      the vertical axis (or the normal to the reflector if ``anglerefl`` is provided.
-      This aperture limitation also avoid including grazing angles whose contributions
-      can introduce aliasing effects. Note that for a homogenous medium and slowly varying
-      heterogenous medium the offset and angle aperture limits may work in the same way;
-    * ``snell``: the maximum allowed snell's angle is expressed as the absolute value of
-      the sum between incident and emerging angles defined as in the ``angleaperture`` case.
-      This aperture limitation is introduced to turn a scattering-based Kirchhoff engine into
-      a reflection-based Kirchhoff engine where each image point is not considered as scatter
-      but as a local horizontal (or straight) reflector.
+      the vertical axis. This aperture limitation also avoid including grazing angles
+      whose contributions can introduce aliasing effects. Note that for a homogenous
+      medium and slowly varying heterogeneous medium the offset and angle aperture limits
+      may work in the same way.
 
     Finally, the adjoint of the demigration operator is a *migration* operator which
     projects data in the model domain creating an image of the subsurface
@@ -252,7 +239,6 @@ def __init__(
         amp: Optional[NDArray] = None,
         aperture: Optional[Tuple[float, float]] = None,
         angleaperture: Union[float, Tuple[float, float]] = 90.0,
-        anglerefl: Optional[NDArray] = None,
         snell: Optional[Tuple[float, float]] = None,
         engine: str = "numpy",
         dtype: DTypeLike = "float64",
@@ -288,17 +274,20 @@ def __init__(
         self.dt = t[1] - t[0]
         self.nt = len(t)
 
-        # store ix-iy locations of sources and receivers
-        dx = x[1] - x[0]
-        if self.ndims == 2:
-            self.six = np.tile((srcs[0] - x[0]) // dx, (nr, 1)).T.astype(int).ravel()
-            self.rix = np.tile((recs[0] - x[0]) // dx, (ns, 1)).astype(int).ravel()
-        elif self.ndims == 3:
-            # TODO: 3D normalized distances
-            pass
-
-        # compute traveltime
+        # store ix-iy locations of sources and receivers for aperture filter
         self.dynamic = dynamic
+        if self.dynamic:
+            dx = x[1] - x[0]
+            if self.ndims == 2:
+                self.six = (
+                    np.tile((srcs[0] - x[0]) // dx, (nr, 1)).T.astype(int).ravel()
+                )
+                self.rix = np.tile((recs[0] - x[0]) // dx, (ns, 1)).astype(int).ravel()
+            elif self.ndims == 3:
+                # TODO: 3D normalized distances
+                raise NotImplementedError("dynamic=True currently not available in 3D")
+
+        # compute traveltime and distances
         self.travsrcrec = True  # use separate tables for src and rec traveltimes
         if mode in ["analytic", "eikonal", "byot"]:
             if mode in ["analytic", "eikonal"]:
@@ -306,63 +295,72 @@ def __init__(
                 (
                     self.trav_srcs,
                     self.trav_recs,
-                    self.dist_srcs,
-                    self.dist_recs,
+                    dist_srcs,
+                    dist_recs,
                     trav_srcs_grad,
                     trav_recs_grad,
                 ) = Kirchhoff._traveltime_table(z, x, srcs, recs, vel, y=y, mode=mode)
                 if self.dynamic:
                     # need to add a scalar in the denominator of amplitude term to avoid
-                    # division by 0, currently set to 1/100 of max distance to avoid having
-                    # to large scaling around the source. This number may change in future
+                    # division by 0, currently set to 1e-4 of max distance to avoid having
+                    # too large scaling around the source. This number may change in future
                     # or left to the user to define
                     epsdist = 1e-2
-                    self.maxdist = epsdist * (
-                        np.max(self.dist_srcs) + np.max(self.dist_recs)
-                    )
-                    # compute angles
+                    self.maxdist = epsdist * (np.max(dist_srcs) + np.max(dist_recs))
                     if self.ndims == 2:
-                        # 2d with vertical
-                        if anglerefl is None:
-                            self.angle_srcs = np.arctan2(
-                                trav_srcs_grad[0], trav_srcs_grad[1]
-                            ).reshape(np.prod(dims), ns)
-                            self.angle_recs = np.arctan2(
-                                trav_recs_grad[0], trav_recs_grad[1]
-                            ).reshape(np.prod(dims), nr)
-                            self.cosangle_srcs = np.cos(self.angle_srcs)
-                            self.cosangle_recs = np.cos(self.angle_recs)
-                        else:
-                            # TODO: 2D with normal
-                            raise NotImplementedError(
-                                "angle scaling with anglerefl currently not available"
-                            )
+                        self.amp_srcs, self.amp_recs = 1.0 / np.sqrt(
+                            dist_srcs + self.maxdist
+                        ), 1.0 / np.sqrt(dist_recs + self.maxdist)
                     else:
-                        # TODO: 3D
-                        raise NotImplementedError(
-                            "dynamic=True currently not available in 3D"
-                        )
+                        self.amp_srcs, self.amp_recs = 1.0 / (
+                            dist_srcs + self.maxdist
+                        ), 1.0 / (dist_recs + self.maxdist)
             else:
                 if isinstance(trav, tuple):
                     self.trav_srcs, self.trav_recs = trav
                 else:
                     self.travsrcrec = False
                     self.trav = trav
+
+                if self.dynamic and not self.travsrcrec:
+                    raise NotImplementedError(
+                        "separate traveltime tables must be provided "
+                        "when selecting mode=dynamic"
+                    )
                 if self.dynamic:
                     if isinstance(amp, tuple):
                         self.amp_srcs, self.amp_recs = amp
                     else:
-                        self.amp = amp
-                    # in byot mode, angleaperture and snell checks are not performed
-                    self.angle_srcs = np.ones(
-                        (self.ny * self.nx * self.nz, ns), dtype=dtype
-                    )
-                    self.angle_recs = np.ones(
-                        (self.ny * self.nx * self.nz, nr), dtype=dtype
-                    )
+                        raise NotImplementedError(
+                            "separate amplitude tables must be provided "
+                        )
+
+                    if self.travsrcrec:
+                        # compute traveltime gradients at image points
+                        trav_srcs_grad = np.gradient(
+                            self.trav_srcs.reshape(*dims, ns),
+                            axis=np.arange(self.ndims),
+                        )
+                        trav_recs_grad = np.gradient(
+                            self.trav_recs.reshape(*dims, nr),
+                            axis=np.arange(self.ndims),
+                        )
         else:
             raise NotImplementedError("method must be analytic, eikonal or byot")
 
+        # compute angles with vertical
+        if self.dynamic:
+            if self.ndims == 2:
+                self.angle_srcs = np.arctan2(
+                    trav_srcs_grad[0], trav_srcs_grad[1]
+                ).reshape(np.prod(dims), ns)
+                self.angle_recs = np.arctan2(
+                    trav_recs_grad[0], trav_recs_grad[1]
+                ).reshape(np.prod(dims), nr)
+            else:
+                # TODO: 3D
+                raise NotImplementedError("dynamic=True currently not available in 3D")
+
         # pre-compute traveltime indices if total traveltime is used
         if not self.travsrcrec:
             self.itrav = (self.trav / self.dt).astype("int32")
@@ -383,31 +381,25 @@ def __init__(
         self.aperturetap = taper(41, 20, "hanning")[20:]
 
         # define aperture
+        # if aperture=None, we want to ensure the check is always matched (no aperture limits...)
         if aperture is not None:
             warnings.warn(
                 "Aperture is currently defined as ratio of offset over depth, "
-                "and may be not ideal for highly heterogenous media"
+                "and may be not ideal for highly heterogeneous media"
             )
         self.aperture = (
-            (2 * self.nx / self.nz,)
-            if aperture is None
-            else _value_or_sized_to_array(aperture)
+            (self.nx + 1,) if aperture is None else _value_or_sized_to_array(aperture)
         )
         if len(self.aperture) == 1:
             self.aperture = np.array([0.8 * self.aperture[0], self.aperture[0]])
 
-        # define angle aperture and snell law
+        # define angle aperture
         angleaperture = [0.0, 1000.0] if angleaperture is None else angleaperture
         self.angleaperture = np.deg2rad(_value_or_sized_to_array(angleaperture))
         if len(self.angleaperture) == 1:
             self.angleaperture = np.array(
                 [0.8 * self.angleaperture[0], self.angleaperture[0]]
             )
-        self.snell = (
-            (np.pi,) if snell is None else np.deg2rad(_value_or_sized_to_array(snell))
-        )
-        if len(self.snell) == 1:
-            self.snell = np.array([0.8 * self.snell[0], self.snell[0]])
 
         # dimensions
         self.ns, self.nr = ns, nr
@@ -635,18 +627,19 @@ def _wavelet_reshaping(
         dimrec: int,
         dimv: int,
     ) -> NDArray:
-        """Apply wavelet reshaping as from theory in [1]_"""
+        """Apply wavelet reshaping to account for omega scaling factor
+        originating from the wave equation"""
         f = np.fft.rfftfreq(len(wav), dt)
         W = np.fft.rfft(wav, n=len(wav))
         if dimsrc == 2 and dimv == 2:
             # 2D
-            Wfilt = W * (2 * np.pi * f)
+            Wfilt = W * np.sqrt(1j * 2 * np.pi * f)
         elif (dimsrc == 2 or dimrec == 2) and dimv == 3:
             # 2.5D
             raise NotImplementedError("2.D wavelet currently not available")
         elif dimsrc == 3 and dimrec == 3 and dimv == 3:
             # 3D
-            Wfilt = W * (-1j * 2 * np.pi * f)
+            Wfilt = W * (1j * 2 * np.pi * f)
         wavfilt = np.fft.irfft(Wfilt, n=len(wav))
         return wavfilt
 
@@ -750,225 +743,6 @@ def _travsrcrec_kirch_rmatvec(
                         )
         return y
 
-    @staticmethod
-    def _amp_kirch_matvec(
-        x: NDArray,
-        y: NDArray,
-        nsnr: int,
-        nt: int,
-        ni: int,
-        itrav: NDArray,
-        travd: NDArray,
-        amp: NDArray,
-        aperturemin: float,
-        aperturemax: float,
-        aperturetap: NDArray,
-        nz: int,
-        six: NDArray,
-        rix: NDArray,
-        angleaperturemin: float,
-        angleaperturemax: float,
-        angles_srcs: NDArray,
-        angles_recs: NDArray,
-        snellmin: float,
-        snellmax: float,
-    ) -> NDArray:
-        nr = angles_recs.shape[-1]
-        daperture = aperturemax - aperturemin
-        dangleaperture = angleaperturemax - angleaperturemin
-        dsnell = snellmax - snellmin
-        for isrcrec in prange(nsnr):
-            # extract traveltime, amplitude, src/rec coordinates at given src/pair
-            itravisrcrec = itrav[:, isrcrec]
-            travdisrcrec = travd[:, isrcrec]
-            ampisrcrec = amp[:, isrcrec]
-            sixisrcrec = six[isrcrec]
-            rixisrcrec = rix[isrcrec]
-            # extract source and receiver angles
-            angles_src = angles_srcs[:, isrcrec // nr]
-            angles_rec = angles_recs[:, isrcrec % nr]
-            for ii in range(ni):
-                # extract traveltime, amplitude at given image point
-                itravii = itravisrcrec[ii]
-                travdii = travdisrcrec[ii]
-                damp = ampisrcrec[ii]
-                # extract source and receiver angle
-                angle_src = angles_src[ii]
-                angle_rec = angles_rec[ii]
-                abs_angle_src = abs(angle_src)
-                abs_angle_rec = abs(angle_rec)
-                abs_angle_src_rec = abs(angle_src + angle_rec)
-                aptap = 1.0
-                # angle apertures checks
-                if (
-                    abs_angle_src < angleaperturemax
-                    and abs_angle_rec < angleaperturemax
-                    and abs_angle_src_rec < snellmax
-                ):
-                    if abs_angle_src >= angleaperturemin:
-                        # extract source angle aperture taper value
-                        aptap = (
-                            aptap
-                            * aperturetap[
-                                int(
-                                    20
-                                    * (abs_angle_src - angleaperturemin)
-                                    // dangleaperture
-                                )
-                            ]
-                        )
-                    if abs_angle_rec >= angleaperturemin:
-                        # extract receiver angle aperture taper value
-                        aptap = (
-                            aptap
-                            * aperturetap[
-                                int(
-                                    20
-                                    * (abs_angle_rec - angleaperturemin)
-                                    // dangleaperture
-                                )
-                            ]
-                        )
-                    if abs_angle_src_rec >= snellmin:
-                        # extract snell taper value
-                        aptap = (
-                            aptap
-                            * aperturetap[
-                                int(20 * (abs_angle_src_rec - snellmin) // dsnell)
-                            ]
-                        )
-
-                    # identify x-index of image point
-                    iz = ii % nz
-                    # aperture check
-                    aperture = abs(sixisrcrec - rixisrcrec) / iz
-                    if aperture < aperturemax:
-                        if aperture >= aperturemin:
-                            # extract aperture taper value
-                            aptap = (
-                                aptap
-                                * aperturetap[
-                                    int(20 * ((aperture - aperturemin) // daperture))
-                                ]
-                            )
-                        # time limit check
-                        if 0 <= itravii < nt - 1:
-                            # assign values
-                            y[isrcrec, itravii] += x[ii] * (1 - travdii) * damp * aptap
-                            y[isrcrec, itravii + 1] += x[ii] * travdii * damp * aptap
-        return y
-
-    @staticmethod
-    def _amp_kirch_rmatvec(
-        x: NDArray,
-        y: NDArray,
-        nsnr: int,
-        nt: int,
-        ni: int,
-        itrav: NDArray,
-        travd: NDArray,
-        amp: NDArray,
-        aperturemin: float,
-        aperturemax: float,
-        aperturetap: NDArray,
-        nz: int,
-        six: NDArray,
-        rix: NDArray,
-        angleaperturemin: float,
-        angleaperturemax: float,
-        angles_srcs: NDArray,
-        angles_recs: NDArray,
-        snellmin: float,
-        snellmax: float,
-    ) -> NDArray:
-        nr = angles_recs.shape[-1]
-        daperture = aperturemax - aperturemin
-        dangleaperture = angleaperturemax - angleaperturemin
-        dsnell = snellmax - snellmin
-        for ii in prange(ni):
-            itravii = itrav[ii]
-            travdii = travd[ii]
-            ampii = amp[ii]
-            # extract source and receiver angles
-            angle_srcs = angles_srcs[ii]
-            angle_recs = angles_recs[ii]
-            # identify x-index of image point
-            iz = ii % nz
-            for isrcrec in range(nsnr):
-                itravisrcrecii = itravii[isrcrec]
-                travdisrcrecii = travdii[isrcrec]
-                sixisrcrec = six[isrcrec]
-                rixisrcrec = rix[isrcrec]
-                # extract source and receiver angle
-                angle_src = angle_srcs[isrcrec // nr]
-                angle_rec = angle_recs[isrcrec % nr]
-                abs_angle_src = abs(angle_src)
-                abs_angle_rec = abs(angle_rec)
-                abs_angle_src_rec = abs(angle_src + angle_rec)
-                aptap = 1.0
-                # angle apertures checks
-                if (
-                    abs_angle_src < angleaperturemax
-                    and abs_angle_rec < angleaperturemax
-                    and abs_angle_src_rec < snellmax
-                ):
-                    if abs_angle_src >= angleaperturemin:
-                        # extract source angle aperture taper value
-                        aptap = (
-                            aptap
-                            * aperturetap[
-                                int(
-                                    20
-                                    * (abs_angle_src - angleaperturemin)
-                                    // dangleaperture
-                                )
-                            ]
-                        )
-                    if abs_angle_rec >= angleaperturemin:
-                        # extract receiver angle aperture taper value
-                        aptap = (
-                            aptap
-                            * aperturetap[
-                                int(
-                                    20
-                                    * (abs_angle_rec - angleaperturemin)
-                                    // dangleaperture
-                                )
-                            ]
-                        )
-                    if abs_angle_src_rec >= snellmin:
-                        # extract snell taper value
-                        aptap = (
-                            aptap
-                            * aperturetap[
-                                int(20 * (abs_angle_src_rec - snellmin) // dsnell)
-                            ]
-                        )
-
-                    # aperture check
-                    aperture = abs(sixisrcrec - rixisrcrec) / iz
-                    if aperture < aperturemax:
-                        if aperture >= aperturemin:
-                            # extract aperture taper value
-                            aptap = (
-                                aptap
-                                * aperturetap[
-                                    int(20 * ((aperture - aperturemin) // daperture))
-                                ]
-                            )
-                        # time limit check
-                        if 0 <= itravisrcrecii < nt - 1:
-                            # assign values
-                            y[ii] += (
-                                (
-                                    x[isrcrec, itravisrcrecii] * (1 - travdisrcrecii)
-                                    + x[isrcrec, itravisrcrecii + 1] * travdisrcrecii
-                                )
-                                * ampii[isrcrec]
-                                * aptap
-                            )
-        return y
-
     @staticmethod
     def _ampsrcrec_kirch_matvec(
         x: NDArray,
@@ -981,8 +755,8 @@ def _ampsrcrec_kirch_matvec(
         vel: NDArray,
         trav_srcs: NDArray,
         trav_recs: NDArray,
-        dist_srcs: NDArray,
-        dist_recs: NDArray,
+        amp_srcs: NDArray,
+        amp_recs: NDArray,
         aperturemin: float,
         aperturemax: float,
         aperturetap: NDArray,
@@ -993,44 +767,39 @@ def _ampsrcrec_kirch_matvec(
         angleaperturemax: float,
         angles_srcs: NDArray,
         angles_recs: NDArray,
-        snellmin: float,
-        snellmax: float,
-        maxdist: float,
     ) -> NDArray:
         daperture = aperturemax - aperturemin
         dangleaperture = angleaperturemax - angleaperturemin
-        dsnell = snellmax - snellmin
         for isrc in prange(ns):
             travisrc = trav_srcs[:, isrc]
-            distisrc = dist_srcs[:, isrc]
+            ampisrc = amp_srcs[:, isrc]
             angleisrc = angles_srcs[:, isrc]
             for irec in range(nr):
                 travirec = trav_recs[:, irec]
                 trav = travisrc + travirec
                 itrav = (trav / dt).astype("int32")
                 travd = trav / dt - itrav
-                distirec = dist_recs[:, irec]
+                ampirec = amp_recs[:, irec]
                 angleirec = angles_recs[:, irec]
-                dist = distisrc + distirec
-                amp = np.abs(np.cos(angleisrc) + np.cos(angleirec)) / (dist + maxdist)
                 sixisrcrec = six[isrc * nr + irec]
                 rixisrcrec = rix[isrc * nr + irec]
+                # compute cosine of half opening angle and total amplitude scaling
+                cosangle = np.cos((angleisrc - angleirec) / 2.0)
+                amp = 2.0 * cosangle * ampisrc * ampirec / vel
                 for ii in range(ni):
                     itravii = itrav[ii]
                     travdii = travd[ii]
-                    damp = amp[ii] / vel[ii]
+                    damp = amp[ii]
                     # extract source and receiver angle at given image point
                     angle_src = angleisrc[ii]
                     angle_rec = angleirec[ii]
                     abs_angle_src = abs(angle_src)
                     abs_angle_rec = abs(angle_rec)
-                    abs_angle_src_rec = abs(angle_src + angle_rec)
-                    aptap = 1.0
                     # angle apertures checks
+                    aptap = 1.0
                     if (
                         abs_angle_src < angleaperturemax
                         and abs_angle_rec < angleaperturemax
-                        and abs_angle_src_rec < snellmax
                     ):
                         if abs_angle_src >= angleaperturemin:
                             # extract source angle aperture taper value
@@ -1056,19 +825,11 @@ def _ampsrcrec_kirch_matvec(
                                     )
                                 ]
                             )
-                        if abs_angle_src_rec >= snellmin:
-                            # extract snell taper value
-                            aptap = (
-                                aptap
-                                * aperturetap[
-                                    int(20 * (abs_angle_src_rec - snellmin) // dsnell)
-                                ]
-                            )
 
                         # identify x-index of image point
                         iz = ii % nz
                         # aperture check
-                        aperture = abs(sixisrcrec - rixisrcrec) / iz
+                        aperture = abs(sixisrcrec - rixisrcrec) / (iz + 1)
                         if aperture < aperturemax:
                             if aperture >= aperturemin:
                                 # extract aperture taper value
@@ -1102,8 +863,8 @@ def _ampsrcrec_kirch_rmatvec(
         vel: NDArray,
         trav_srcs: NDArray,
         trav_recs: NDArray,
-        dist_srcs: NDArray,
-        dist_recs: NDArray,
+        amp_srcs: NDArray,
+        amp_recs: NDArray,
         aperturemin: float,
         aperturemax: float,
         aperturetap: NDArray,
@@ -1114,18 +875,14 @@ def _ampsrcrec_kirch_rmatvec(
         angleaperturemax: float,
         angles_srcs: NDArray,
         angles_recs: NDArray,
-        snellmin: float,
-        snellmax: float,
-        maxdist: float,
     ) -> NDArray:
         daperture = aperturemax - aperturemin
         dangleaperture = angleaperturemax - angleaperturemin
-        dsnell = snellmax - snellmin
         for ii in prange(ni):
             trav_srcsii = trav_srcs[ii]
             trav_recsii = trav_recs[ii]
-            dist_srcsii = dist_srcs[ii]
-            dist_recsii = dist_recs[ii]
+            amp_srcsii = amp_srcs[ii]
+            amp_recsii = amp_recs[ii]
             velii = vel[ii]
             angle_srcsii = angles_srcs[ii]
             angle_recsii = angles_recs[ii]
@@ -1133,31 +890,27 @@ def _ampsrcrec_kirch_rmatvec(
             iz = ii % nz
             for isrc in range(ns):
                 trav_srcii = trav_srcsii[isrc]
-                dist_srcii = dist_srcsii[isrc]
+                amp_srcii = amp_srcsii[isrc]
                 angle_src = angle_srcsii[isrc]
                 for irec in range(nr):
                     trav_recii = trav_recsii[irec]
                     travii = trav_srcii + trav_recii
                     itravii = int(travii / dt)
                     travdii = travii / dt - itravii
-                    dist_recii = dist_recsii[irec]
+                    amp_recii = amp_recsii[irec]
                     angle_rec = angle_recsii[irec]
-                    dist = dist_srcii + dist_recii
-                    ampii = np.abs(np.cos(angle_src) + np.cos(angle_rec)) / (
-                        dist + maxdist
-                    )
-                    damp = ampii / velii
                     sixisrcrec = six[isrc * nr + irec]
                     rixisrcrec = rix[isrc * nr + irec]
                     abs_angle_src = abs(angle_src)
                     abs_angle_rec = abs(angle_rec)
-                    abs_angle_src_rec = abs(angle_src + angle_rec)
-                    aptap = 1.0
+                    # compute cosine of half opening angle and total amplitude scaling
+                    cosangle = np.cos((angle_src - angle_rec) / 2.0)
+                    damp = 2.0 * cosangle * amp_srcii * amp_recii / velii
                     # angle apertures checks
+                    aptap = 1.0
                     if (
                         abs_angle_src < angleaperturemax
                         and abs_angle_rec < angleaperturemax
-                        and abs_angle_src_rec < snellmax
                     ):
                         if abs_angle_src >= angleaperturemin:
                             # extract source angle aperture taper value
@@ -1183,17 +936,9 @@ def _ampsrcrec_kirch_rmatvec(
                                     )
                                 ]
                             )
-                        if abs_angle_src_rec >= snellmin:
-                            # extract snell taper value
-                            aptap = (
-                                aptap
-                                * aperturetap[
-                                    int(20 * (abs_angle_src_rec - snellmin) // dsnell)
-                                ]
-                            )
 
                         # aperture check
-                        aperture = abs(sixisrcrec - rixisrcrec) / iz
+                        aperture = abs(sixisrcrec - rixisrcrec) / (iz + 1)
                         if aperture < aperturemax:
                             if aperture >= aperturemin:
                                 # extract aperture taper value
@@ -1229,9 +974,6 @@ def _register_multiplications(self, engine: str) -> None:
             if self.dynamic and self.travsrcrec:
                 self._kirch_matvec = jit(**numba_opts)(self._ampsrcrec_kirch_matvec)
                 self._kirch_rmatvec = jit(**numba_opts)(self._ampsrcrec_kirch_rmatvec)
-            elif self.dynamic and not self.travsrcrec:
-                self._kirch_matvec = jit(**numba_opts)(self._amp_kirch_matvec)
-                self._kirch_rmatvec = jit(**numba_opts)(self._amp_kirch_rmatvec)
             elif self.travsrcrec:
                 self._kirch_matvec = jit(**numba_opts)(self._travsrcrec_kirch_matvec)
                 self._kirch_rmatvec = jit(**numba_opts)(self._travsrcrec_kirch_rmatvec)
@@ -1245,9 +987,6 @@ def _register_multiplications(self, engine: str) -> None:
             if self.dynamic and self.travsrcrec:
                 self._kirch_matvec = self._ampsrcrec_kirch_matvec
                 self._kirch_rmatvec = self._ampsrcrec_kirch_rmatvec
-            elif self.dynamic and not self.travsrcrec:
-                self._kirch_matvec = self._amp_kirch_matvec
-                self._kirch_rmatvec = self._amp_kirch_rmatvec
             elif self.travsrcrec:
                 self._kirch_matvec = self._travsrcrec_kirch_matvec
                 self._kirch_rmatvec = self._travsrcrec_kirch_rmatvec
@@ -1270,32 +1009,8 @@ def _matvec(self, x: NDArray) -> NDArray:
                 self.vel,
                 self.trav_srcs,
                 self.trav_recs,
-                self.dist_srcs,
-                self.dist_recs,
-                self.aperture[0],
-                self.aperture[1],
-                self.aperturetap,
-                self.nz,
-                self.six,
-                self.rix,
-                self.angleaperture[0],
-                self.angleaperture[1],
-                self.angle_srcs,
-                self.angle_recs,
-                self.snell[0],
-                self.snell[1],
-                self.maxdist,
-            )
-        elif self.dynamic and not self.travsrcrec:
-            inputs = (
-                x.ravel(),
-                y,
-                self.nsnr,
-                self.nt,
-                self.ni,
-                self.itrav,
-                self.travd,
-                self.amp,
+                self.amp_srcs,
+                self.amp_recs,
                 self.aperture[0],
                 self.aperture[1],
                 self.aperturetap,
@@ -1306,10 +1021,8 @@ def _matvec(self, x: NDArray) -> NDArray:
                 self.angleaperture[1],
                 self.angle_srcs,
                 self.angle_recs,
-                self.snell[0],
-                self.snell[1],
             )
-        elif not self.dynamic and self.travsrcrec:
+        elif self.travsrcrec:
             inputs = (
                 x.ravel(),
                 y,
@@ -1321,7 +1034,7 @@ def _matvec(self, x: NDArray) -> NDArray:
                 self.trav_srcs,
                 self.trav_recs,
             )
-        elif not self.dynamic and not self.travsrcrec:
+        elif not self.travsrcrec:
             inputs = (x.ravel(), y, self.nsnr, self.nt, self.ni, self.itrav, self.travd)
 
         y = self._kirch_matvec(*inputs)
@@ -1345,32 +1058,8 @@ def _rmatvec(self, x: NDArray) -> NDArray:
                 self.vel,
                 self.trav_srcs,
                 self.trav_recs,
-                self.dist_srcs,
-                self.dist_recs,
-                self.aperture[0],
-                self.aperture[1],
-                self.aperturetap,
-                self.nz,
-                self.six,
-                self.rix,
-                self.angleaperture[0],
-                self.angleaperture[1],
-                self.angle_srcs,
-                self.angle_recs,
-                self.snell[0],
-                self.snell[1],
-                self.maxdist,
-            )
-        elif self.dynamic and not self.travsrcrec:
-            inputs = (
-                x,
-                y,
-                self.nsnr,
-                self.nt,
-                self.ni,
-                self.itrav,
-                self.travd,
-                self.amp,
+                self.amp_srcs,
+                self.amp_recs,
                 self.aperture[0],
                 self.aperture[1],
                 self.aperturetap,
@@ -1381,10 +1070,8 @@ def _rmatvec(self, x: NDArray) -> NDArray:
                 self.angleaperture[1],
                 self.angle_srcs,
                 self.angle_recs,
-                self.snell[0],
-                self.snell[1],
             )
-        elif not self.dynamic and self.travsrcrec:
+        elif self.travsrcrec:
             inputs = (
                 x,
                 y,
@@ -1396,7 +1083,7 @@ def _rmatvec(self, x: NDArray) -> NDArray:
                 self.trav_srcs,
                 self.trav_recs,
             )
-        elif not self.dynamic and not self.travsrcrec:
+        elif not self.travsrcrec:
             inputs = (x, y, self.nsnr, self.nt, self.ni, self.itrav, self.travd)
 
         y = self._kirch_rmatvec(*inputs)
diff --git a/pytests/test_kirchhoff.py b/pytests/test_kirchhoff.py
index 317821c3..b210a95b 100755
--- a/pytests/test_kirchhoff.py
+++ b/pytests/test_kirchhoff.py
@@ -287,7 +287,7 @@ def test_kirchhoff3d(par):
 )
 def test_kirchhoff2d_trav_vs_travsrcrec(par):
     """Compare 2D Kirchhoff operator forward and adjoint when using trav (original behavior)
-    or trav_src and trav_rec (new reccomended behaviour)"""
+    or trav_src and trav_rec (new reccommended behaviour)"""
 
     # new behaviour
     Dop = Kirchhoff(
@@ -311,21 +311,7 @@ def test_kirchhoff2d_trav_vs_travsrcrec(par):
     ) + Dop.trav_recs.reshape(PAR["nx"] * PAR["nz"], 1, PAR["nrx"])
     trav = trav.reshape(PAR["nx"] * PAR["nz"], PAR["nsx"] * PAR["nrx"])
     if par["dynamic"]:
-        dist = Dop.dist_srcs.reshape(
-            PAR["nx"] * PAR["nz"], PAR["nsx"], 1
-        ) + Dop.dist_recs.reshape(PAR["nx"] * PAR["nz"], 1, PAR["nrx"])
-        dist = dist.reshape(PAR["nx"] * PAR["nz"], PAR["nsx"] * PAR["nrx"])
-
-        cosangle = np.cos(Dop.angle_srcs).reshape(
-            PAR["nx"] * PAR["nz"], PAR["nsx"], 1
-        ) + np.cos(Dop.angle_recs).reshape(PAR["nx"] * PAR["nz"], 1, PAR["nrx"])
-        cosangle = cosangle.reshape(PAR["nx"] * PAR["nz"], PAR["nsx"] * PAR["nrx"])
-
-        epsdist = 1e-2
-        amp = 1 / (dist + epsdist * np.max(dist))
-
-        amp *= np.abs(cosangle)
-        amp /= v0
+        amp = (Dop.amp_srcs, Dop.amp_recs)
 
     D1op = Kirchhoff(
         z,

From 0e5fb0b9fbeef25649d3ff78a76df08c66753831 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sat, 30 Sep 2023 23:10:15 +0300
Subject: [PATCH 02/16] tmp: provisionally force setuptools_scm<8.0.0

---
 requirements-dev.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/requirements-dev.txt b/requirements-dev.txt
index 8eabb62c..59bce051 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -12,7 +12,7 @@ matplotlib
 ipython
 pytest
 pytest-runner
-setuptools_scm
+setuptools_scm<8.0.0
 docutils<0.18
 Sphinx
 pydata-sphinx-theme

From c586d49da39ae67e56bff637e83ce5216a4e2e15 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sat, 30 Sep 2023 23:42:30 +0300
Subject: [PATCH 03/16] build: switch from setup.py to pyproject.toml

---
 pyproject.toml       | 47 ++++++++++++++++++++++++++++++++++++++
 requirements-dev.txt |  2 +-
 setup.py             | 54 --------------------------------------------
 3 files changed, 48 insertions(+), 55 deletions(-)
 create mode 100644 pyproject.toml
 delete mode 100755 setup.py

diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 00000000..173ac171
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,47 @@
+[build-system]
+requires = [
+    "setuptools>=60",
+    "setuptools-scm>=8.0"]
+
+[project]
+name = "pylops"
+authors = [
+    {name = "Matteo Ravasi", email = "matteoravasi@gmail.com"},
+]
+description = "Python library implementing linear operators to allow solving large-scale optimization problems without requiring to explicitly create a dense (or sparse) matrix."
+readme = "README.md"
+requires-python = ">=3.7"
+keywords = ["algebra", "inverse problems", "large-scale optimization"]
+license.file = "LICENSE.md"
+classifiers = [
+        "Development Status :: 5 - Production/Stable",
+        "Intended Audience :: Developers",
+        "Intended Audience :: Science/Research",
+        "License :: OSI Approved :: GNU Lesser General Public License v2 or later (LGPLv2+)",
+        "Natural Language :: English",
+        "Programming Language :: Python :: 3.8",
+        "Programming Language :: Python :: 3.9",
+        "Topic :: Scientific/Engineering :: Mathematics",
+    ]
+dynamic = ["version", ]
+dependencies = [
+        "numpy >= 1.21.0",
+        "scipy >= 1.4.0",
+        "pylops >= 1.17.0",
+        "setuptools >= 61.2.0",
+    ]
+[project.optional-dependencies]
+advanced = ["llvmlite",
+            "numba",
+            "pyfftw",
+            "PyWavelets",
+            "scikit-fmm",
+            "spgl1",
+]
+
+[tool.setuptools.packages.find]
+where = ["pylops"]
+namespaces = false
+
+[tool.setuptools_scm]
+version_file = "pylops/version.py"
diff --git a/requirements-dev.txt b/requirements-dev.txt
index 59bce051..8eabb62c 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -12,7 +12,7 @@ matplotlib
 ipython
 pytest
 pytest-runner
-setuptools_scm<8.0.0
+setuptools_scm
 docutils<0.18
 Sphinx
 pydata-sphinx-theme
diff --git a/setup.py b/setup.py
deleted file mode 100755
index 8b82afa2..00000000
--- a/setup.py
+++ /dev/null
@@ -1,54 +0,0 @@
-import os
-
-from setuptools import find_packages, setup
-
-
-def src(pth):
-    return os.path.join(os.path.dirname(__file__), pth)
-
-
-# Project description
-descr = (
-    "Python library implementing linear operators to allow solving large-scale optimization "
-    "problems without requiring to explicitly create a dense (or sparse) matrix."
-)
-
-# Setup
-setup(
-    name="pylops",
-    description=descr,
-    long_description=open(src("README.md")).read(),
-    long_description_content_type="text/markdown",
-    keywords=["algebra", "inverse problems", "large-scale optimization"],
-    classifiers=[
-        "Development Status :: 5 - Production/Stable",
-        "Intended Audience :: Developers",
-        "Intended Audience :: Science/Research",
-        "License :: OSI Approved :: GNU Lesser General Public License v2 or later (LGPLv2+)",
-        "Natural Language :: English",
-        "Programming Language :: Python :: 3.8",
-        "Programming Language :: Python :: 3.9",
-        "Topic :: Scientific/Engineering :: Mathematics",
-    ],
-    author="mrava",
-    author_email="matteoravasi@gmail.com",
-    install_requires=["numpy >= 1.21.0", "scipy >= 1.4.0"],
-    extras_require={
-        "advanced": [
-            "llvmlite",
-            "numba",
-            "pyfftw",
-            "PyWavelets",
-            "scikit-fmm",
-            "spgl1",
-        ]
-    },
-    packages=find_packages(exclude=["pytests"]),
-    use_scm_version=dict(
-        root=".", relative_to=__file__, write_to=src("pylops/version.py")
-    ),
-    setup_requires=["pytest-runner", "setuptools_scm"],
-    test_suite="pytests",
-    tests_require=["pytest"],
-    zip_safe=True,
-)

From bdd79d59e3fb2b97da3805e756b520cd4d1cc859 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sat, 30 Sep 2023 23:55:34 +0300
Subject: [PATCH 04/16] minor: fixed pyproject.toml

---
 pyproject.toml | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 173ac171..d525eb86 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -23,13 +23,14 @@ classifiers = [
         "Programming Language :: Python :: 3.9",
         "Topic :: Scientific/Engineering :: Mathematics",
     ]
-dynamic = ["version", ]
+dynamic = ["version"]
 dependencies = [
         "numpy >= 1.21.0",
         "scipy >= 1.4.0",
         "pylops >= 1.17.0",
         "setuptools >= 61.2.0",
     ]
+
 [project.optional-dependencies]
 advanced = ["llvmlite",
             "numba",
@@ -41,6 +42,7 @@ advanced = ["llvmlite",
 
 [tool.setuptools.packages.find]
 where = ["pylops"]
+exclude = ["pytests"]
 namespaces = false
 
 [tool.setuptools_scm]

From 40e28cf6687aa04574fd878bf41ddc5cc776841a Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sun, 1 Oct 2023 08:55:44 +0300
Subject: [PATCH 05/16] minor: simplify pyproject.toml

---
 .github/workflows/build.yaml |  2 +-
 azure-pipelines.yml          |  4 ++--
 pyproject.toml               | 19 +++++--------------
 3 files changed, 8 insertions(+), 17 deletions(-)

diff --git a/.github/workflows/build.yaml b/.github/workflows/build.yaml
index f6852e4c..63b4a4e3 100644
--- a/.github/workflows/build.yaml
+++ b/.github/workflows/build.yaml
@@ -7,7 +7,7 @@ jobs:
     strategy:
       matrix:
         platform: [ ubuntu-latest, macos-latest ]
-        python-version: ["3.7", "3.8", "3.9", "3.10"]
+        python-version: ["3.8", "3.9", "3.10"]
 
     runs-on: ${{ matrix.platform }}
     steps:
diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index d778e70a..661541a1 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -55,7 +55,7 @@ jobs:
   steps:
   - task: UsePythonVersion@0
     inputs:
-      versionSpec: '3.7'
+      versionSpec: '3.8'
       architecture: 'x64'
 
   - script: |
@@ -84,7 +84,7 @@ jobs:
   steps:
   - task: UsePythonVersion@0
     inputs:
-      versionSpec: '3.7'
+      versionSpec: '3.8'
       architecture: 'x64'
 
   - script: |
diff --git a/pyproject.toml b/pyproject.toml
index d525eb86..1b36d1f1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,16 +1,15 @@
 [build-system]
-requires = [
-    "setuptools>=60",
-    "setuptools-scm>=8.0"]
+requires = ["setuptools", "setuptools-scm"]
+build-backend = "setuptools.build_meta"
 
 [project]
 name = "pylops"
 authors = [
     {name = "Matteo Ravasi", email = "matteoravasi@gmail.com"},
 ]
-description = "Python library implementing linear operators to allow solving large-scale optimization problems without requiring to explicitly create a dense (or sparse) matrix."
+description = "Python library implementing linear operators to allow solving large-scale optimization problems."
 readme = "README.md"
-requires-python = ">=3.7"
+requires-python = ">=3.8"
 keywords = ["algebra", "inverse problems", "large-scale optimization"]
 license.file = "LICENSE.md"
 classifiers = [
@@ -23,13 +22,13 @@ classifiers = [
         "Programming Language :: Python :: 3.9",
         "Topic :: Scientific/Engineering :: Mathematics",
     ]
-dynamic = ["version"]
 dependencies = [
         "numpy >= 1.21.0",
         "scipy >= 1.4.0",
         "pylops >= 1.17.0",
         "setuptools >= 61.2.0",
     ]
+dynamic = ["version"]
 
 [project.optional-dependencies]
 advanced = ["llvmlite",
@@ -39,11 +38,3 @@ advanced = ["llvmlite",
             "scikit-fmm",
             "spgl1",
 ]
-
-[tool.setuptools.packages.find]
-where = ["pylops"]
-exclude = ["pytests"]
-namespaces = false
-
-[tool.setuptools_scm]
-version_file = "pylops/version.py"

From 8913d5678b281b2db6c035fb1a2b0f661fb921bd Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sun, 1 Oct 2023 09:17:06 +0300
Subject: [PATCH 06/16] minor: playing with dynamic version setting

---
 pyproject.toml | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 1b36d1f1..e10a9027 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -28,7 +28,7 @@ dependencies = [
         "pylops >= 1.17.0",
         "setuptools >= 61.2.0",
     ]
-dynamic = ["version"]
+dynamic = ["version", ]
 
 [project.optional-dependencies]
 advanced = ["llvmlite",
@@ -38,3 +38,10 @@ advanced = ["llvmlite",
             "scikit-fmm",
             "spgl1",
 ]
+
+[tool.setuptools.packages.find]
+exclude = ["pytests"]
+namespaces = false
+
+[tool.setuptools_scm]
+version_file = "pylops/version.py"

From d332e6338afd1942ae83fba82c9ef609ccba3097 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sun, 1 Oct 2023 09:36:41 +0300
Subject: [PATCH 07/16] minor: remove pylops from dependencies (added by
 mistake)

---
 pyproject.toml | 2 --
 1 file changed, 2 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index e10a9027..90d8e517 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -25,8 +25,6 @@ classifiers = [
 dependencies = [
         "numpy >= 1.21.0",
         "scipy >= 1.4.0",
-        "pylops >= 1.17.0",
-        "setuptools >= 61.2.0",
     ]
 dynamic = ["version", ]
 

From 4bbd0019ca98f0aa10efcc0dc4c6d540ce33605c Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sun, 1 Oct 2023 09:53:32 +0300
Subject: [PATCH 08/16] minor: moved metadata to setup.cfg

---
 pyproject.toml | 43 +---------------------------------------
 setup.cfg      | 54 +++++++++++++++++++++++++++++++++++++++++++++++++-
 2 files changed, 54 insertions(+), 43 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index 90d8e517..727c79dd 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,45 +1,4 @@
 [build-system]
 requires = ["setuptools", "setuptools-scm"]
 build-backend = "setuptools.build_meta"
-
-[project]
-name = "pylops"
-authors = [
-    {name = "Matteo Ravasi", email = "matteoravasi@gmail.com"},
-]
-description = "Python library implementing linear operators to allow solving large-scale optimization problems."
-readme = "README.md"
-requires-python = ">=3.8"
-keywords = ["algebra", "inverse problems", "large-scale optimization"]
-license.file = "LICENSE.md"
-classifiers = [
-        "Development Status :: 5 - Production/Stable",
-        "Intended Audience :: Developers",
-        "Intended Audience :: Science/Research",
-        "License :: OSI Approved :: GNU Lesser General Public License v2 or later (LGPLv2+)",
-        "Natural Language :: English",
-        "Programming Language :: Python :: 3.8",
-        "Programming Language :: Python :: 3.9",
-        "Topic :: Scientific/Engineering :: Mathematics",
-    ]
-dependencies = [
-        "numpy >= 1.21.0",
-        "scipy >= 1.4.0",
-    ]
-dynamic = ["version", ]
-
-[project.optional-dependencies]
-advanced = ["llvmlite",
-            "numba",
-            "pyfftw",
-            "PyWavelets",
-            "scikit-fmm",
-            "spgl1",
-]
-
-[tool.setuptools.packages.find]
-exclude = ["pytests"]
-namespaces = false
-
-[tool.setuptools_scm]
-version_file = "pylops/version.py"
+x
diff --git a/setup.cfg b/setup.cfg
index 1691b3bc..67636c69 100755
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,3 +1,56 @@
+[metadata]
+name = pylops
+fullname = PyLops
+description = Python library implementing linear operators to allow solving large-scale optimization problems
+long_description = file: README.md
+long_description_content_type = text/markdown
+author = The PyLops Development Team
+author_email = matteoravasi@gmail.com
+maintainer = "Matteo Ravasi"
+maintainer_email = matteoravasi@gmail.com
+license = LGPL-3.0 License
+license_file = LICENSE.md
+platform = any
+keywords = algebra, inverse problems, large-scale optimization
+classifiers =
+    Development Status :: 5 - Production/Stable
+    Intended Audience :: Developers
+    Intended Audience :: Science/Research
+    Intended Audience :: Education
+    License :: OSI Approved :: GNU Lesser General Public License v3 (LGPLv3),
+    Natural Language :: English
+    Operating System :: OS Independent
+    Programming Language :: Python :: 3 :: Only
+    Programming Language :: Python :: 3.8
+    Programming Language :: Python :: 3.9
+    Programming Language :: Python :: 3.10
+    Topic :: Scientific/Engineering : Mathematics
+url = https://github.com/pylops/pylops
+
+project_urls =
+    Documentation = https://pylops.readthedocs.io/
+    Release Notes = https://github.com/pylops/pylops/releases
+    Bug Tracker = https://github.com/pylops/pylops/issues
+    Source Code = https://github.com/pylops/pylops
+
+[options]
+zip_safe = True
+include_package_data = True
+packages = find:
+python_requires = >=3.8
+install_requires =
+    numpy >= 1.21.0
+    scipy >= 1.4.0
+
+[options.extras_require]
+advanced =
+    llvmlite
+    numba
+    pyfftw
+    PyWavelets
+    scikit-fmm
+    spgl1
+
 [aliases]
 test=pytest
 
@@ -11,7 +64,6 @@ per-file-ignores =
     __init__.py: F401, F403, F405
 max-line-length = 88
 
-
 # mypy global options
 [mypy]
 plugins = numpy.typing.mypy_plugin

From d206560e56f628143d0e76f0e821678d30714493 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sun, 1 Oct 2023 09:54:36 +0300
Subject: [PATCH 09/16] minor: change way to run tests in azure pipelines

---
 azure-pipelines.yml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index 661541a1..ef5306bb 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -65,7 +65,7 @@ jobs:
     displayName: 'Install prerequisites and library'
 
   - script: |
-      python setup.py test
+      pytest
     condition: succeededOrFailed()
     displayName: 'Run tests'
 
@@ -94,6 +94,6 @@ jobs:
     displayName: 'Install prerequisites and library'
 
   - script: |
-      python setup.py test
+      pytest
     condition: succeededOrFailed()
     displayName: 'Run tests'

From 541c12f772f2bd890b35062e5114c0290d7cfb99 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sun, 1 Oct 2023 10:08:43 +0300
Subject: [PATCH 10/16] minor: fix typo in pyproject.toml

---
 pyproject.toml | 1 -
 1 file changed, 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 727c79dd..8fd8d67e 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,4 +1,3 @@
 [build-system]
 requires = ["setuptools", "setuptools-scm"]
 build-backend = "setuptools.build_meta"
-x

From a37f0ef159861375794d8cad764799c7d9861670 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sun, 1 Oct 2023 10:21:55 +0300
Subject: [PATCH 11/16] minor: move from setuptools to pip in rdt file

---
 .readthedocs.yaml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index 146206bb..a74e3480 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -19,5 +19,5 @@ sphinx:
 python:
   install:
     - requirements: requirements-dev.txt
-    - method: setuptools
+    - method: pip
       path: .

From 70fdaae692415078ea90cdf8d706fa71816a2562 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sun, 1 Oct 2023 10:34:15 +0300
Subject: [PATCH 12/16] minor: fix coverage action

---
 .github/workflows/codacy-coverage-reporter.yaml | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/.github/workflows/codacy-coverage-reporter.yaml b/.github/workflows/codacy-coverage-reporter.yaml
index b0a8f90d..3e390ee5 100644
--- a/.github/workflows/codacy-coverage-reporter.yaml
+++ b/.github/workflows/codacy-coverage-reporter.yaml
@@ -17,8 +17,9 @@ jobs:
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
-        pip install flake8 pytest coverage
+        pip install flake8 pytest
         if [ -f requirements.txt ]; then pip install -r requirements-dev.txt; fi
+        pip install coverage
     - name: Install pylops
       run: |
         pip install .

From 99ab70b365c6d5ff86a0841838c683ceb5af0017 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sun, 1 Oct 2023 10:55:22 +0300
Subject: [PATCH 13/16] minor: test coverage action

---
 .github/workflows/codacy-coverage-reporter.yaml | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/.github/workflows/codacy-coverage-reporter.yaml b/.github/workflows/codacy-coverage-reporter.yaml
index 3e390ee5..6932e565 100644
--- a/.github/workflows/codacy-coverage-reporter.yaml
+++ b/.github/workflows/codacy-coverage-reporter.yaml
@@ -13,13 +13,12 @@ jobs:
     - name: Set up Python
       uses: actions/setup-python@v4
       with:
-        python-version: 3.9
+        python-version: 3.8
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
         pip install flake8 pytest
         if [ -f requirements.txt ]; then pip install -r requirements-dev.txt; fi
-        pip install coverage
     - name: Install pylops
       run: |
         pip install .

From 14d8fc2a1feb7ee888682808c88eea5ae174155d Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sun, 1 Oct 2023 20:09:56 +0300
Subject: [PATCH 14/16] fix: trying to fix coverage

---
 .github/workflows/codacy-coverage-reporter.yaml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/codacy-coverage-reporter.yaml b/.github/workflows/codacy-coverage-reporter.yaml
index 6932e565..837e1478 100644
--- a/.github/workflows/codacy-coverage-reporter.yaml
+++ b/.github/workflows/codacy-coverage-reporter.yaml
@@ -17,7 +17,7 @@ jobs:
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
-        pip install flake8 pytest
+        pip install flake8 pytest coverage
         if [ -f requirements.txt ]; then pip install -r requirements-dev.txt; fi
     - name: Install pylops
       run: |

From d69d238c46caceecd26451fec64221668bf77953 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sun, 1 Oct 2023 20:15:18 +0300
Subject: [PATCH 15/16] fix: testing empty fix coverage

---
 .github/workflows/codacy-coverage-reporter.yaml | 12 +++---------
 1 file changed, 3 insertions(+), 9 deletions(-)

diff --git a/.github/workflows/codacy-coverage-reporter.yaml b/.github/workflows/codacy-coverage-reporter.yaml
index 837e1478..188c7574 100644
--- a/.github/workflows/codacy-coverage-reporter.yaml
+++ b/.github/workflows/codacy-coverage-reporter.yaml
@@ -17,17 +17,11 @@ jobs:
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
-        pip install flake8 pytest coverage
+        pip install flake8 pytest
         if [ -f requirements.txt ]; then pip install -r requirements-dev.txt; fi
     - name: Install pylops
       run: |
         pip install .
-    - name: Code coverage with coverage
+    - name: Test with pytest
       run: |
-        coverage run -m pytest
-        coverage xml
-    - name: Run codacy-coverage-reporter
-      uses: codacy/codacy-coverage-reporter-action@v1
-      with:
-        project-token: ${{ secrets.CODACY_PROJECT_TOKEN }}
-        coverage-reports: coverage.xml
+        pytest

From 0f73b7997439484a0d4924e17cf4e8a90e5b4442 Mon Sep 17 00:00:00 2001
From: mrava87 <matteoravasi@gmail.com>
Date: Sun, 1 Oct 2023 20:31:06 +0300
Subject: [PATCH 16/16] fix: testing empty fix coverage

---
 .github/workflows/codacy-coverage-reporter.yaml | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/.github/workflows/codacy-coverage-reporter.yaml b/.github/workflows/codacy-coverage-reporter.yaml
index 188c7574..837e1478 100644
--- a/.github/workflows/codacy-coverage-reporter.yaml
+++ b/.github/workflows/codacy-coverage-reporter.yaml
@@ -17,11 +17,17 @@ jobs:
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
-        pip install flake8 pytest
+        pip install flake8 pytest coverage
         if [ -f requirements.txt ]; then pip install -r requirements-dev.txt; fi
     - name: Install pylops
       run: |
         pip install .
-    - name: Test with pytest
+    - name: Code coverage with coverage
       run: |
-        pytest
+        coverage run -m pytest
+        coverage xml
+    - name: Run codacy-coverage-reporter
+      uses: codacy/codacy-coverage-reporter-action@v1
+      with:
+        project-token: ${{ secrets.CODACY_PROJECT_TOKEN }}
+        coverage-reports: coverage.xml