Changes in trunk/src/STGrid.cpp [2371:2461]

File:

• trunk/src/STGrid.cpp (modified) (29 diffs)

trunk/src/STGrid.cpp

@@ -1 @@
-#include <iostream>
-#include <fstream>
-
+//
+// C++ Implementation: STGrid
+//
+// Description:
+//
+//
+// Author: Takeshi Nakazato <takeshi.nakazato@nao.ac.jp>, (C) 2011
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
 #include <casa/BasicSL/String.h>
 #include <casa/Arrays/Vector.h>
-#include <casa/Arrays/Matrix.h>
-#include <casa/Arrays/Cube.h>
 #include <casa/Arrays/ArrayMath.h>
-#include <casa/Arrays/ArrayPartMath.h>
 #include <casa/Quanta/Quantum.h>
 #include <casa/Quanta/QuantumHolder.h>
@@ -15 +20 @@
 #include <tables/Tables/Table.h>
 #include <tables/Tables/TableRecord.h>
+#include <tables/Tables/TableRow.h>
 #include <tables/Tables/ExprNode.h>
 #include <tables/Tables/ScalarColumn.h>
 #include <tables/Tables/ArrayColumn.h>
+#include <tables/Tables/TableCopy.h>
 
 #include <measures/Measures/MDirection.h>
 
-#include <MathUtils.h>
+#include "MathUtils.h"
+#include <atnf/PKSIO/SrcType.h>
 
 #include "STGrid.h"
 
 using namespace std ;
+using namespace concurrent ;
 using namespace casa ;
 using namespace asap ;
@@ -31 +40 @@
 namespace asap {
 
+// for performance check
+double eToInt = 0.0 ;
+double eGetWeight = 0.0 ;
+
 // constructor
 STGrid::STGrid()
+  : vshape_( 1 ), wshape_( 2 ), dshape_( 2 )
 {
   init() ;
@@ -38 +52 @@
 
 STGrid::STGrid( const string infile )
+  : vshape_( 1 ), wshape_( 2 ), dshape_( 2 )
 {
   init() ;
 
   setFileIn( infile ) ;
+}
+
+STGrid::STGrid( const vector<string> infile )
+{
+  init() ;
+
+  setFileList( infile ) ;
 }
 
@@ -60 +82 @@
   userSupport_ = -1 ;
   convSampling_ = 100 ;
+  nprocessed_ = 0 ;
+  nchunk_ = 0 ;
+
+  // initialize user input 
+  nxUI_ = -1 ;
+  nyUI_ = -1 ;
+  cellxUI_ = "" ;
+  cellyUI_ = "" ;
+  centerUI_ = "" ;
+  doclip_ = False ;
 }
 
 void STGrid::setFileIn( const string infile )
 {
+  nfile_ = 1 ;
   String name( infile ) ;
-  if ( infile_.compare( name ) != 0 ) {
-    infile_ = String( infile ) ;
-    tab_ = Table( infile_ ) ;
+  infileList_.resize( nfile_ ) ;
+  infileList_[0] = String(infile) ;
+}
+
+void STGrid::setFileList( const vector<string> infile )
+{
+  nfile_ = infile.size() ;
+  infileList_.resize( nfile_ ) ;
+  for ( uInt i = 0 ; i < nfile_ ; i++ ) {
+    infileList_[i] = infile[i] ;
   }
 }
@@ -74 +114 @@
 {
   pollist_.assign( Vector<uInt>( pols ) ) ;
-  cout << "pollist_ = " << pollist_ << endl ;
 }
 
@@ -80 +119 @@
 {
   scanlist_.assign( Vector<uInt>( scans ) ) ;
-  cout << "scanlist_ = " << scanlist_ << endl ;
 }
 
@@ -87 +125 @@
   wtype_ = String( wType ) ;
   wtype_.upcase() ;
-  cout << "wtype_ = " << wtype_ << endl ; 
 }
 
@@ -96 +133 @@
                           string scenter ) 
 {
-  ROArrayColumn<Double> dirCol( tab_, "DIRECTION" ) ;
-  Matrix<Double> direction = dirCol.getColumn() ;
-  Double rmax, rmin, dmax, dmin ;
-  minMax( rmin, rmax, direction.row( 0 ) ) ;
-  minMax( dmin, dmax, direction.row( 1 ) ) ;
-
-  Int npx = (Int)nx ;
-  Int npy = (Int)ny ;
-  String cellx( scellx ) ;
-  String celly( scelly ) ;
-  String center( scenter ) ;
-  setupGrid( npx, npy, 
-             cellx, celly, 
-             rmin, rmax, 
-             dmin, dmax, 
-             center ) ;
+  nxUI_ = (Int)nx ;
+  nyUI_ = (Int)ny ;
+  cellxUI_ = String( scellx ) ;
+  cellyUI_ = String( scelly ) ;
+  centerUI_ = String( scenter ) ;
 }
   
@@ -150 +176 @@
                 Double*);
 }
-void STGrid::grid() 
+void STGrid::call_ggridsd( Array<Double> &xypos,
+                           Array<Complex> &spectra,
+                           Int &nvispol,
+                           Int &nvischan,
+                           Array<Int> &flagtra,
+                           Array<Int> &flagrow,
+                           Array<Float> &weight,
+                           Int &nrow,
+                           Int &irow,
+                           Array<Complex> &gdata,
+                           Array<Float> &gwgt,
+                           Int &nx,
+                           Int &ny,
+                           Int &npol,
+                           Int &nchan,
+                           Int &support,
+                           Int &sampling,
+                           Vector<Float> &convFunc,
+                           Int *chanMap,
+                           Int *polMap ) 
+{
+  // parameters for gridding
+  Int idopsf = 0 ;
+  Int len = npol*nchan ;
+  Double *sumw_p = new Double[len] ;
+  {
+    Double *work_p = sumw_p ;
+    for ( Int i = 0 ; i < len ; i++ ) {
+      *work_p = 0.0 ;
+      work_p++ ;
+    }
+  }
+
+  // prepare pointer
+  Bool deletePos, deleteData, deleteWgt, deleteFlag, deleteFlagR, deleteConv, deleteDataG, deleteWgtG ;
+  Double *xy_p = xypos.getStorage( deletePos ) ;
+  const Complex *values_p = spectra.getStorage( deleteData ) ;
+  const Int *flag_p = flagtra.getStorage( deleteFlag ) ;
+  const Int *rflag_p = flagrow.getStorage( deleteFlagR ) ;
+  const Float *wgt_p = weight.getStorage( deleteWgt ) ;
+  Complex *grid_p = gdata.getStorage( deleteDataG ) ;
+  Float *wgrid_p = gwgt.getStorage( deleteWgtG ) ;
+  Float *conv_p = convFunc.getStorage( deleteConv ) ;
+
+  // pass copy of irow to ggridsd since it will be modified in theroutine
+  Int irowCopy = irow ;
+      
+  // call ggridsd
+  ggridsd( xy_p,
+           values_p,
+           &nvispol,
+           &nvischan,
+           &idopsf,
+           flag_p,
+           rflag_p,
+           wgt_p,
+           &nrow,
+           &irowCopy,
+           grid_p,
+           wgrid_p,
+           &nx,
+           &ny,
+           &npol,
+           &nchan,
+           &support,
+           &sampling,
+           conv_p,
+           chanMap,
+           polMap,
+           sumw_p ) ;
+
+  // finalization
+  xypos.putStorage( xy_p, deletePos ) ;
+  spectra.freeStorage( values_p, deleteData ) ;
+  flagtra.freeStorage( flag_p, deleteFlag ) ;
+  flagrow.freeStorage( rflag_p, deleteFlagR ) ;
+  weight.freeStorage( wgt_p, deleteWgt ) ;
+  gdata.putStorage( grid_p, deleteDataG ) ;
+  gwgt.putStorage( wgrid_p, deleteWgtG ) ;
+  convFunc.putStorage( conv_p, deleteConv ) ;
+  delete sumw_p ;
+}
+
+#define NEED_UNDERSCORES
+#if defined(NEED_UNDERSCORES)
+#define ggridsd2 ggridsd2_
+#endif
+extern "C" { 
+   void ggridsd2(Double*,
+                 const Complex*,
+                 Int*,
+                 Int*,
+                 Int*,
+                 const Int*,
+                 const Int*,
+                 const Float*,
+                 Int*,
+                 Int*,
+                 Complex*,
+                 Float*,
+                 Int*,
+                 Complex*,
+                 Float*,
+                 Float*,
+                 Complex*,
+                 Float*,
+                 Float*,
+                 Int*,
+                 Int*,
+                 Int *,
+                 Int *,
+                 Int*,
+                 Int*,
+                 Float*,
+                 Int*,
+                 Int*,
+                 Double*);
+}
+void STGrid::call_ggridsd2( Array<Double> &xypos,
+                            Array<Complex> &spectra,
+                            Int &nvispol,
+                            Int &nvischan,
+                            Array<Int> &flagtra,
+                            Array<Int> &flagrow,
+                            Array<Float> &weight,
+                            Int &nrow,
+                            Int &irow,
+                            Array<Complex> &gdata,
+                            Array<Float> &gwgt,
+                            Array<Int> &npoints,
+                            Array<Complex> &clipmin,
+                            Array<Float> &clipwmin,
+                            Array<Float> &clipcmin,
+                            Array<Complex> &clipmax,
+                            Array<Float> &clipwmax,
+                            Array<Float> &clipcmax,
+                            Int &nx,
+                            Int &ny,
+                            Int &npol,
+                            Int &nchan,
+                            Int &support,
+                            Int &sampling,
+                            Vector<Float> &convFunc,
+                            Int *chanMap,
+                            Int *polMap ) 
+{
+  // parameters for gridding
+  Int idopsf = 0 ;
+  Int len = npol*nchan ;
+  Double *sumw_p = new Double[len] ;
+  {
+    Double *work_p = sumw_p ;
+    for ( Int i = 0 ; i < len ; i++ ) {
+      *work_p = 0.0 ;
+      work_p++ ;
+    }
+  }
+
+  // prepare pointer
+  Bool deletePos, deleteData, deleteWgt, deleteFlag, deleteFlagR, deleteConv, deleteDataG, deleteWgtG, deleteNpts, deleteCMin, deleteCWMin, deleteCCMin, deleteCMax, deleteCWMax, deleteCCMax ;
+  Double *xy_p = xypos.getStorage( deletePos ) ;
+  const Complex *values_p = spectra.getStorage( deleteData ) ;
+  const Int *flag_p = flagtra.getStorage( deleteFlag ) ;
+  const Int *rflag_p = flagrow.getStorage( deleteFlagR ) ;
+  const Float *wgt_p = weight.getStorage( deleteWgt ) ;
+  Complex *grid_p = gdata.getStorage( deleteDataG ) ;
+  Float *wgrid_p = gwgt.getStorage( deleteWgtG ) ;
+  Float *conv_p = convFunc.getStorage( deleteConv ) ;
+  Int *npts_p = npoints.getStorage( deleteNpts ) ;
+  Complex *cmin_p = clipmin.getStorage( deleteCMin ) ;
+  Float *cwmin_p = clipwmin.getStorage( deleteCWMin ) ;
+  Float *ccmin_p = clipcmin.getStorage( deleteCCMin ) ;
+  Complex *cmax_p = clipmax.getStorage( deleteCMax ) ;
+  Float *cwmax_p = clipwmax.getStorage( deleteCWMax ) ;
+  Float *ccmax_p = clipcmax.getStorage( deleteCCMax ) ;
+
+  // pass copy of irow to ggridsd since it will be modified in theroutine
+  Int irowCopy = irow ;
+      
+  // call ggridsd
+  ggridsd2( xy_p,
+            values_p,
+            &nvispol,
+            &nvischan,
+            &idopsf,
+            flag_p,
+            rflag_p,
+            wgt_p,
+            &nrow,
+            &irowCopy,
+            grid_p,
+            wgrid_p,
+            npts_p,
+            cmin_p,
+            cwmin_p,
+            ccmin_p,
+            cmax_p,
+            cwmax_p,
+            ccmax_p,
+            &nx,
+            &ny,
+            &npol,
+            &nchan,
+            &support,
+            &sampling,
+            conv_p,
+            chanMap,
+            polMap,
+            sumw_p ) ;
+
+  // finalization
+  xypos.putStorage( xy_p, deletePos ) ;
+  spectra.freeStorage( values_p, deleteData ) ;
+  flagtra.freeStorage( flag_p, deleteFlag ) ;
+  flagrow.freeStorage( rflag_p, deleteFlagR ) ;
+  weight.freeStorage( wgt_p, deleteWgt ) ;
+  gdata.putStorage( grid_p, deleteDataG ) ;
+  gwgt.putStorage( wgrid_p, deleteWgtG ) ;
+  convFunc.putStorage( conv_p, deleteConv ) ;
+  clipmin.putStorage( cmin_p, deleteCMin ) ;
+  clipwmin.putStorage( cwmin_p, deleteCWMin ) ;
+  clipcmin.putStorage( ccmin_p, deleteCCMin ) ;
+  clipmax.putStorage( cmax_p, deleteCMax ) ;
+  clipwmax.putStorage( cwmax_p, deleteCWMax ) ;
+  clipcmax.putStorage( ccmax_p, deleteCCMax ) ;
+  delete sumw_p ;
+}
+
+void STGrid::grid()
 {
   LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
-
-  // retrieve data
-  Cube<Float> spectra ;
-  Matrix<Double> direction ;
-  Cube<uChar> flagtra ;
-  Matrix<uInt> rflag ;
-  Matrix<Float> weight ;
-  double t0, t1 ;
+  double t0,t1 ;
+
+  // data selection
   t0 = mathutil::gettimeofday_sec() ;
-  getData( spectra, direction, flagtra, rflag, weight ) ;
+  selectData() ;
   t1 = mathutil::gettimeofday_sec() ;
-  os << "getData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ; 
-  IPosition sshape = spectra.shape() ;
-  //os << "spectra.shape()=" << spectra.shape() << LogIO::POST ;
-  //os << "max(spectra) = " << max(spectra) << LogIO::POST ;
-  //os << "weight = " << weight << LogIO::POST ;
-
-  // flagtra: uChar -> Int
-  // rflag: uInt -> Int
-  Cube<Int> flagI ;
-  Matrix<Int> rflagI ;
-  t0 = mathutil::gettimeofday_sec() ;
-  toInt( &flagtra, &flagI ) ;
-  toInt( &rflag, &rflagI ) ;
-  t1 = mathutil::gettimeofday_sec() ;
-  os << "toInt: elapsed time is " << t1-t0 << " sec." << LogIO::POST ; 
-  
+  os << LogIO::DEBUGGING << "selectData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
+
+  setupGrid() ;
+  setupArray() ;
+
+  if ( wtype_.compare("UNIFORM") != 0 &&
+       wtype_.compare("TINT") != 0 && 
+       wtype_.compare("TSYS") != 0 &&
+       wtype_.compare("TINTSYS") != 0 ) {
+    LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
+    os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
+    wtype_ = "UNIFORM" ;
+  }
+
   // grid parameter
   os << LogIO::DEBUGGING ;
+  os << "----------" << endl ;
+  os << "Data selection summary" << endl ;
+  os << "   ifno = " << ifno_ << endl ;
+  os << "   pollist = " << pollist_ << endl ;
+  os << "   scanlist = " << scanlist_ << endl ;
   os << "----------" << endl ;
   os << "Grid parameter summary" << endl ;
@@ -187 +439 @@
   os << "   (cellx,celly) = (" << cellx_ << "," << celly_ << ")" << endl ;
   os << "   center = " << center_ << endl ;
+  os << "   weighting = " << wtype_ << endl ;
+  os << "   convfunc = " << convType_ << " with support " << convSupport_ << endl ; 
+  os << "   doclip = " << (doclip_?"True":"False") << endl ;
   os << "----------" << LogIO::POST ;
   os << LogIO::NORMAL ;
 
+  if ( doclip_ )
+    gridPerRowWithClipping() ;
+  else 
+    gridPerRow() ;
+}
+
+void STGrid::updateChunkShape()
+{
+  // TODO: nchunk_ must be determined from nchan_, npol_, and (nx_,ny_) 
+  //       by considering data size to be allocated for ggridsd input/output
+  nchunk_ = 400 ;
+  nchunk_ = min( nchunk_, nrow_ ) ;
+  vshape_ = IPosition( 1, nchunk_ ) ;
+  wshape_ = IPosition( 2, nchan_, nchunk_ ) ;
+  dshape_ = IPosition( 2, 2, nchunk_ ) ;
+}
+
+struct STGChunk {
+  Int nrow ;
+  Array<Complex> spectra;
+  Array<Int> flagtra;
+  Array<Int> rflag;
+  Array<Float> weight;
+  Array<Double> direction;
+  STGChunk(IPosition const &wshape, IPosition const &vshape,
+           IPosition const &dshape)
+    : spectra(wshape), flagtra(wshape), rflag(vshape), weight(wshape),
+      direction(dshape)
+  { }
+};
+
+struct STCommonData {
+  Int gnx;
+  Int gny;
+  Int *chanMap;
+  Vector<Float> convFunc ;
+  Array<Complex> gdataArrC;
+  Array<Float> gwgtArr;
+  STCommonData(IPosition const &gshape, Array<Float> const &data)
+    : gdataArrC(gshape, 0.0), gwgtArr(data) {}
+};
+
+struct STCommonDataWithClipping {
+  Int gnx;
+  Int gny;
+  Int *chanMap;
+  Vector<Float> convFunc ;
+  Array<Complex> gdataArrC;
+  Array<Float> gwgtArr;
+  Array<Int> npoints ;
+  Array<Complex> clipMin ;
+  Array<Float> clipWMin ;
+  Array<Float> clipCMin ;
+  Array<Complex> clipMax ;
+  Array<Float> clipWMax ;
+  Array<Float> clipCMax ;  
+  STCommonDataWithClipping(IPosition const &gshape, 
+                           IPosition const &pshape, 
+                           Array<Float> const &data)
+    : gdataArrC(gshape, 0.0), 
+      gwgtArr(data), 
+      npoints(pshape, 0),
+      clipMin(gshape, Complex(FLT_MAX,0.0)),
+      clipWMin(gshape, 0.0),
+      clipCMin(gshape, 0.0),
+      clipMax(gshape, Complex(-FLT_MAX,0.0)),
+      clipWMax(gshape, 0.0),
+      clipCMax(gshape, 0.0)
+  {}
+};
+
+#define DO_AHEAD 3
+
+struct STContext {
+  STCommonData &common;
+  FIFO<STGChunk *, DO_AHEAD> queue;
+  STGrid *const self;
+  const Int pol;
+  STContext(STGrid *obj, STCommonData &common, Int pol)
+    : self(obj), common(common), pol(pol) {}
+};
+
+struct STContextWithClipping {
+  STCommonDataWithClipping &common;
+  FIFO<STGChunk *, DO_AHEAD> queue;
+  STGrid *const self;
+  const Int pol;
+  STContextWithClipping(STGrid *obj, STCommonDataWithClipping &common, Int pol)
+    : self(obj), common(common), pol(pol) {}
+};
+
+
+bool STGrid::produceChunk(void *ctx) throw(PCException)
+{
+  STContext &context = *(STContext *)ctx;
+  if ( context.self->nprocessed_ >= context.self->nrow_ ) {
+    return false;
+  }
+  STGChunk *chunk = new STGChunk(context.self->wshape_,
+                                 context.self->vshape_,
+                                 context.self->dshape_);
+
+  double t0 = mathutil::gettimeofday_sec() ;
+  chunk->nrow = context.self->getDataChunk(
+        context.self->wshape_, context.self->vshape_, context.self->dshape_,
+        chunk->spectra, chunk->direction,
+        chunk->flagtra, chunk->rflag, chunk->weight);
+  double t1 = mathutil::gettimeofday_sec() ;
+  context.self->eGetData_ += t1-t0 ;
+
+  context.queue.lock();
+  context.queue.put(chunk);
+  context.queue.unlock();
+  return true;
+}
+
+void STGrid::consumeChunk(void *ctx) throw(PCException)
+{
+  STContext &context = *(STContext *)ctx;
+  STGChunk *chunk = NULL;
+  try {
+    context.queue.lock();
+    chunk = context.queue.get();
+    context.queue.unlock();
+  } catch (FullException &e) {
+    context.queue.unlock();
+    // TODO: log error
+    throw PCException();
+  }
+
+  double t0, t1 ;
+  // world -> pixel
+  Array<Double> xypos( context.self->dshape_ ) ;
+  t0 = mathutil::gettimeofday_sec() ;
+  context.self->toPixel( chunk->direction, xypos ) ;
+  t1 = mathutil::gettimeofday_sec() ;
+  context.self->eToPixel_ += t1-t0 ;
+   
+  // call ggridsd
+  Int nvispol = 1 ;
+  Int irow = -1 ;
+  t0 = mathutil::gettimeofday_sec() ;
+  context.self->call_ggridsd( xypos,
+                chunk->spectra,
+                nvispol,
+                context.self->nchan_,
+                chunk->flagtra,
+                chunk->rflag,
+                chunk->weight,
+                chunk->nrow,
+                irow,
+                context.common.gdataArrC,
+                context.common.gwgtArr,
+                context.common.gnx,
+                context.common.gny,
+                context.self->npol_,
+                context.self->nchan_,
+                context.self->convSupport_,
+                context.self->convSampling_,
+                context.common.convFunc,
+                context.common.chanMap,
+                (Int*)&context.pol ) ;
+  t1 = mathutil::gettimeofday_sec() ;
+  context.self->eGGridSD_ += t1-t0 ;
+  
+  delete chunk;
+}
+
+void STGrid::gridPerRow()
+{
+  LogIO os( LogOrigin("STGrid", "gridPerRow", WHERE) ) ;
+  double t0, t1 ;
+
+
+  // grid data
+  // Extend grid plane with convSupport_
+  //   Int gnx = nx_+convSupport_*2 ;
+  //   Int gny = ny_+convSupport_*2 ;
+  Int gnx = nx_;
+  Int gny = ny_;
+
+  IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
+  // 2011/12/20 TN
+  // data_ and gwgtArr share storage
+  data_.resize( gshape ) ;
+  data_ = 0.0 ;
+  STCommonData common = STCommonData(gshape, data_);
+  common.gnx = gnx ;
+  common.gny = gny ;
+
+  // parameters for gridding
+  Int *chanMap = new Int[nchan_] ;
+  for ( Int i = 0 ; i < nchan_ ; i++ ) {
+    chanMap[i] = i ;
+  }
+  common.chanMap = chanMap;
+
   // convolution kernel
-  Vector<Float> convFunc ;
   t0 = mathutil::gettimeofday_sec() ;
-  setConvFunc( convFunc ) ;
+  setConvFunc( common.convFunc ) ;
   t1 = mathutil::gettimeofday_sec() ;
-  os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ; 
-  //cout << "convSupport=" << convSupport_ << endl ;
-  //cout << "convFunc=" << convFunc << endl ;
-
+  os << LogIO::DEBUGGING << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ; 
+
+  // for performance check
+  eGetData_ = 0.0 ;
+  eToPixel_ = 0.0 ;
+  eGGridSD_ = 0.0 ;
+  double eInitPol = 0.0 ;
+
+  for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) {
+    initTable( ifile ) ;
+
+    os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ;   
+    Broker broker = Broker(produceChunk, consumeChunk);
+    for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
+      t0 = mathutil::gettimeofday_sec() ;
+      initPol( ipol ) ; // set ptab_ and attach()
+      t1 = mathutil::gettimeofday_sec() ;
+      eInitPol += t1-t0 ;
+      
+      STContext context(this, common, ipol);
+      
+      os << "start pol " << ipol << LogIO::POST ;
+      
+      nprocessed_ = 0 ;
+#if 1
+      broker.runProducerAsMasterThread(&context, DO_AHEAD);
+#else
+      for (;;) {
+        bool produced = produceChunk(&context);
+        if (! produced) {
+          break;
+        }
+        consumeChunk(&context);
+      }
+#endif
+
+      os << "end pol " << ipol << LogIO::POST ;
+
+    }
+    os << "end table " << ifile << LogIO::POST ;   
+  }
+  os << LogIO::DEBUGGING << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ; 
+  os << LogIO::DEBUGGING << "getData: elapsed time is " << eGetData_-eToInt-eGetWeight << " sec." << LogIO::POST ; 
+  os << LogIO::DEBUGGING << "toPixel: elapsed time is " << eToPixel_ << " sec." << LogIO::POST ; 
+  os << LogIO::DEBUGGING << "ggridsd: elapsed time is " << eGGridSD_ << " sec." << LogIO::POST ; 
+  os << LogIO::DEBUGGING << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
+  os << LogIO::DEBUGGING << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
+  
+  delete chanMap ;
+
+  // set data
+  setData( common.gdataArrC, common.gwgtArr ) ;
+
+}
+
+void STGrid::consumeChunkWithClipping(void *ctx) throw(PCException)
+{
+  STContextWithClipping &context = *(STContextWithClipping *)ctx;
+  STGChunk *chunk = NULL;
+  try {
+    context.queue.lock();
+    chunk = context.queue.get();
+    context.queue.unlock();
+  } catch (FullException &e) {
+    context.queue.unlock();
+    // TODO: log error
+    throw PCException();
+  }
+
+  double t0, t1 ;
   // world -> pixel
-  Matrix<Double> xypos( direction.shape(), 0.0 ) ;
+  Array<Double> xypos( context.self->dshape_ ) ;
   t0 = mathutil::gettimeofday_sec() ;
-  toPixel( direction, xypos ) ;  
+  context.self->toPixel( chunk->direction, xypos ) ;
   t1 = mathutil::gettimeofday_sec() ;
-  os << "toPixel: elapsed time is " << t1-t0 << " sec." << LogIO::POST ; 
+  context.self->eToPixel_ += t1-t0 ;
+   
+  // call ggridsd
+  Int nvispol = 1 ;
+  Int irow = -1 ;
+  t0 = mathutil::gettimeofday_sec() ;
+  context.self->call_ggridsd2( xypos,
+                chunk->spectra,
+                nvispol,
+                context.self->nchan_,
+                chunk->flagtra,
+                chunk->rflag,
+                chunk->weight,
+                chunk->nrow,
+                irow,
+                context.common.gdataArrC,
+                context.common.gwgtArr,
+                context.common.npoints,
+                context.common.clipMin,
+                context.common.clipWMin,
+                context.common.clipCMin,
+                context.common.clipMax,
+                context.common.clipWMax,
+                context.common.clipCMax,
+                context.common.gnx,
+                context.common.gny,
+                context.self->npol_,
+                context.self->nchan_,
+                context.self->convSupport_,
+                context.self->convSampling_,
+                context.common.convFunc,
+                context.common.chanMap,
+                (Int*)&context.pol ) ;
+  t1 = mathutil::gettimeofday_sec() ;
+  context.self->eGGridSD_ += t1-t0 ;
   
-  // call ggridsd
-  Bool deletePos, deleteData, deleteWgt, deleteFlag, deleteFlagR, deleteConv, deleteDataG, deleteWgtG ;
-  Double *xypos_p = xypos.getStorage( deletePos ) ;
-  Cube<Complex> dataC( spectra.shape(), 0.0 ) ;
-  setReal( dataC, spectra ) ;
-  const Complex *data_p = dataC.getStorage( deleteData ) ;
-  const Float *wgt_p = weight.getStorage( deleteWgt ) ;
-  const Int *flag_p = flagI.getStorage( deleteFlag ) ;
-  const Int *rflag_p = rflagI.getStorage( deleteFlagR ) ;
-  Float *conv_p = convFunc.getStorage( deleteConv ) ;
+  delete chunk;
+}
+
+void STGrid::gridPerRowWithClipping()
+{
+  LogIO os( LogOrigin("STGrid", "gridPerRowWithClipping", WHERE) ) ;
+  double t0, t1 ;
+
+
+  // grid data
   // Extend grid plane with convSupport_
-  Int gnx = nx_ ;
-  Int gny = ny_ ;
-//   Int gnx = nx_+convSupport_*2 ;
-//   Int gny = ny_+convSupport_*2 ;
+  //   Int gnx = nx_+convSupport_*2 ;
+  //   Int gny = ny_+convSupport_*2 ;
+  Int gnx = nx_;
+  Int gny = ny_;
+
   IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
-  Array<Complex> gdataArrC( gshape, 0.0 ) ;
-  Array<Float> gwgtArr( gshape, 0.0 ) ;
-  Complex *gdata_p = gdataArrC.getStorage( deleteDataG ) ;
-  Float *wdata_p = gwgtArr.getStorage( deleteWgtG ) ;
-  Int idopsf = 0 ;
+  IPosition pshape( 3, gnx, gny, npol_ ) ;
+  // 2011/12/20 TN
+  // data_ and gwgtArr share storage
+  data_.resize( gshape ) ;
+  data_ = 0.0 ;
+  STCommonDataWithClipping common = STCommonDataWithClipping( gshape,
+                                                              pshape, 
+                                                              data_ ) ;
+  common.gnx = gnx ;
+  common.gny = gny ;
+
+  // parameters for gridding
   Int *chanMap = new Int[nchan_] ;
-  {
-    Int *work_p = chanMap ;
-    for ( Int i = 0 ; i < nchan_ ; i++ ) {
-      *work_p = i ;
-      work_p++ ;
-    }
-  }
-  Int *polMap = new Int[npol_] ;
-  {
-    Int *work_p = polMap ;
-    for ( Int i = 0 ; i < npol_ ; i++ ) {
-      *work_p = i ;
-      work_p++ ;
-    }
-  }
-  Double *sumw_p = new Double[npol_*nchan_] ;
-  {
-    Double *work_p = sumw_p ;
-    for ( Int i = 0 ; i < npol_*nchan_ ; i++ ) {
-      *work_p = 0.0 ;
-      work_p++ ;
-    }
-  }
+  for ( Int i = 0 ; i < nchan_ ; i++ ) {
+    chanMap[i] = i ;
+  }
+  common.chanMap = chanMap;
+
+  // convolution kernel
   t0 = mathutil::gettimeofday_sec() ;
-  Int irow = -1 ;
-  ggridsd( xypos_p,
-           data_p,
-           &npol_,
-           &nchan_,
-           &idopsf,
-           flag_p,
-           rflag_p,
-           wgt_p,
-           &nrow_,
-           &irow,
-           gdata_p,
-           wdata_p, 
-           &gnx,
-           &gny,
-           &npol_,
-           &nchan_,
-           &convSupport_,
-           &convSampling_,
-           conv_p,
-           chanMap,
-           polMap,
-           sumw_p ) ;
+  setConvFunc( common.convFunc ) ;
   t1 = mathutil::gettimeofday_sec() ;
-  os << "ggridsd: elapsed time is " << t1-t0 << " sec." << LogIO::POST ; 
-  xypos.putStorage( xypos_p, deletePos ) ;
-  dataC.freeStorage( data_p, deleteData ) ;
-  weight.freeStorage( wgt_p, deleteWgt ) ;
-  flagI.freeStorage( flag_p, deleteFlag ) ;
-  rflagI.freeStorage( rflag_p, deleteFlagR ) ;
-  convFunc.putStorage( conv_p, deleteConv ) ;
-  delete polMap ;
+  os << LogIO::DEBUGGING << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ; 
+
+  // for performance check
+  eGetData_ = 0.0 ;
+  eToPixel_ = 0.0 ;
+  eGGridSD_ = 0.0 ;
+  double eInitPol = 0.0 ;
+
+  for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) {
+    initTable( ifile ) ;
+
+    os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ;   
+    Broker broker = Broker(produceChunk, consumeChunkWithClipping);
+    for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
+      t0 = mathutil::gettimeofday_sec() ;
+      initPol( ipol ) ; // set ptab_ and attach()
+      t1 = mathutil::gettimeofday_sec() ;
+      eInitPol += t1-t0 ;
+      
+      STContextWithClipping context(this, common, ipol);
+      
+      os << "start pol " << ipol << LogIO::POST ;
+      
+      nprocessed_ = 0 ;
+#if 1
+      broker.runProducerAsMasterThread(&context, DO_AHEAD);
+#else
+      for (;;) {
+        bool produced = produceChunk(&context);
+        if (! produced) {
+          break;
+        }
+        consumeChunkWithClipping(&context);
+      }
+#endif
+
+      os << "end pol " << ipol << LogIO::POST ;
+
+    }
+    os << "end table " << ifile << LogIO::POST ;   
+  }
+  os << LogIO::DEBUGGING << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ; 
+  os << LogIO::DEBUGGING << "getData: elapsed time is " << eGetData_-eToInt-eGetWeight << " sec." << LogIO::POST ; 
+  os << LogIO::DEBUGGING << "toPixel: elapsed time is " << eToPixel_ << " sec." << LogIO::POST ; 
+  os << LogIO::DEBUGGING << "ggridsd2: elapsed time is " << eGGridSD_ << " sec." << LogIO::POST ; 
+  os << LogIO::DEBUGGING << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
+  os << LogIO::DEBUGGING << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
+  
   delete chanMap ;
-  gdataArrC.putStorage( gdata_p, deleteDataG ) ;
-  gwgtArr.putStorage( wdata_p, deleteWgtG ) ;
-  Array<Float> gdataArr = real( gdataArrC ) ;
-  setData( data_, gdataArr, gwgtArr ) ;
-  //Matrix<Double> sumWeight( IPosition( 2, npol_, nchan_ ), sumw_p, TAKE_OVER ) ;
-  delete sumw_p ;
-  //cout << "sumWeight = " << sumWeight << endl ;
-//   os << "gdataArr = " << gdataArr << LogIO::POST ;
-//   os << "gwgtArr = " << gwgtArr << LogIO::POST ;
-//   os << "data_ " << data_ << LogIO::POST ;
-}
-
-void STGrid::setData( Array<Float> &data,
-                      Array<Float> &gdata,
+
+  // clip min and max in each grid
+//   os << "BEFORE CLIPPING" << LogIO::POST ;
+//   os << "gdataArrC=" << common.gdataArrC << LogIO::POST ;
+//   os << "gwgtArr=" << common.gwgtArr << LogIO::POST ;
+  t0 = mathutil::gettimeofday_sec() ;
+  clipMinMax( common.gdataArrC, 
+              common.gwgtArr,
+              common.npoints,
+              common.clipMin,
+              common.clipWMin,
+              common.clipCMin,
+              common.clipMax,
+              common.clipWMax,
+              common.clipCMax ) ;
+  t1 = mathutil::gettimeofday_sec() ;
+  os << LogIO::DEBUGGING << "clipMinMax: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
+//   os << "AFTER CLIPPING" << LogIO::POST ;
+//   os << "gdataArrC=" << common.gdataArrC << LogIO::POST ;
+//   os << "gwgtArr=" << common.gwgtArr << LogIO::POST ;
+
+  // set data
+  setData( common.gdataArrC, common.gwgtArr ) ;
+
+}
+
+void STGrid::clipMinMax( Array<Complex> &grid,
+                         Array<Float> &weight,
+                         Array<Int> &npoints,
+                         Array<Complex> &clipmin,
+                         Array<Float> &clipwmin,
+                         Array<Float> &clipcmin,
+                         Array<Complex> &clipmax,
+                         Array<Float> &clipwmax,
+                         Array<Float> &clipcmax )
+{
+  //LogIO os( LogOrigin("STGrid","clipMinMax",WHERE) ) ;
+
+  // prepare pointers
+  Bool delG, delW, delNP, delCMin, delCWMin, delCCMin, delCMax, delCWMax, delCCMax ;
+  Complex *grid_p = grid.getStorage( delG ) ;
+  Float *wgt_p = weight.getStorage( delW ) ;
+  const Int *npts_p = npoints.getStorage( delNP ) ;
+  const Complex *cmin_p = clipmin.getStorage( delCMin ) ;
+  const Float *cwmin_p = clipwmin.getStorage( delCWMin ) ;
+  const Float *ccmin_p = clipcmin.getStorage( delCCMin ) ;
+  const Complex *cmax_p = clipmax.getStorage( delCMax ) ;
+  const Float *cwmax_p = clipwmax.getStorage( delCWMax ) ;
+  const Float *ccmax_p = clipcmax.getStorage( delCCMax ) ;
+
+  const IPosition &gshape = grid.shape() ;
+  long offset = gshape[0] * gshape[1] * gshape[2] ; // nx * ny * npol
+  Int nchan = gshape[3] ;
+  long origin = nchan * offset ;
+  for ( long i = 0 ; i < offset ; i++ ) {
+    if ( *npts_p > 2 ) {
+      for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
+        // clip minimum and maximum
+        *grid_p -= (*cmin_p)*(*cwmin_p)*(*ccmin_p)
+          + (*cmax_p)*(*cwmax_p)*(*ccmax_p) ;
+        *wgt_p -= (*cwmin_p)*(*ccmin_p)
+          + (*cwmax_p)*(*ccmax_p) ;
+        
+        grid_p += offset ;
+        wgt_p += offset ;
+        cmin_p += offset ;
+        cwmin_p += offset ;
+        ccmin_p += offset ;
+        cmax_p += offset ;
+        cwmax_p += offset ;
+        ccmax_p += offset ;
+      }
+      grid_p -= origin ;
+      wgt_p -= origin ;
+      cmin_p -= origin ;
+      cwmin_p -= origin ;
+      ccmin_p -= origin ;
+      cmax_p -= origin ;
+      cwmax_p -= origin ;
+      ccmax_p -= origin ;
+    }
+    grid_p++ ;
+    wgt_p++ ;
+    npts_p++ ;
+    cmin_p++ ;
+    cwmin_p++ ;
+    ccmin_p++ ;
+    cmax_p++ ;
+    cwmax_p++ ;
+    ccmax_p++ ;
+  }
+  grid_p -= offset ;
+  wgt_p -= offset ;
+  npts_p -= offset ;
+  cmin_p -= offset ;
+  cwmin_p -= offset ;
+  ccmin_p -= offset ;
+  cmax_p -= offset ;
+  cwmax_p -= offset ;
+  ccmax_p -= offset ;  
+
+  // finalization
+  grid.putStorage( grid_p, delG ) ;
+  weight.putStorage( wgt_p, delW ) ;
+  npoints.freeStorage( npts_p, delNP ) ;
+  clipmin.freeStorage( cmin_p, delCMin ) ;
+  clipwmin.freeStorage( cwmin_p, delCWMin ) ;
+  clipcmin.freeStorage( ccmin_p, delCCMin ) ;
+  clipmax.freeStorage( cmax_p, delCMax ) ;
+  clipwmax.freeStorage( cwmax_p, delCWMax ) ;
+  clipcmax.freeStorage( ccmax_p, delCCMax ) ;
+}
+
+void STGrid::initPol( Int ipol ) 
+{
+  LogIO os( LogOrigin("STGrid","initPol",WHERE) ) ;
+  if ( npolOrg_ == 1 ) {
+    os << "single polarization data." << LogIO::POST ;
+    ptab_ = tab_ ;
+  }
+  else 
+    ptab_ = tab_( tab_.col("POLNO") == pollist_[ipol] ) ;
+
+  attach( ptab_ ) ;
+}
+
+void STGrid::initTable( uInt idx ) 
+{
+  tab_ = tableList_[idx] ;
+  nrow_ = rows_[idx] ;
+  updateChunkShape() ;
+}
+
+void STGrid::setData( Array<Complex> &gdata,
                       Array<Float> &gwgt )
 {
+  // 2011/12/20 TN
+  // gwgt and data_ share storage
   LogIO os( LogOrigin("STGrid","setData",WHERE) ) ;
   double t0, t1 ;
   t0 = mathutil::gettimeofday_sec() ;
-  data.resize( gdata.shape() ) ;
-  uInt len = data.nelements() ;
-  Float *w0_p ;
-  const Float *w1_p, *w2_p ;
-  Bool b0, b1, b2 ;
-  Float *data_p = data.getStorage( b0 ) ;
-  const Float *gdata_p = gdata.getStorage( b1 ) ;
-  const Float *gwgt_p = gwgt.getStorage( b2 ) ;
-  w0_p = data_p ;
+  uInt len = data_.nelements() ;
+  const Complex *w1_p ;
+  Float *w2_p ;
+  Bool b1, b2 ;
+  const Complex *gdata_p = gdata.getStorage( b1 ) ;
+  Float *gwgt_p = data_.getStorage( b2 ) ;
   w1_p = gdata_p ;
   w2_p = gwgt_p ;
   for ( uInt i = 0 ; i < len ; i++ ) {
-    *w0_p = (*w2_p > 0.0) ? (*w1_p / *w2_p) : 0.0 ;
-    w0_p++ ;
+    if ( *w2_p > 0.0 ) *w2_p = (*w1_p).real() / *w2_p ;
     w1_p++ ;
     w2_p++ ;
   }
-  data.putStorage( data_p, b0 ) ;
   gdata.freeStorage( gdata_p, b1 ) ;
-  gwgt.freeStorage( gwgt_p, b2 ) ;
+  data_.putStorage( gwgt_p, b2 ) ;
   t1 = mathutil::gettimeofday_sec() ;
-  os << "setData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ; 
+  os << LogIO::DEBUGGING << "setData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ; 
+}
+
+void STGrid::setupGrid() 
+{
+  Double xmin,xmax,ymin,ymax ;
+  mapExtent( xmin, xmax, ymin, ymax ) ;
+  
+  setupGrid( nxUI_, nyUI_, cellxUI_, cellyUI_, 
+             xmin, xmax, ymin, ymax, centerUI_ ) ;
 }
 
@@ -338 +1024 @@
                         String &center )
 {
+  LogIO os( LogOrigin("STGrid","setupGrid",WHERE) ) ;
   //cout << "nx=" << nx << ", ny=" << ny << endl ;
 
@@ -371 +1058 @@
     center_(0) = xcen.getValue( "rad" ) ;
     center_(1) = ycen.getValue( "rad" ) ;
-  }
-
+    double base = 0.5 * (xmin + xmax) ;
+    int maxrotate = 1 ;
+    int nelem = 2 * maxrotate + 1 ;
+    double *sep = new double[nelem] ;
+    for ( int i = 0 ; i < nelem ; i++ )
+      sep[i] = abs(base - center_[0] - (i-maxrotate) * C::_2pi) ;
+//     os << "sep[0]=" << sep[0] << endl  
+//        << "sep[1]=" << sep[1] << endl
+//        << "sep[2]=" << sep[2] << LogIO::POST ;
+    int idx = 0 ;
+    base = sep[0] ;
+    int nrotate = 0 ;
+    while ( idx < nelem ) {
+      if ( base > sep[idx] ) {
+        base = sep[idx] ;
+        nrotate = idx ;
+      }
+      idx++ ;
+    }
+    delete sep ;
+    nrotate -= maxrotate ;
+//     os << "nrotate = " << nrotate << LogIO::POST ;
+    center_[0] += nrotate * C::_2pi ;
+  }
+//   os << "xmin=" << xmin << LogIO::POST ;
+//   os << "center_=" << center_ << LogIO::POST ;
+
+
+  nx_ = nx ;
+  ny_ = ny ;
+  if ( nx < 0 && ny > 0 ) {
+    nx_ = ny ;
+    ny_ = ny ;
+  }
+  if ( ny < 0 && nx > 0 ) {
+    nx_ = nx ;
+    ny_ = nx ;
+  }
 
   //Double wx = xmax - xmin ;
@@ -381 +1104 @@
   wx *= 1.10 ;
   wy *= 1.10 ;
+
   Quantum<Double> qcellx ;
   Quantum<Double> qcelly ;
-  nx_ = nx ;
-  ny_ = ny ;
-  if ( nx < 0 && ny > 0 ) {
-    nx_ = ny ;
-    ny_ = ny ;
-  }
-  if ( ny < 0 && nx > 0 ) {
-    nx_ = nx ;
-    ny_ = nx ;
-  }
   //cout << "nx_ = " << nx_ << ",  ny_ = " << ny_ << endl ;
   if ( cellx.size() != 0 && celly.size() != 0 ) {
@@ -399 +1113 @@
   }
   else if ( celly.size() != 0 ) {
-    cout << "Using celly to x-axis..." << endl ;
+    os << "Using celly to x-axis..." << LogIO::POST ;
     readQuantity( qcelly, celly ) ;
     qcellx = qcelly ;
   }
   else if ( cellx.size() != 0 ) {
-    cout << "Using cellx to y-axis..." << endl ;
+    os << "Using cellx to y-axis..." << LogIO::POST ;
     readQuantity( qcellx, cellx ) ;
     qcelly = qcellx ;
@@ -410 +1124 @@
   else {
     if ( nx_ < 0 ) {
-      cout << "No user preference in grid setting. Using default..." << endl ;
+      os << "No user preference in grid setting. Using default..." << LogIO::POST ;
       readQuantity( qcellx, "1.0arcmin" ) ;
       qcelly = qcellx ;
     }
     else {
+      if ( wx == 0.0 ) {
+        os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
+        wx = 0.00290888 ;
+      }
+      if ( wy == 0.0 ) {
+        os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
+        wy = 0.00290888 ;
+      }
       qcellx = Quantum<Double>( wx/nx_, "rad" ) ;
       qcelly = Quantum<Double>( wy/ny_, "rad" ) ;
@@ -420 +1142 @@
   }
   cellx_ = qcellx.getValue( "rad" ) ;
+  // DEC correction
+  cellx_ /= cos( center_[1] ) ;
   celly_ = qcelly.getValue( "rad" ) ;
+  //os << "cellx_=" << cellx_ << ", celly_=" << celly_ << ", cos("<<center_(1)<<")=" << cos(center_(1)) << LogIO::POST ;
   if ( nx_ < 0 ) {
+    if ( wx == 0.0 ) {
+      os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
+      wx = 0.00290888 ;
+    }
+    if ( wy == 0.0 ) {
+      os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
+      wy = 0.00290888 ;
+    }
     nx_ = Int( ceil( wx/cellx_ ) ) ;
     ny_ = Int( ceil( wy/celly_ ) ) ;
@@ -427 +1160 @@
 }
 
-void STGrid::selectData( Table &tab )
-{
+void STGrid::mapExtent( Double &xmin, Double &xmax, 
+                        Double &ymin, Double &ymax ) 
+{
+  //LogIO os( LogOrigin("STGrid","mapExtent",WHERE) ) ;
+  directionCol_.attach( tableList_[0], "DIRECTION" ) ;
+  Matrix<Double> direction = directionCol_.getColumn() ;
+  //os << "dirCol.nrow() = " << dirCol.nrow() << LogIO::POST ;
+  minMax( xmin, xmax, direction.row( 0 ) ) ;
+  minMax( ymin, ymax, direction.row( 1 ) ) ;
+  Double amin, amax, bmin, bmax ;
+  for ( uInt i = 1 ; i < nfile_ ; i++ ) {
+    directionCol_.attach( tableList_[i], "DIRECTION" ) ;
+    direction.assign( directionCol_.getColumn() ) ;
+    //os << "dirCol.nrow() = " << dirCol.nrow() << LogIO::POST ;
+    minMax( amin, amax, direction.row( 0 ) ) ;
+    minMax( bmin, bmax, direction.row( 1 ) ) ;
+    xmin = min( xmin, amin ) ;
+    xmax = max( xmax, amax ) ;
+    ymin = min( ymin, bmin ) ;
+    ymax = max( ymax, bmax ) ;
+  }
+  //os << "(xmin,xmax)=(" << xmin << "," << xmax << ")" << LogIO::POST ; 
+  //os << "(ymin,ymax)=(" << ymin << "," << ymax << ")" << LogIO::POST ; 
+}
+
+void STGrid::selectData()
+{
+  LogIO os( LogOrigin("STGrid","selectData",WHERE) ) ;    
   Int ifno = ifno_ ;
-  Table taborg( infile_ ) ;
-  if ( ifno == -1 ) {
-    LogIO os( LogOrigin("STGrid","selectData",WHERE) ) ;
-//     os << LogIO::SEVERE
-//        << "Please set IFNO before actual gridding" 
-//        << LogIO::EXCEPTION ;
+  tableList_.resize( nfile_ ) ;
+  if ( ifno_ == -1 ) {
+    Table taborg( infileList_[0] ) ;
     ROScalarColumn<uInt> ifnoCol( taborg, "IFNO" ) ;
-    ifno = ifnoCol( 0 ) ;
+    ifno_ = ifnoCol( 0 ) ;
     os << LogIO::WARN
-       << "IFNO is not given. Using default IFNO: " << ifno << LogIO::POST ;
-  }
-//   tab = taborg( taborg.col("IFNO") == ifno ) ;
-  TableExprNode node ;
-  node = taborg.col("IFNO") == ifno ;
-  if ( scanlist_.size() > 0 ) {
-    node = node && taborg.col("SCANNO").in( scanlist_ ) ;
-  }
-  tab = taborg( node ) ;
-  if ( tab.nrow() == 0 ) {
-    LogIO os( LogOrigin("STGrid","selectData",WHERE) ) ;
-    os << LogIO::SEVERE
-       << "No corresponding rows for given selection: IFNO " << ifno 
-       << " SCANNO " << scanlist_ 
-       << LogIO::EXCEPTION ;
-  }
-}
-
-void STGrid::getData( Cube<Float> &spectra,
-                      Matrix<Double> &direction,
-                      Cube<uChar> &flagtra,
-                      Matrix<uInt> &rflag,
-                      Matrix<Float> &weight ) 
-{
-  Table tab ;
-  selectData( tab ) ;
-  updatePolList( tab ) ;
-//   cout << "npol_ = " << npol_ << endl ;
-//   cout << "nchan_ = " << nchan_ << endl ;
-//   cout << "nrow_ = " << nrow_ << endl ;
-  spectra.resize( npol_, nchan_, nrow_ ) ;
-  flagtra.resize( npol_, nchan_, nrow_ ) ;
-  rflag.resize( npol_, nrow_ ) ;
-  Cube<Float> tsys( npol_, nchan_, nrow_ ) ;
-  Matrix<Double> tint( npol_, nrow_ ) ;
-  // boolean for pointer access
-  Bool bsp, bfl, bfr, bts, bti ;
-  // pointer to the data
-  Float *sp_p = spectra.getStorage( bsp ) ;
-  uChar *fl_p = flagtra.getStorage( bfl ) ;
-  uInt *fr_p = rflag.getStorage( bfr ) ;
-  Float *ts_p = tsys.getStorage( bts ) ;
-  Double *ti_p = tint.getStorage( bti ) ;
-  // working pointer
-  Float *wsp_p = sp_p ;
-  uChar *wfl_p = fl_p ;
-  uInt *wfr_p = fr_p ;
-  Float *wts_p = ts_p ;
-  Double *wti_p = ti_p ;
-  uInt len = nchan_ * nrow_ ;
-  IPosition mshape( 2, nchan_, nrow_ ) ;
-  IPosition vshape( 1, nrow_ ) ;
-  for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
-    Table subt = tab( tab.col("POLNO") == pollist_[ipol] ) ;
-    ROArrayColumn<Float> spectraCol( subt, "SPECTRA" ) ;
-    ROArrayColumn<Double> directionCol( subt, "DIRECTION" ) ;
-    ROArrayColumn<uChar> flagtraCol( subt, "FLAGTRA" ) ;
-    ROScalarColumn<uInt> rflagCol( subt, "FLAGROW" ) ;
-    ROArrayColumn<Float> tsysCol( subt, "TSYS" ) ;
-    ROScalarColumn<Double> tintCol( subt, "INTERVAL" ) ;
-    Matrix<Float> spSlice( mshape, wsp_p, SHARE ) ; 
-    Matrix<uChar> flSlice( mshape, wfl_p, SHARE ) ;
-    Vector<uInt> frSlice( vshape, wfr_p, SHARE ) ;
-    spectraCol.getColumn( spSlice ) ;
-    flagtraCol.getColumn( flSlice ) ;
-    rflagCol.getColumn( frSlice ) ;
-    if ( ipol == 0 )
-      directionCol.getColumn( direction ) ;
-    Vector<Float> tmpF = tsysCol( 0 ) ;
-    Vector<Double> tmpD( vshape, wti_p, SHARE ) ;
-    Matrix<Float> tsSlice( mshape, wts_p, SHARE ) ;
-    if ( tmpF.nelements() == (uInt)nchan_ ) {
-      tsysCol.getColumn( tsSlice ) ;
+       << "IFNO is not given. Using default IFNO: " << ifno_ << LogIO::POST ;
+  }
+  for ( uInt i = 0 ; i < nfile_ ; i++ ) {
+    Table taborg( infileList_[i] ) ;
+    TableExprNode node ;
+    if ( ifno != -1 || isMultiIF( taborg ) ) {
+      os << "apply selection on IFNO" << LogIO::POST ;
+      node = taborg.col("IFNO") == ifno_ ;
+    }
+    if ( scanlist_.size() > 0 ) {
+      os << "apply selection on SCANNO" << LogIO::POST ;
+      node = node && taborg.col("SCANNO").in( scanlist_ ) ;
+    }
+    if ( node.isNull() ) {
+      tableList_[i] = taborg ;
     }
     else {
-      tsSlice = tmpF( 0 ) ;
-    }
-    tintCol.getColumn( tmpD ) ;
-
-    wsp_p += len ;
-    wfl_p += len ;
-    wfr_p += nrow_ ;
-    wts_p += len ;
-    wti_p += nrow_ ;
-  }
-  spectra.putStorage( sp_p, bsp ) ;
-  flagtra.putStorage( fl_p, bfl ) ;
-  rflag.putStorage( fr_p, bfr ) ;
-  tsys.putStorage( ts_p, bts ) ;
-  tint.putStorage( ti_p, bti ) ;
-
+      tableList_[i] = taborg( node ) ;
+    }
+    os << LogIO::DEBUGGING << "tableList_[" << i << "].nrow()=" << tableList_[i].nrow() << LogIO::POST ;
+    if ( tableList_[i].nrow() == 0 ) {
+      os << LogIO::SEVERE
+         << "No corresponding rows for given selection: IFNO " << ifno_ ;
+      if ( scanlist_.size() > 0 )
+        os << " SCANNO " << scanlist_ ;
+      os << LogIO::EXCEPTION ;
+    }
+  }
+}
+
+Bool STGrid::isMultiIF( Table &tab ) 
+{
+  ROScalarColumn<uInt> ifnoCol( tab, "IFNO" ) ;
+  Vector<uInt> ifnos = ifnoCol.getColumn() ;
+  return anyNE( ifnos, ifnos[0] ) ;
+}
+
+void STGrid::attach( Table &tab )
+{
+  // attach to table
+  spectraCol_.attach( tab, "SPECTRA" ) ;
+  flagtraCol_.attach( tab, "FLAGTRA" ) ;
+  directionCol_.attach( tab, "DIRECTION" ) ;
+  flagRowCol_.attach( tab, "FLAGROW" ) ;
+  tsysCol_.attach( tab, "TSYS" ) ;
+  intervalCol_.attach( tab, "INTERVAL" ) ;
+}
+
+Int STGrid::getDataChunk(
+                         IPosition const &wshape,
+                         IPosition const &vshape,
+                         IPosition const &dshape,
+                         Array<Complex> &spectra,
+                         Array<Double> &direction,
+                         Array<Int> &flagtra,
+                         Array<Int> &rflag,
+                         Array<Float> &weight ) 
+{
+  LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
+
+  Array<Float> spectraF_(wshape);
+  Array<uChar> flagtraUC_(wshape);
+  Array<uInt> rflagUI_(vshape);
+  Int nrow = getDataChunk( spectraF_, direction, flagtraUC_, rflagUI_, weight ) ;
+  if ( nrow < nchunk_ ) {
+    spectra.resize( spectraF_.shape() ) ;
+    flagtra.resize( flagtraUC_.shape() ) ;
+    rflag.resize( rflagUI_.shape() ) ;
+  }
+  double t0, t1 ;
+  t0 = mathutil::gettimeofday_sec() ;
+  convertArray( spectra, spectraF_ ) ;
+  toInt( flagtraUC_, flagtra ) ;
+  toInt( rflagUI_, rflag ) ;
+  t1 = mathutil::gettimeofday_sec() ;
+  eToInt = t1 - t0 ;
+  
+  return nrow ;
+}
+
+#if 0
+Int STGrid::getDataChunk( Array<Complex> &spectra,
+                          Array<Double> &direction,
+                          Array<Int> &flagtra,
+                          Array<Int> &rflag,
+                          Array<Float> &weight ) 
+{
+  LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
+  Int nrow = getDataChunk( spectraF_, direction, flagtraUC_, rflagUI_, weight ) ;
+  if ( nrow < nchunk_ ) {
+    spectra.resize( spectraF_.shape() ) ;
+    flagtra.resize( flagtraUC_.shape() ) ;
+    rflag.resize( rflagUI_.shape() ) ;
+  }
+  double t0, t1 ;
+  t0 = mathutil::gettimeofday_sec() ;
+  convertArray( spectra, spectraF_ ) ;
+  toInt( flagtraUC_, flagtra ) ;
+  toInt( rflagUI_, rflag ) ;
+  t1 = mathutil::gettimeofday_sec() ;
+  eToInt = t1 - t0 ;
+  
+  return nrow ;
+}
+#endif
+
+Int STGrid::getDataChunk( Array<Float> &spectra,
+                          Array<Double> &direction,
+                          Array<uChar> &flagtra,
+                          Array<uInt> &rflag,
+                          Array<Float> &weight ) 
+{
+  LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
+  Int nrow = spectra.shape()[1] ;
+  Int remainingRow = nrow_ - nprocessed_ ;
+  if ( remainingRow < nrow ) {
+    nrow = remainingRow ;
+    IPosition mshape( 2, nchan_, nrow ) ;
+    IPosition vshape( 1, nrow ) ;
+    spectra.resize( mshape ) ;
+    flagtra.resize( mshape ) ;
+    direction.resize( IPosition(2,2,nrow) ) ;
+    rflag.resize( vshape ) ;
+    weight.resize( mshape ) ;
+  }
+  // 2011/12/22 TN
+  // tsys shares its storage with weight
+  Array<Float> tsys( weight ) ;
+  Array<Double> tint( rflag.shape() ) ;
+
+  Vector<uInt> rflagVec( rflag ) ;
+  Vector<Double> tintVec( tint ) ;
+
+  RefRows rows( nprocessed_, nprocessed_+nrow-1, 1 ) ;
+  //os<<LogIO::DEBUGGING<<"nprocessed_="<<nprocessed_<<": rows.nrows()="<<rows.nrows()<<LogIO::POST ;
+  spectraCol_.getColumnCells( rows, spectra ) ;
+  flagtraCol_.getColumnCells( rows, flagtra ) ;
+  directionCol_.getColumnCells( rows, direction ) ;
+  flagRowCol_.getColumnCells( rows, rflagVec ) ;
+  intervalCol_.getColumnCells( rows, tintVec ) ;
+  Vector<Float> tsysTemp = tsysCol_( nprocessed_ ) ;
+  if ( tsysTemp.nelements() == (uInt)nchan_ )
+    tsysCol_.getColumnCells( rows, tsys ) ;
+  else
+    tsys = tsysTemp[0] ;
+
+  double t0,t1 ;
+  t0 = mathutil::gettimeofday_sec() ;
   getWeight( weight, tsys, tint ) ;
-}
-
-void STGrid::updatePolList( Table &tab ) 
-{
-  ROScalarColumn<uInt> polnoCol( tab, "POLNO" ) ;
+  t1 = mathutil::gettimeofday_sec() ;
+  eGetWeight += t1-t0 ;
+
+  nprocessed_ += nrow ;
+  
+  return nrow ;
+}
+
+void STGrid::setupArray() 
+{
+  LogIO os( LogOrigin("STGrid","setupArray",WHERE) ) ;
+  ROScalarColumn<uInt> polnoCol( tableList_[0], "POLNO" ) ;
   Vector<uInt> pols = polnoCol.getColumn() ;
+  //os << pols << LogIO::POST ;
   Vector<uInt> pollistOrg ;
-  uInt npolOrg = 0 ;
+  npolOrg_ = 0 ;
   uInt polno ;
   for ( uInt i = 0 ; i < polnoCol.nrow() ; i++ ) {
@@ -544 +1365 @@
     polno = pols( i ) ; 
     if ( allNE( pollistOrg, polno ) ) {
-      pollistOrg.resize( npolOrg+1, True ) ;
-      pollistOrg[npolOrg] = polno ;
-      npolOrg++ ;
+      pollistOrg.resize( npolOrg_+1, True ) ;
+      pollistOrg[npolOrg_] = polno ;
+      npolOrg_++ ;
     }
   }
@@ -565 +1386 @@
   npol_ = pollist_.size() ;
   if ( npol_ == 0 ) {
-    LogIO os( LogOrigin("STGrid","updatePolList",WHERE) ) ;
     os << LogIO::SEVERE << "Empty pollist" << LogIO::EXCEPTION ;
   }
-  nrow_ = tab.nrow() / npolOrg ;
-  ROArrayColumn<uChar> tmpCol( tab, "FLAGTRA" ) ;
-  nchan_ = tmpCol( 0 ).nelements() ;
-//   LogIO os( LogOrigin("STGrid","updatePolList",WHERE) ) ;
+  rows_.resize( nfile_ ) ;
+  for ( uInt i = 0 ; i < nfile_ ; i++ ) {
+    rows_[i] = tableList_[i].nrow() / npolOrg_ ;
+    //if ( nrow_ < rows_[i] ) 
+    //  nrow_ = rows_[i] ;
+  }
+  flagtraCol_.attach( tableList_[0], "FLAGTRA" ) ;
+  nchan_ = flagtraCol_( 0 ).nelements() ;
 //   os << "npol_ = " << npol_ << "(" << pollist_ << ")" << endl 
 //      << "nchan_ = " << nchan_ << endl 
@@ -577 +1401 @@
 }
 
-void STGrid::getWeight( Matrix<Float> &w,
-                        Cube<Float> &tsys,
-                        Matrix<Double> &tint ) 
+void STGrid::getWeight( Array<Float> &w,
+                              Array<Float> &tsys,
+                              Array<Double> &tint ) 
 {
   LogIO os( LogOrigin("STGrid","getWeight",WHERE) ) ;
-  double t0, t1 ;
-  t0 = mathutil::gettimeofday_sec() ;
-  // resize
-  w.resize( nchan_, nrow_ ) ;
+
+  // 2011/12/22 TN
+  // w (weight) and tsys share storage
+  IPosition refShape = tsys.shape() ;
+  Int nchan = refShape[0] ;
+  Int nrow = refShape[1] ;
+//   os << "nchan=" << nchan << ", nrow=" << nrow << LogIO::POST ;
+//   os << "w.shape()=" << w.shape() << endl
+//      << "tsys.shape()=" << tsys.shape() << endl
+//      << "tint.shape()=" << tint.shape() << LogIO::POST ;
 
   // set weight
-  Bool warn = False ;
   if ( wtype_.compare( "UNIFORM" ) == 0 ) {
     w = 1.0 ;
   }
   else if ( wtype_.compare( "TINT" ) == 0 ) {
-    if ( npol_ > 1 ) warn = True ;
     Bool b0, b1 ;
     Float *w_p = w.getStorage( b0 ) ;
@@ -599 +1427 @@
     const Double *ti_p = tint.getStorage( b1 ) ;
     const Double *w1_p = ti_p ;
-    for ( Int irow = 0 ; irow < nrow_ ; irow++ ) {
-      Float val = (Float)(polMean( w1_p )) ;
-      for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
-        *w0_p = val ;
+    for ( Int irow = 0 ; irow < nrow ; irow++ ) {
+      for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
+        *w0_p = *w1_p ;
         w0_p++ ;
       }
+      w1_p++ ;
     }
     w.putStorage( w_p, b0 ) ;
@@ -610 +1438 @@
   }
   else if ( wtype_.compare( "TSYS" ) == 0 ) {
-    if ( npol_ > 1 ) warn = True ;
-    Bool b0, b1 ;
+    Bool b0 ;
     Float *w_p = w.getStorage( b0 ) ;
     Float *w0_p = w_p ;
-    const Float *ts_p = tsys.getStorage( b1 ) ;
-    const Float *w1_p = ts_p ;
-    for ( Int irow = 0 ; irow < nrow_ ; irow++ ) {
-      for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
-        Float val = polMean( w1_p ) ;
-        *w0_p = 1.0 / ( val * val ) ;
+    for ( Int irow = 0 ; irow < nrow ; irow++ ) {
+      for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
+        Float temp = *w0_p ;
+        *w0_p = 1.0 / ( temp * temp ) ;
         w0_p++ ;
       }
     }
     w.putStorage( w_p, b0 ) ;
-    tsys.freeStorage( ts_p, b1 ) ;
   }
   else if ( wtype_.compare( "TINTSYS" ) == 0 ) {
-    if ( npol_ > 1 ) warn = True ;
-    Bool b0, b1, b2 ;
+    Bool b0, b1 ;
     Float *w_p = w.getStorage( b0 ) ;
     Float *w0_p = w_p ;
     const Double *ti_p = tint.getStorage( b1 ) ;
     const Double *w1_p = ti_p ;
-    const Float *ts_p = tsys.getStorage( b2 ) ;
-    const Float *w2_p = ts_p ;
-    for ( Int irow = 0 ; irow < nrow_ ; irow++ ) {
-      Float interval = (Float)(polMean( w1_p )) ;
-      for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
-        Float temp = polMean( w2_p ) ;
+    for ( Int irow = 0 ; irow < nrow ; irow++ ) {
+      Float interval = *w1_p ;
+      for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
+        Float temp = *w0_p ;
         *w0_p = interval / ( temp * temp ) ;
         w0_p++ ;
       }
+      w1_p++ ;
     }
     w.putStorage( w_p, b0 ) ;
     tint.freeStorage( ti_p, b1 ) ;
-    tsys.freeStorage( ts_p, b2 ) ;
   }
   else {
     //LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
-    os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
+    //os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
     w = 1.0 ;
   }
-
-  if ( npol_ > 1 ) {
-    //LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
-    os << LogIO::WARN << "STGrid doesn't support assigning polarization-dependent weight. Use averaged weight over polarization." << LogIO::POST ;
-  }
-  t1 = mathutil::gettimeofday_sec() ;
-  os << "getWeight: elapsed time is " << t1-t0 << " sec" << LogIO::POST ;
-}
-
-Float STGrid::polMean( const Float *p ) 
-{
-  Float v = 0.0 ;
-  for ( Int i = 0 ; i < npol_ ; i++ ) {
-    v += *p ;
-    p++ ;
-  }
-  v /= npol_ ;
-  return v ;
-}
-
-Double STGrid::polMean( const Double *p ) 
-{
-  Double v = 0.0 ;
-  for ( Int i = 0 ; i < npol_ ; i++ ) {
-    v += *p ;
-    p++ ;
-  }
-  v /= npol_ ;
-  return v ;
-}
-
-void STGrid::toInt( Array<uChar> *u, Array<Int> *v ) 
-{
-  uInt len = u->nelements() ;
+}
+
+void STGrid::toInt( Array<uChar> &u, Array<Int> &v ) 
+{
+  uInt len = u.nelements() ;
   Int *int_p = new Int[len] ;
   Bool deleteIt ;
-  const uChar *data_p = u->getStorage( deleteIt ) ;
+  const uChar *data_p = u.getStorage( deleteIt ) ;
   Int *i_p = int_p ;
   const uChar *u_p = data_p ;
@@ -696 +1488 @@
     u_p++ ;
   }
-  u->freeStorage( data_p, deleteIt ) ;
-  v->takeStorage( u->shape(), int_p, TAKE_OVER ) ;
-}
-
-void STGrid::toInt( Array<uInt> *u, Array<Int> *v ) 
-{
-  uInt len = u->nelements() ;
+  u.freeStorage( data_p, deleteIt ) ;
+  v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
+}
+
+void STGrid::toInt( Array<uInt> &u, Array<Int> &v ) 
+{
+  uInt len = u.nelements() ;
   Int *int_p = new Int[len] ;
   Bool deleteIt ;
-  const uInt *data_p = u->getStorage( deleteIt ) ;
+  const uInt *data_p = u.getStorage( deleteIt ) ;
   Int *i_p = int_p ;
   const uInt *u_p = data_p ;
@@ -713 +1505 @@
     u_p++ ;
   }
-  u->freeStorage( data_p, deleteIt ) ;
-  v->takeStorage( u->shape(), int_p, TAKE_OVER ) ;
-}
-
-void STGrid::toPixel( Matrix<Double> &world, Matrix<Double> &pixel )
+  u.freeStorage( data_p, deleteIt ) ;
+  v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
+}
+
+void STGrid::toPixel( Array<Double> &world, Array<Double> &pixel )
 {
   // gridding will be done on (nx_+2*convSupport_) x (ny_+2*convSupport_) 
   // grid plane to avoid unexpected behavior on grid edge
-  Vector<Double> pixc( 2 ) ;
-  pixc(0) = Double( nx_-1 ) * 0.5 ;
-  pixc(1) = Double( ny_-1 ) * 0.5 ;
-//   pixc(0) = Double( nx_+2*convSupport_-1 ) * 0.5 ;
-//   pixc(1) = Double( ny_+2*convSupport_-1 ) * 0.5 ;
+  Block<Double> pixc( 2 ) ;
+  pixc[0] = Double( nx_-1 ) * 0.5 ;
+  pixc[1] = Double( ny_-1 ) * 0.5 ;
+//   pixc[0] = Double( nx_+2*convSupport_-1 ) * 0.5 ;
+//   pixc[1] = Double( ny_+2*convSupport_-1 ) * 0.5 ;
   uInt nrow = world.shape()[1] ;
-  Vector<Double> cell( 2 ) ;
-  cell(0) = cellx_ ;
-  cell(1) = celly_ ;
-  for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
-    for ( uInt i = 0 ; i < 2 ; i++ ) {
-      pixel( i, irow ) = pixc(i) + ( world(i, irow) - center_(i) ) / cell(i) ;
-    }
-  }
-//   String gridfile = "grid."+convType_+"."+String::toString(convSupport_)+".dat" ; 
-//   ofstream ofs( gridfile.c_str(), ios::out ) ; 
-//   ofs << "center " << center_(0) << " " << pixc(0) 
-//       << " " << center_(1) << " " << pixc(1) << endl ;
-//   for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
-//     ofs << irow ;
-//     for ( uInt i = 0 ; i < 2 ; i++ ) {
-//       ofs << " " << world(i, irow) << " " << pixel(i, irow) ;
-//     }
-//     ofs << endl ;
-//   }
-//   ofs.close() ;
+  Bool bw, bp ;
+  const Double *w_p = world.getStorage( bw ) ;
+  Double *p_p = pixel.getStorage( bp ) ;
+  const Double *ww_p = w_p ;
+  Double *wp_p = p_p ;
+  for ( uInt i = 0 ; i < nrow ; i++ ) {
+    *wp_p = pixc[0] + ( *ww_p - center_[0] ) / cellx_ ;
+    wp_p++ ;
+    ww_p++ ;
+    *wp_p = pixc[1] + ( *ww_p - center_[1] ) / celly_ ;
+    wp_p++ ;
+    ww_p++ ;
+  }
+  world.freeStorage( w_p, bw ) ;
+  pixel.putStorage( p_p, bp ) ;  
 }
 
@@ -812 +1599 @@
     // to take into account Gaussian tail
     if ( convSupport_ < 0 )
-      convSupport_ = 12 ; // 3 * 4
+      convSupport_ = 4 ; // 1 * 4
     else {
       convSupport_ = userSupport_ * 4 ;
@@ -839 +1626 @@
   String outfile_ ;
   if ( outfile.size() == 0 ) {
-    if ( infile_.lastchar() == '/' ) {
-      outfile_ = infile_.substr( 0, infile_.size()-1 ) ;
+    if ( infileList_[0].lastchar() == '/' ) {
+      outfile_ = infileList_[0].substr( 0, infileList_[0].size()-1 ) ;
     }
     else {
-      outfile_ = infile_ ;
+      outfile_ = infileList_[0] ;
     }
     outfile_ += ".grid" ;
@@ -897 +1684 @@
 
   t1 = mathutil::gettimeofday_sec() ;
-  os << "saveData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ; 
-  
+  os << LogIO::DEBUGGING << "saveData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ; 
+
+  fillMainColumns( tab ) ;
+
   return outfile_ ;
 }
@@ -904 +1693 @@
 void STGrid::prepareTable( Table &tab, String &name ) 
 {
-  Table t( infile_, Table::Old ) ;
+  Table t( infileList_[0], Table::Old ) ;
   t.deepCopy( name, Table::New, False, t.endianFormat(), True ) ;
   tab = Table( name, Table::Update ) ;
-}
-}
+  // 2012/02/13 TN
+  // explicitly copy subtables since no rows including subtables are 
+  // copied by Table::deepCopy with noRows=True
+  TableCopy::copySubTables( tab, t ) ;
+}
+
+void STGrid::fillMainColumns( Table &tab ) 
+{
+  // values for fill
+  Table t( infileList_[0], Table::Old ) ;
+  Table tsel = t( t.col( "IFNO" ) == (uInt)ifno_, 1 ) ;
+  ROTableRow row( tsel ) ;
+  row.get( 0 ) ;
+  const TableRecord &rec = row.record() ;
+  uInt freqId = rec.asuInt( "FREQ_ID" ) ;
+  uInt molId = rec.asuInt( "MOLECULE_ID" ) ;
+  uInt tcalId = rec.asuInt( "TCAL_ID" ) ;
+  uInt focusId = rec.asuInt( "FOCUS_ID" ) ;
+  uInt weatherId = rec.asuInt( "WEATHER_ID" ) ;
+  String srcname = rec.asString( "SRCNAME" ) ;
+  String fieldname = rec.asString( "FIELDNAME" ) ;
+  Vector<Float> defaultTsys( 1, 1.0 ) ;
+  // @todo how to set flagtra for gridded spectra?
+  Vector<uChar> flagtra = rec.asArrayuChar( "FLAGTRA" ) ;
+  flagtra = (uChar)0 ;
+  Float opacity = rec.asFloat( "OPACITY" ) ;
+  Double srcvel = rec.asDouble( "SRCVELOCITY" ) ;
+  Vector<Double> srcpm = rec.asArrayDouble( "SRCPROPERMOTION" ) ;
+  Vector<Double> srcdir = rec.asArrayDouble( "SRCDIRECTION" ) ;
+  Vector<Double> scanrate = rec.asArrayDouble( "SCANRATE" ) ;
+  Double time = rec.asDouble( "TIME" ) ;
+  Double interval = rec.asDouble( "INTERVAL" ) ;
+
+  // fill columns
+  Int nrow = tab.nrow() ;
+  ScalarColumn<uInt> scannoCol( tab, "SCANNO" ) ;
+  ScalarColumn<uInt> ifnoCol( tab, "IFNO" ) ;
+  ScalarColumn<uInt> freqIdCol( tab, "FREQ_ID" ) ;
+  ScalarColumn<uInt> molIdCol( tab, "MOLECULE_ID" ) ;
+  ScalarColumn<uInt> tcalidCol( tab, "TCAL_ID" ) ;
+  ScalarColumn<Int> fitidCol( tab, "FIT_ID" ) ;
+  ScalarColumn<uInt> focusidCol( tab, "FOCUS_ID" ) ;
+  ScalarColumn<uInt> weatheridCol( tab, "WEATHER_ID" ) ;
+  ArrayColumn<uChar> flagtraCol( tab, "FLAGTRA" ) ;
+  ScalarColumn<uInt> rflagCol( tab, "FLAGROW" ) ;
+  ArrayColumn<Float> tsysCol( tab, "TSYS" ) ;
+  ScalarColumn<String> srcnameCol( tab, "SRCNAME" ) ;
+  ScalarColumn<String> fieldnameCol( tab, "FIELDNAME" ) ;
+  ScalarColumn<Int> srctypeCol( tab, "SRCTYPE" ) ;
+  ScalarColumn<Float> opacityCol( tab, "OPACITY" ) ;
+  ScalarColumn<Double> srcvelCol( tab, "SRCVELOCITY" ) ;
+  ArrayColumn<Double> srcpmCol( tab, "SRCPROPERMOTION" ) ;
+  ArrayColumn<Double> srcdirCol( tab, "SRCDIRECTION" ) ;
+  ArrayColumn<Double> scanrateCol( tab, "SCANRATE" ) ;
+  ScalarColumn<Double> timeCol( tab, "TIME" ) ;
+  ScalarColumn<Double> intervalCol( tab, "INTERVAL" ) ;
+  for ( Int i = 0 ; i < nrow ; i++ ) {
+    scannoCol.put( i, (uInt)i ) ;
+    ifnoCol.put( i, (uInt)ifno_ ) ;
+    freqIdCol.put( i, freqId ) ;
+    molIdCol.put( i, molId ) ;
+    tcalidCol.put( i, tcalId ) ;
+    fitidCol.put( i, -1 ) ;
+    focusidCol.put( i, focusId ) ;
+    weatheridCol.put( i, weatherId ) ;
+    flagtraCol.put( i, flagtra ) ;
+    rflagCol.put( i, 0 ) ;
+    tsysCol.put( i, defaultTsys ) ;
+    srcnameCol.put( i, srcname ) ;
+    fieldnameCol.put( i, fieldname ) ;
+    srctypeCol.put( i, (Int)SrcType::PSON ) ;
+    opacityCol.put( i, opacity ) ;
+    srcvelCol.put( i, srcvel ) ;
+    srcpmCol.put( i, srcpm ) ;
+    srcdirCol.put( i, srcdir ) ;
+    scanrateCol.put( i, scanrate ) ;
+    timeCol.put( i, time ) ;
+    intervalCol.put( i, interval ) ;
+  }
+}
+
+}