adds gauss-bessel kernel

cygrid/cygrid.pyx

@@ -94,7 +94,8 @@ __all__ = ['Cygrid', 'WcsGrid', 'SlGrid', 'ShapeError']
 
 
 from .kernels cimport (
-    gaussian_1d_kernel, gaussian_2d_kernel, tapered_sinc_1d_kernel
+    gaussian_1d_kernel, gaussian_2d_kernel, tapered_sinc_1d_kernel,
+    gauss_bessel_kernel
     )
 from .hphashtab cimport HpxHashTable
 from .helpers cimport (
@@ -424,6 +425,10 @@ cdef class Cygrid(object):
                     '1D tapered-sinc kernel', 3, <bint> False,
                     # ('kernel_sigma', 'param_a', 'param_b')
                     ),
+                'gaussbessel': (
+                    '1D Gauss-Bessel kernel', 3, <bint> False,
+                    # ('kernel_sigma', 'param_a', 'param_b')
+                    ),
                 }
 
         try:
@@ -456,6 +461,8 @@ cdef class Cygrid(object):
             self.kernel_func_ptr = gaussian_2d_kernel
         elif kernel_type_u == 'tapered_sinc':
             self.kernel_func_ptr = tapered_sinc_1d_kernel
+        elif kernel_type_u == 'gaussbessel':
+            self.kernel_func_ptr = gauss_bessel_kernel
 
         self.kernel_set = <bint> True
         self.kernel_params_arr = kernel_params_arr

cygrid/kernels.pxd

@@ -35,6 +35,7 @@ from numpy cimport (
     int8_t, int16_t, int32_t, int64_t, uint8_t, uint16_t,
     uint32_t, uint64_t, float32_t, float64_t
     )
+cimport scipy.special.cython_special as csc
 
 
 cdef:
@@ -50,4 +51,7 @@ cdef:
     double tapered_sinc_1d_kernel(
         double distance, double bearing, double[::1] kernel_params
         ) nogil
+    double gauss_bessel_kernel(
+        double distance, double bearing, double[::1] kernel_params
+        ) nogil
 

cygrid/kernels.pyx

@@ -145,3 +145,40 @@ cdef double tapered_sinc_1d_kernel(
     arg = PI * distance / kernel_params[0]
 
     return sinc(arg / kernel_params[2]) * exp(-(arg / kernel_params[1]) ** 2)
+
+
+cdef double gauss_bessel_kernel(
+        double distance, double bearing, double[::1] kernel_params
+        ) nogil:
+    '''
+    Gaussian tapered Bessel kernel following the definition
+    of Mangum et al. 2007 A&A 474 679M.
+
+    Parameters
+    ----------
+    distance : double
+        Radial distance/separation
+    bearing : double
+        unused - this is only in the call signature to allow function pointers
+    kernel_params : double[::1]
+        mem-view on kernel-parameters array
+        must contain:
+        kernel_params[0] : A third of the Half Power Beam Width
+        kernel_params[1] : 2.52
+        kernel_params[2] : 1.55
+
+    Returns
+    -------
+    Kernel weight : double
+    '''
+
+    cdef:
+        double arg, Earg, gauss, bessel
+
+    arg = PI * fabs(distance) / kernel_params[0] / kernel_params[2]
+    bessel = csc.j1(arg) / arg
+
+    Earg = distance * distance / kernel_params[0] / kernel_params[0] / kernel_params[1] / kernel_params[1] / 2.0
+    gauss = exp(-Earg)
+
+    return gauss * bessel

pyproject.toml

@@ -13,6 +13,7 @@ requires = [
     "numpy==1.23.4; python_version=='3.11'",
     # "oldest-supported-numpy",
     "extension-helpers",
+    "scipy"
     ]
 build-backend = "setuptools.build_meta"