diff --git a/gridder_cxx.h b/gridder_cxx.h index cabe41c199c70a32e1c8f6fb109a03f9b8a108c3..0183ad9e384baeaf11592afd2e58d61680ed233a 100644 --- a/gridder_cxx.h +++ b/gridder_cxx.h @@ -295,6 +295,13 @@ template<typename T> struct UVW UVW (T u_, T v_, T w_) : u(u_), v(v_), w(w_) {} UVW operator* (T fct) const { return UVW(u*fct, v*fct, w*fct); } + void Flip() { u=-u; v=-v; w=-w; } + bool FixW() + { + bool flip = w<0; + if (flip) Flip(); + return flip; + } }; template<typename T> class Baselines @@ -314,7 +321,10 @@ template<typename T> class Baselines myassert(nrows*nchan<(size_t(1)<<32), "too many entries in MS"); f_over_c.resize(nchan); for (size_t i=0; i<nchan; ++i) + { + myassert(freq[i]>0, "negative channel frequency encountered"); f_over_c[i] = freq(i)/speedOfLight; + } coord.resize(nrows); for (size_t i=0; i<coord.size(); ++i) coord[i] = UVW<T>(coord_(i,0), coord_(i,1), coord_(i,2)); @@ -741,9 +751,11 @@ template<typename T, typename Serv> void x2grid_c for (size_t ipart=0; ipart<np; ++ipart) { UVW<T> coord = srv.getCoord(ipart); + auto flip = coord.FixW(); hlp.prep(coord.u, coord.v); auto * ptr = hlp.p0w; auto v(srv.getVis(ipart)*emb); + if (flip) v=conj(v); for (size_t cu=0; cu<supp; ++cu) { complex<T> tmp(v*ku[cu]); @@ -795,6 +807,7 @@ template<typename T, typename Serv> void grid2x_c for (size_t ipart=0; ipart<np; ++ipart) { UVW<T> coord = srv.getCoord(ipart); + auto flip = coord.FixW(); hlp.prep(coord.u, coord.v); complex<T> r = 0; const auto * ptr = hlp.p0r; @@ -806,6 +819,7 @@ template<typename T, typename Serv> void grid2x_c r += tmp*ku[cu]; ptr += jump; } + if (flip) r=conj(r); srv.setVis(ipart, r*emb); } } @@ -1067,7 +1081,7 @@ template<typename T, typename Serv> void wstack_common( #pragma omp parallel for num_threads(nthreads) reduction(min:wmin) reduction(max:wmax) for (size_t ipart=0; ipart<nvis; ++ipart) { - wval[ipart] = srv.getCoord(ipart).w; + wval[ipart] = abs(srv.getCoord(ipart).w); wmin = min(wmin,wval[ipart]); wmax = max(wmax,wval[ipart]); } @@ -1283,8 +1297,8 @@ template<typename T> size_t getIdxSize(const Baselines<T> &baselines, for (int ichan=chbegin; ichan<chend; ++ichan) if (!flags(irow,ichan)) { - auto uvw = baselines.effectiveCoord(RowChan{irow,size_t(ichan)}); - if ((uvw.w>=wmin) && (uvw.w<wmax)) + auto w = abs(baselines.effectiveCoord(RowChan{irow,size_t(ichan)}).w); + if ((w>=wmin) && (w<wmax)) ++res; } return res; @@ -1313,6 +1327,7 @@ template<typename T> void fillIdx(const Baselines<T> &baselines, if (!flags(irow, ichan)) { auto uvw = baselines.effectiveCoord(RowChan{irow,size_t(ichan)}); + if (uvw.w<0) uvw.Flip(); if ((uvw.w>=wmin) && (uvw.w<wmax)) { T u, v; @@ -1332,8 +1347,8 @@ template<typename T> void fillIdx(const Baselines<T> &baselines, for (int ichan=chbegin; ichan<chend; ++ichan) if (!flags(irow, ichan)) { - auto uvw = baselines.effectiveCoord(RowChan{irow,size_t(ichan)}); - if ((uvw.w>=wmin) && (uvw.w<wmax)) + auto w = abs(baselines.effectiveCoord(RowChan{irow,size_t(ichan)}).w); + if ((w>=wmin) && (w<wmax)) res[acc[tmp[idx++]]++] = irow*nchan+ichan; } } @@ -1357,6 +1372,7 @@ template<typename T> vector<uint32_t> getWgtIndices(const Baselines<T> &baseline if ((!have_wgt) || (wgt(irow,ichan)!=0)) { auto uvw = baselines.effectiveCoord(RowChan{irow,size_t(ichan)}); + if (uvw.w<0) uvw.Flip(); T u, v; int iu0, iv0; gconf.getpix(uvw.u, uvw.v, u, v, iu0, iv0);