Commit 4a97f909 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

better, but not OK yet

parent eacea104
......@@ -594,6 +594,9 @@ template<typename T> class GridderConfig
complex<T> wscreen(double x, double y, double w, bool adjoint) const
{
constexpr double pi = 3.141592653589793238462643383279502884197;
// PhD : sqrt (1-x*x-y*y)) - 1
//here 1 - sin(sqrt(x*x+y*y))**2/(1.+cos(sqrt(x*x+y*y)))
// == (cos(sqrt(x*x+y*y))**2 +cos(sqrt(x*x+y*y))) / (1.+cos(sqrt(x*x+y*y)))
double eps = sqrt(x+y);
double s = sin(eps);
double nm1 = -s*s/(1.+cos(eps));
......@@ -1539,12 +1542,14 @@ 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);
auto dirty = dirty_.data();
for (ptrdiff_t i=0; i<accum_.size(); ++i)
accum[i] += dirty[i];
}
// correct for w gridding
{
py::gil_scoped_release release;
......@@ -1565,16 +1570,21 @@ cout << "applying correction for gridding in w direction" << endl;
#pragma omp for schedule(static)
for (size_t i=0; i<=nx_dirty/2; ++i)
{
double fx = x0+i*psx;
// double fx = x0+i*psx;
double fx = -0.5+i*1./nx_dirty;
fx *= fx;
for (size_t j=0; j<=ny_dirty/2; ++j)
{
double fy = y0+j*psy;
// double fy = y0+j*psy;
double fy = -0.5+j*1./ny_dirty;
fy*=fy;
auto n=sqrt(1.-fx-fy)-1.;
// auto n=sqrt(1.-fx-fy)-1.;
auto n=cos(sqrt(fx+fy)) -1;
// cout << n << endl;
double fct = 0.;
for (int i=0; i<p; ++i)
fct += wgt[i]*psi[i]*cos(pi*w_supp*n*x[i]);
// cout << fct << endl;
fct = 1./(w_supp*fct);
if ((i==nx_dirty/2) &&(j==ny_dirty/2))
cout << i << " " << j << " " << fct << endl;
......
Supports Markdown
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