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

better plane selection algorithm

parent 318424a6
......@@ -1503,13 +1503,14 @@ cout << "data w range: " << wmin << " to " << wmax << endl;
wmax += 0.5*w_supp*dw;
size_t nplanes = size_t((wmax-wmin)/dw)+2;
cout << "nplanes: " << nplanes << endl;
double dwmax=0.5*w_supp*dw;
vector<size_t> nvis_plane(nplanes,0);
vector<int> minplane(nvis);
for (size_t ipart=0; ipart<nvis; ++ipart)
{
int iplane = int(wval[ipart]-wmin)/dw;
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];
int plane0 = int((wval[ipart]-wmin)/dw+0.5*w_supp)-int(w_supp-1);
minplane[ipart]=plane0;
for (int i=max<int>(0,plane0); i<min<int>(nplanes,plane0+w_supp); ++i)
++nvis_plane[i];
}
auto accum_ = makeArray<complex<T>>({nx_dirty, ny_dirty});
......@@ -1528,10 +1529,10 @@ cout << "working on w plane #" << iw << " containing " << nvis_plane[iw]
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)
if ((iw>=minplane[ipart]) && (iw<minplane[ipart]+w_supp))
{
myassert(cnt<nvis_plane[iw],"must not happen");
double x=2./(w_supp*dw)*abs(wcur-wval[ipart]);
double x=min(1.,2./(w_supp*dw)*abs(wcur-wval[ipart]));
vis_loc[cnt] = vis[ipart]*exp(beta*(sqrt(1.-x*x)-1.));
idx_loc[cnt] = idx[ipart];
++cnt;
......
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