diff --git a/w-stacking-fftw.c b/w-stacking-fftw.c index f13b3c188bee33b4dc5769f8888f4b63de73805b..ef56aa48017689fca0a30fe290593ee7def22a55 100644 --- a/w-stacking-fftw.c +++ b/w-stacking-fftw.c @@ -643,8 +643,8 @@ if(rank == 0){ dx, dw, w_support, - xaxis, - yaxis, + xaxis, + yaxis, gridss, num_threads); diff --git a/w-stacking.cu b/w-stacking.cu index d5e8e8ca5aa849f7601ae2566f0ec799b558c3da..8728fdab1bb6a374e81fa91a728aa90c2a5baba1 100644 --- a/w-stacking.cu +++ b/w-stacking.cu @@ -22,6 +22,22 @@ gauss_kernel_norm(double norm, double std22, double u_dist, double v_dist) return conv_weight; } +void makeGaussKernel(double * kernel, + int KernelLen, + double std22) +{ + + double norm = std22/PI; + int n = 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); + } + + for (int i = 0; i != mid; i++) kernel[i] = kernel[n - 1 - i]; + +} + // Kaiser-Bessel Kernel: it is adapted from WSClean double bessel0(double x, double precision) { // Calculate I_0 = SUM of m 0 -> inf [ (x/2)^(2m) ] @@ -188,7 +204,9 @@ void wstack( 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); + // Loop over visibilities. // Switch between CUDA and GPU versions #ifdef __CUDACC__ @@ -309,11 +327,14 @@ 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); - double conv_weight = gauss_kernel_norm(norm,std22,u_dist,v_dist); + //double conv_weight = gauss_kernel_norm(norm,std22,u_dist,v_dist); + double conv_weight = convkernel[jKer]*convkernel[kKer]; // Loops over frequencies and polarizations double add_term_real = 0.0; double add_term_img = 0.0;