deal more gracefully with corner cases

gridder_cxx.h

@@ -399,7 +399,9 @@ template<typename T> class GridderConfig
     complex<T> wscreen(double x, double y, double w, bool adjoint) const
       {
       constexpr double pi = 3.141592653589793238462643383279502884197;
-      double nm1 = (-x-y)/(sqrt(1.-x-y)+1); // more accurate form of sqrt(1-x-y)
+      double tmp = 1-x-y;
+      if (tmp<0) return 0.;
+      double nm1 = (-x-y)/(sqrt(tmp)+1); // more accurate form of sqrt(1-x-y)
       double n = nm1+1., xn = 1./n;
       double phase = 2*pi*w*nm1;
       if (adjoint) phase *= -1;
@@ -967,8 +969,13 @@ template<typename T, typename T2> void apply_wcorr(const GridderConfig<T> &gconf
       {
       double fy = y0+j*psy;
       fy*=fy;
-      auto nm1 = (-fx-fy)/(sqrt(1.-fx-fy)+1.); // accurate form of sqrt(1-x-y)
-      T fct = kernel.corfac(nm1*dw);
+      T fct = 0;
+      double tmp = 1.-fx-fy;
+      if (tmp>=0)
+        {
+        auto nm1 = (-fx-fy)/(sqrt(1.-fx-fy)+1.); // accurate form of sqrt(1-x-y)
+        fct = (nm1<=-1) ? 0. : kernel.corfac(nm1*dw);
+        }
       size_t i2 = nx_dirty-i, j2 = ny_dirty-j;
       dirty(i,j)*=fct;
       if ((i>0)&&(i<i2))
@@ -1068,7 +1075,7 @@ template<typename T, typename Serv> void x2dirty_wstack(const Baselines<T> &base
   for (size_t iw=0; iw<nplanes; ++iw)
     {
     if (nvis_plane[iw]==0) continue;
-    double wcur = wmin+iw*dw;
+    T wcur = wmin+iw*dw;
     size_t cnt=0;
     vis_loc_.resize(nvis_plane[iw]);
     auto vis_loc = mav<complex<T>,1>(vis_loc_.data(), {nvis_plane[iw]});
@@ -1078,13 +1085,13 @@ template<typename T, typename Serv> void x2dirty_wstack(const Baselines<T> &base
     for (size_t ipart=0; ipart<nvis; ++ipart)
       if ((int(iw)>=minplane[ipart]) && (iw<minplane[ipart]+w_supp))
         mapper[cnt++] = ipart;
-myassert(cnt==nvis_plane[iw],"must not happen 2");
 
+    T tmpval=T(2)/(w_supp*dw);
 #pragma omp parallel for num_threads(nthreads)
     for (size_t i=0; i<nvis_plane[iw]; ++i)
       {
       auto ipart = mapper[i];
-      double x=min(1.,2./(w_supp*dw)*abs(wcur-wval[ipart]));
+      T x=min(T(1), tmpval*abs(wcur-wval[ipart]));
       vis_loc[i] = srv.getVis(ipart)*kernel(x);
       idx_loc[i] = srv.getIdx(ipart);
       }
@@ -1148,7 +1155,7 @@ template<typename T, typename Serv> void dirty2x_wstack(const Baselines<T> &base
   for (size_t iw=0; iw<nplanes; ++iw)
     {
     if (nvis_plane[iw]==0) continue;
-    double wcur = wmin+iw*dw;
+    T wcur = wmin+iw*dw;
     size_t cnt=0;
     vis_loc_.resize(nvis_plane[iw]);
     auto vis_loc = mav<complex<T>,1>(vis_loc_.data(), {nvis_plane[iw]});
@@ -1158,8 +1165,6 @@ template<typename T, typename Serv> void dirty2x_wstack(const Baselines<T> &base
     for (size_t ipart=0; ipart<nvis; ++ipart)
       if ((int(iw)>=minplane[ipart]) && (iw<minplane[ipart]+w_supp))
         mapper[cnt++] = ipart;
-myassert(cnt==nvis_plane[iw],"must not happen 2");
-
 #pragma omp parallel for num_threads(nthreads)
     for (size_t i=0; i<nvis_plane[iw]; ++i)
         idx_loc[i] = srv.getIdx(mapper[i]);
@@ -1168,11 +1173,12 @@ myassert(cnt==nvis_plane[iw],"must not happen 2");
     gconf.dirty2grid_c(tdirty2, grid);
     grid2vis_c(baselines, gconf, const_mav<uint32_t,1>(idx_loc), const_mav<complex<T>,2>(grid), vis_loc, const_mav<T,1>(nullptr,{0}));
 
+    T tmpval=T(2)/(w_supp*dw);
 #pragma omp parallel for num_threads(nthreads)
     for (size_t i=0; i<nvis_plane[iw]; ++i)
       {
       auto ipart = mapper[i];
-      double x=min(1.,2./(w_supp*dw)*abs(wcur-wval[ipart]));
+      T x=min(T(1), tmpval*abs(wcur-wval[ipart]));
       srv.addVis(ipart, vis_loc[i]*kernel(x));
       }
     }