diff --git a/nifty_gridder.cc b/nifty_gridder.cc index 1828fbb42f1f9052472c70de0112a0a711eb2fd0..7d065c36798602e9072243a80ebdc2479b4e4402 100644 --- a/nifty_gridder.cc +++ b/nifty_gridder.cc @@ -1471,9 +1471,9 @@ template<typename T> pyarr_c<complex<T>> vis2dirty_wstack(const Baselines<T> &ba auto psx=gconf.Pixsize_x(); auto psy=gconf.Pixsize_y(); checkArray(vis_, "vis", {0}); + auto vis=vis_.template unchecked<1>(); size_t nvis = size_t(vis_.shape(0)); checkArray(idx_, "idx", {nvis}); - auto vis=vis_.template unchecked<1>(); auto idx=idx_.template unchecked<1>(); // determine w values for every visibility, and min/max w; @@ -1487,19 +1487,13 @@ template<typename T> pyarr_c<complex<T>> vis2dirty_wstack(const Baselines<T> &ba } cout << "data w range: " << wmin << " to " << wmax << endl; -//FIXME temporary -double dw; -{ double x0 = -0.5*nx_dirty*psx, y0 = -0.5*ny_dirty*psy; double nmin = sqrt(max(1.-x0*x0-y0*y0,0.))-1.; - dw = 0.25/abs(nmin); -} -cout << "delta w: " << dw << endl; + double dw = 0.25/abs(nmin); - double w_eps=1e-12; // FIXME + double w_eps=1e-7; // FIXME auto w_supp = get_supp(w_eps); -cout << "w_supp: " << w_supp << endl; auto beta=2.3*w_supp; wmin -= 0.5*w_supp*dw; wmax += 0.5*w_supp*dw; @@ -1513,16 +1507,15 @@ cout << "nplanes: " << nplanes << endl; for (int i=max<int>(0,iplane-w_supp); i<min<int>(nplanes,iplane+w_supp+1); ++i) if (abs(wval[ipart]-(wmin+i*dw))<dwmax) ++nvis_plane[i]; } -cout << "nvis/plane: "; -for (auto nv:nvis_plane) cout << nv << " "; -cout << endl; + 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.; for (size_t iw=0; iw<nplanes; ++iw) { -cout << "working on w plane #" << iw << endl; +cout << "working on w plane #" << iw << " containing " << nvis_plane[iw] + << " visibilities" << endl; if (nvis_plane[iw]==0) continue; double wcur = wmin+iw*dw; size_t cnt=0; @@ -1531,7 +1524,6 @@ cout << "working on w plane #" << iw << endl; pyarr<uint32_t> idx_loc_ = makeArray<uint32_t>({nvis_plane[iw]}); auto idx_loc = idx_loc_.mutable_data(); for (size_t ipart=0; ipart<nvis; ++ipart) -{ if (abs(wval[ipart]-wcur) < dwmax) { myassert(cnt<nvis_plane[iw],"must not happen"); @@ -1540,7 +1532,6 @@ myassert(cnt<nvis_plane[iw],"must not happen"); idx_loc[cnt] = idx[ipart]; ++cnt; } -} myassert(cnt==nvis_plane[iw],"must not happen 2"); auto dirty_ = gconf.apply_wscreen( gconf.grid2dirty_c(vis2grid_c(baselines, gconf, idx_loc_, vis_loc_, None, None)), wcur, false); @@ -1561,8 +1552,6 @@ cout << "applying correction for gridding in w direction" << endl; auto psi = x; for (auto &v:psi) v = exp(beta*(sqrt(1-v*v)-1.)); - double x0 = -0.5*nx_dirty*psx, - y0 = -0.5*ny_dirty*psy; #pragma omp parallel num_threads(nthreads) { #pragma omp for schedule(static)