Tabulated gaussian kernel with increased precision added (only CPU version)

Makefile

@@ -40,6 +40,8 @@ OPT += -DWRITE_IMAGE
 #OPT += -DPARALLELIO
 # Normalize uvw in case it is not done in the binMS
 OPT += -DNORMALIZE_UVW
+# Gridding kernel: GAUSS, GAUSS_HI_PRECISION
+OPT += -DGAUSS_HI_PRECISION
 
 ifeq (FITSIO,$(findstring FITSIO,$(OPT)))
 	LIBS += -L$(FITSIO_LIB) -lcfitsio

w-stacking.cu

@@ -24,11 +24,12 @@ gauss_kernel_norm(double norm, double std22, double u_dist, double v_dist)
 
 void makeGaussKernel(double * kernel,
 		     int KernelLen,
+		     int increaseprecision,
 		     double std22)
 {
 
   double norm = std22/PI;
-  int n = KernelLen, mid = n / 2;
+  int n = increaseprecision*KernelLen, mid = n / 2;
   for (int i = 0; i != mid + 1; i++) {
       double term = (double)i / mid;
       kernel[mid + i] = sqrt(norm) * exp(-(term*term)*std22);
@@ -201,11 +202,12 @@ void wstack(
     // initialize the convolution kernel
     // gaussian:
     int KernelLen = (w_support-1)/2;
+    int increaseprecision = 21; // this number must be odd: increaseprecison*w_support must be odd (w_support must be odd)
     double std = 1.0;
     double std22 = 1.0/(2.0*std*std);
     double norm = std22/PI;
-    double * convkernel = malloc(w_support*sizeof(*convkernel));
-    makeGaussKernel(convkernel,w_support,std22);
+    double * convkernel = malloc(increaseprecision*w_support*sizeof(*convkernel));
+    makeGaussKernel(convkernel,w_support,increaseprecision,std22);
     
     // Loop over visibilities.
 // Switch between CUDA and GPU versions
@@ -327,14 +329,17 @@ void wstack(
 
             for (j = jmin; j <= jmax; j++)
             {
-		int jKer = (int)j-pos_u+KernelLen;
-		int kKer = (int)k-pos_v+KernelLen;
                 double u_dist = (double)j+0.5 - pos_u;
 		long iKer = 2 * (j + k*grid_size_x + grid_w*grid_size_x*grid_size_y);
+                int jKer = (int)(increaseprecision*u_dist) + KernelLen;
+                int kKer = (int)(increaseprecision*v_dist) + KernelLen;
 
-
-		//double conv_weight = gauss_kernel_norm(norm,std22,u_dist,v_dist);
+                #ifdef GAUSS_HI_PRECISION 
+		double conv_weight = gauss_kernel_norm(norm,std22,u_dist,v_dist);
+                #endif
+                #ifdef GAUSS
 		double conv_weight = convkernel[jKer]*convkernel[kKer];
+                #endif
 		// Loops over frequencies and polarizations
 		double add_term_real = 0.0;
 		double add_term_img = 0.0;