flip negative w values automatically

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);