Commit a26e60b3 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

more cleanup

parent 85c43c0e
......@@ -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)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment