add w correction

nifty_gridder.cc

@@ -1466,8 +1466,10 @@ template<typename T> pyarr_c<T> vis2dirty_wstack(const Baselines<T> &baselines,
   const GridderConfig<T> &gconf, const pyarr<uint32_t> &idx_,
   const pyarr<complex<T>> &vis_)
   {
-  auto nxdirty=gconf.Nxdirty();
-  auto nydirty=gconf.Nydirty();
+  auto nx_dirty=gconf.Nxdirty();
+  auto ny_dirty=gconf.Nydirty();
+  auto psx=gconf.Pixsize_x();
+  auto psy=gconf.Pixsize_y();
   checkArray(vis_, "vis", {0});
   size_t nvis = size_t(vis_.shape(0));
   checkArray(idx_, "idx", {nvis});
@@ -1508,7 +1510,7 @@ cout << "nplanes: " << nplanes << endl;
 cout << "nvis/plane: ";
 for (auto nv:nvis_plane) cout << nv << " ";
 cout << endl;
-  auto accum_ = makeArray<complex<T>>({nxdirty, nydirty});
+  auto accum_ = makeArray<complex<T>>({nx_dirty, ny_dirty});
   auto accum = accum_.mutable_data();
   for (ptrdiff_t i=0; i<accum_.size(); ++i)
     accum[i] = 0.;
@@ -1536,6 +1538,44 @@ cout << endl;
       accum[i] += dirty[i];
     }
   // correct for w gridding
+  {
+  py::gil_scoped_release release;
+
+  constexpr double pi = 3.141592653589793238462643383279502884197;
+  double x0 = -0.5*nx_dirty*psx,
+         y0 = -0.5*ny_dirty*psy;
+  auto beta = 2.3*w_supp;
+  auto p = int(1.5*w_supp+2);
+  vector<double> x, wgt;
+  legendre_prep(2*p,x,wgt);
+  auto psi = x;
+  for (auto &v:psi)
+    v = exp(beta*(sqrt(1-v*v)-1.));
+#pragma omp parallel num_threads(nthreads)
+{
+#pragma omp for schedule(static)
+  for (size_t i=0; i<=nx_dirty/2; ++i)
+    {
+    double fx = x0+i*psx;
+    fx *= fx;
+    for (size_t j=0; j<=ny_dirty/2; ++j)
+      {
+      double fy = y0+j*psy;
+      fy*=fy;
+      auto n=sqrt(1.-fx-fy);
+      double fct = 0.;
+      for (int i=0; i<p; ++i)
+        fct += wgt[i]*psi[i]*cos(pi*w_supp*n*x[i]);
+      fct = 1./(w_supp*fct);
+      size_t i2 = nx_dirty-i, j2 = ny_dirty-j;
+      accum[ny_dirty*i+j]*=fct;
+      accum[ny_dirty*i2+j]*=fct;
+      accum[ny_dirty*i2+j2]*=fct;
+      accum[ny_dirty*i+j2]*=fct;
+      }
+    }
+}
+  }
   return accum_;
   }