Kaiser-Bessel and tabulated Gaussian available on CPU. Bug fixed in uvw normalization

Makefile

@@ -40,7 +40,7 @@ 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
+# Gridding kernel: GAUSS, GAUSS_HI_PRECISION, KAISERBESSEL
 OPT += -DGAUSS_HI_PRECISION
 
 ifeq (FITSIO,$(findstring FITSIO,$(OPT)))

w-stacking-fftw.c

@@ -354,12 +354,12 @@ if(rank == 0){
 	maxv_all = maxv;
 	maxw_all = maxw;
         #endif
-	double offset = 0.0;
-	double ming = MIN(minu_all,minv_all);
-	double maxg = MAX(maxu_all,maxv_all);
+	double offset = 0.001;
+	double ming = MAX(abs(minu_all),abs(minv_all));
+	double maxg = MAX(abs(maxu_all),abs(maxv_all));
+	maxg = MAX(maxg,ming);
 	minw = minw_all;
 	maxw = maxw_all;
-        ming = ming-offset*ming;
         maxg = maxg+offset*maxg;
         for (long inorm=0; inorm<Nmeasures; inorm++)
 	    {

w-stacking.cu

@@ -31,11 +31,12 @@ void makeGaussKernel(double * kernel,
   double norm = std22/PI;
   int n = increaseprecision*KernelLen, mid = n / 2;
   for (int i = 0; i != mid + 1; i++) {
-      double term = (double)i / mid;
+      double term = (double)i/(double)increaseprecision;
       kernel[mid + i] = sqrt(norm) * exp(-(term*term)*std22);
   }
 
   for (int i = 0; i != mid; i++) kernel[i] = kernel[n - 1 - i];
+//  for (int i = 0; i < n; i++) printf("%f\n",kernel[i]);
 
 }
 
@@ -53,10 +54,11 @@ double bessel0(double x, double precision) {
 }
 void makeKaiserBesselKernel(double * kernel,
 		            int KernelLen,
+			    int increaseprecision,
                             double alpha,
                             double overSamplingFactor,
                             int withSinc) {
-  int n = KernelLen, mid = n / 2;
+  int n = increaseprecision*KernelLen, mid = n / 2;
   double * sincKernel = malloc((mid + 1)*sizeof(*sincKernel));
   const double filterRatio = 1.0 / overSamplingFactor;
   sincKernel[0] = filterRatio;
@@ -73,6 +75,7 @@ void makeKaiserBesselKernel(double * kernel,
                 normFactor;
   }
   for (int i = 0; i != mid; i++) kernel[i] = kernel[n - 1 - i];
+  //for (int i = 0; i < n; i++) printf("%f\n",kernel[i]);
 }
 
 
@@ -202,12 +205,21 @@ 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)
+    int increaseprecision = 5; // 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(increaseprecision*w_support*sizeof(*convkernel));
+    double overSamplingFactor = 1.0;
+    int withSinc = 0;
+    double alpha = 8.6;
+    #ifdef GAUSS
     makeGaussKernel(convkernel,w_support,increaseprecision,std22);
+    #endif
+    #ifdef KAISERBESSEL
+    makeKaiserBesselKernel(convkernel, w_support, increaseprecision, alpha, overSamplingFactor, withSinc);
+    #endif
+
     
     // Loop over visibilities.
 // Switch between CUDA and GPU versions
@@ -325,20 +337,28 @@ void wstack(
         for (k = kmin; k <= kmax; k++)
         {
 
-            double v_dist = (double)k+0.5 - pos_v;
+            //double v_dist = (double)k+0.5 - pos_v;
+            double v_dist = (double)k - pos_v;
 
             for (j = jmin; j <= jmax; j++)
             {
-                double u_dist = (double)j+0.5 - pos_u;
+                //double u_dist = (double)j+0.5 - pos_u;
+                double u_dist = (double)j - 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;
+                int jKer = (int)(increaseprecision * (fabs(u_dist+(double)KernelLen)));
+                int kKer = (int)(increaseprecision * (fabs(v_dist+(double)KernelLen)));
 
                 #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];
+		//if(jKer < 0 || jKer >= 35 || kKer < 0 || kKer >= 35)
+		//	printf("%f %d %f %d\n",fabs(u_dist+(double)KernelLen),jKer,fabs(v_dist+(double)KernelLen),kKer);
+		//printf("%d %d %d %d %f %f %f %f %f\n",jKer, j, kKer, k, pos_u, pos_v, u_dist,v_dist,conv_weight);
+                #endif
+                #ifdef KAISERBESSEL
+		double conv_weight = convkernel[jKer]*convkernel[kKer];
 	        #endif	
 		// Loops over frequencies and polarizations
 		double add_term_real = 0.0;