diff --git a/cygrid/cygrid.pyx b/cygrid/cygrid.pyx index 6f38ba2..52b2b69 100644 --- a/cygrid/cygrid.pyx +++ b/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, False, # ('kernel_sigma', 'param_a', 'param_b') ), + 'gaussbessel': ( + '1D Gauss-Bessel kernel', 3, 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 = True self.kernel_params_arr = kernel_params_arr diff --git a/cygrid/kernels.pxd b/cygrid/kernels.pxd index 705dc19..2f496c5 100644 --- a/cygrid/kernels.pxd +++ b/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 diff --git a/cygrid/kernels.pyx b/cygrid/kernels.pyx index 7c78961..3d081a7 100644 --- a/cygrid/kernels.pyx +++ b/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 diff --git a/pyproject.toml b/pyproject.toml index 489da7e..d092bb2 100644 --- a/pyproject.toml +++ b/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"