factor out min/max computation for w; tweak beta values

gridder_cxx.h

@@ -375,7 +375,7 @@ class ES_Kernel
 
   public:
     ES_Kernel(size_t supp_, size_t nthreads)
-      : beta(2.3*supp_), p(int(1.5*supp_+2)), supp(supp_)
+      : beta(get_beta(supp_)*supp_), p(int(1.5*supp_+2)), supp(supp_)
       {
       legendre_prep(2*p,x,wgt,nthreads);
       psi=x;
@@ -392,16 +392,24 @@ class ES_Kernel
         tmp += wgt[i]*psi[i]*cos(pi*supp*v*x[i]);
       return 1./(supp*tmp);
       }
+    static double get_beta(size_t supp)
+      {
+      static const vector<double> opt_beta {-1, 0.14, 1.70, 2.08, 2.205, 2.26,
+        2.29, 2.307, 2.316, 2.3265, 2.3324, 2.282, 2.294, 2.304, 2.3138, 2.317};
+      myassert(supp<opt_beta.size(), "bad support size");
+      return opt_beta[supp];
+      }
+
     static size_t get_supp(double epsilon)
       {
-      static const vector<double> maxmaperr { 1e8, 0.32, 0.021, 6.2e-4,
-        1.08e-5, 1.25e-7, 8.25e-10, 5.70e-12, 1.22e-13, 2.48e-15, 4.82e-17,
-        6.74e-19, 5.41e-21, 4.41e-23, 7.88e-25, 3.9e-26 };
+      static const vector<double> maxmaperr { 1e8, 0.19, 2.98e-3, 5.98e-5,
+        1.11e-6, 2.01e-8, 3.55e-10, 5.31e-12, 8.81e-14, 1.34e-15, 2.17e-17,
+        2.12e-19, 2.88e-21, 3.92e-23, 8.21e-25, 7.13e-27 };
       double epssq = epsilon*epsilon;
 
       for (size_t i=1; i<maxmaperr.size(); ++i)
         if (epssq>maxmaperr[i]) return i;
-      myfail("requested epsilon too small - minimum is 2e-13");
+      myfail("requested epsilon too small - minimum is 1e-13");
       }
   };
 
@@ -552,7 +560,7 @@ class GridderConfig
         psx(pixsize_x), psy(pixsize_y),
         supp(ES_Kernel::get_supp(epsilon)), nsafe((supp+1)/2),
         nu(max(2*nsafe,2*nx_dirty)), nv(max(2*nsafe,2*ny_dirty)),
-        beta(2.3*supp),
+        beta(ES_Kernel::get_beta(supp)*supp),
         cfu(nx_dirty), cfv(ny_dirty), nthreads(nthreads_),
         ushift(supp*(-0.5)+1+nu), vshift(supp*(-0.5)+1+nv),
         maxiu0((nu+nsafe)-supp), maxiv0((nv+nsafe)-supp)
@@ -1214,16 +1222,14 @@ template<typename T> void apply_wcorr(const GridderConfig &gconf,
     }
   }
 
-template<typename Serv> void wstack_common(
-  const GridderConfig &gconf, const Serv &srv, double &wmin,
-  double &dw, size_t &nplanes, vector<size_t> &nvis_plane,
-  vector<int> &minplane, size_t verbosity)
+template<typename Serv> void wminmax(const GridderConfig &gconf,
+  const Serv &srv, double &wmin, double &wmax)
   {
   size_t nvis = srv.Nvis();
   size_t nthreads = gconf.Nthreads();
 
-  wmin=1e38;
-  double wmax=-1e38;
+  wmin= 1e38;
+  wmax=-1e38;
 // FIXME maybe this can be done more intelligently
 #pragma omp parallel for num_threads(nthreads) reduction(min:wmin) reduction(max:wmax)
   for (size_t ipart=0; ipart<nvis; ++ipart)
@@ -1232,6 +1238,18 @@ template<typename Serv> void wstack_common(
     wmin = min(wmin,wval);
     wmax = max(wmax,wval);
     }
+  }
+
+template<typename Serv> void wstack_common(
+  const GridderConfig &gconf, const Serv &srv, double &wmin,
+  double &dw, size_t &nplanes, vector<size_t> &nvis_plane,
+  vector<int> &minplane, size_t verbosity)
+  {
+  size_t nvis = srv.Nvis();
+  size_t nthreads = gconf.Nthreads();
+  double wmax;
+
+  wminmax(gconf, srv, wmin, wmax);
   if (verbosity>0) cout << "Using " << nthreads << " threads" << endl;
   if (verbosity>0) cout << "W range: " << wmin << " to " << wmax << endl;
   double x0 = -0.5*gconf.Nxdirty()*gconf.Pixsize_x(),