source: trunk/src/STGrid.cpp @ 2397

Last change on this file since 2397 was 2397, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: Yes CAS-2816

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

Updated log messages and documentation.
File size: 52.8 KB
Line 
1#include <iostream>
2#include <fstream>
3#include <cfloat>
4
5#include <casa/BasicSL/String.h>
6#include <casa/Arrays/Vector.h>
7#include <casa/Arrays/Matrix.h>
8#include <casa/Arrays/Cube.h>
9#include <casa/Arrays/ArrayMath.h>
10#include <casa/Arrays/ArrayIter.h>
11#include <casa/Quanta/Quantum.h>
12#include <casa/Quanta/QuantumHolder.h>
13#include <casa/Utilities/CountedPtr.h>
14#include <casa/Logging/LogIO.h>
15
16#include <tables/Tables/Table.h>
17#include <tables/Tables/TableRecord.h>
18#include <tables/Tables/ExprNode.h>
19#include <tables/Tables/ScalarColumn.h>
20#include <tables/Tables/ArrayColumn.h>
21#include <tables/Tables/TableIter.h>
22
23#include <measures/Measures/MDirection.h>
24
25#include <MathUtils.h>
26
27#include "STGrid.h"
28
29using namespace std ;
30using namespace concurrent ;
31using namespace casa ;
32using namespace asap ;
33
34namespace asap {
35
36// for performance check
37double eToInt = 0.0 ;
38double eGetWeight = 0.0 ;
39
40// constructor
41STGrid::STGrid()
42  : vshape_( 1 ), wshape_( 2 ), dshape_( 2 )
43{
44  init() ;
45}
46
47STGrid::STGrid( const string infile )
48  : vshape_( 1 ), wshape_( 2 ), dshape_( 2 )
49{
50  init() ;
51
52  setFileIn( infile ) ;
53}
54
55STGrid::STGrid( const vector<string> infile )
56{
57  init() ;
58
59  setFileList( infile ) ;
60}
61
62void  STGrid::init()
63{
64  ifno_ = -1 ;
65  nx_ = -1 ;
66  ny_ = -1 ;
67  npol_ = 0 ;
68  nchan_ = 0 ;
69  nrow_ = 0 ;
70  cellx_ = 0.0 ;
71  celly_ = 0.0 ;
72  center_ = Vector<Double> ( 2, 0.0 ) ;
73  convType_ = "BOX" ;
74  wtype_ = "UNIFORM" ;
75  convSupport_ = -1 ;
76  userSupport_ = -1 ;
77  convSampling_ = 100 ;
78  nprocessed_ = 0 ;
79  nchunk_ = 0 ;
80
81  // initialize user input
82  nxUI_ = -1 ;
83  nyUI_ = -1 ;
84  cellxUI_ = "" ;
85  cellyUI_ = "" ;
86  centerUI_ = "" ;
87  doclip_ = False ;
88}
89
90void STGrid::setFileIn( const string infile )
91{
92  nfile_ = 1 ;
93  String name( infile ) ;
94  infileList_.resize( nfile_ ) ;
95  infileList_[0] = String(infile) ;
96}
97
98void STGrid::setFileList( const vector<string> infile )
99{
100  nfile_ = infile.size() ;
101  infileList_.resize( nfile_ ) ;
102  for ( uInt i = 0 ; i < nfile_ ; i++ ) {
103    infileList_[i] = infile[i] ;
104  }
105}
106
107void STGrid::setPolList( vector<unsigned int> pols )
108{
109  pollist_.assign( Vector<uInt>( pols ) ) ;
110}
111
112void STGrid::setScanList( vector<unsigned int> scans )
113{
114  scanlist_.assign( Vector<uInt>( scans ) ) ;
115}
116
117void STGrid::setWeight( const string wType )
118{
119  wtype_ = String( wType ) ;
120  wtype_.upcase() ;
121}
122
123void STGrid::defineImage( int nx,
124                          int ny,
125                          string scellx,
126                          string scelly,
127                          string scenter )
128{
129  nxUI_ = (Int)nx ;
130  nyUI_ = (Int)ny ;
131  cellxUI_ = String( scellx ) ;
132  cellyUI_ = String( scelly ) ;
133  centerUI_ = String( scenter ) ;
134}
135 
136void STGrid::setFunc( string convType,
137                      int convSupport )
138{
139  convType_ = String( convType ) ;
140  convType_.upcase() ;
141  userSupport_ = (Int)convSupport ;
142}
143
144#define NEED_UNDERSCORES
145#if defined(NEED_UNDERSCORES)
146#define ggridsd ggridsd_
147#endif
148extern "C" {
149   void ggridsd(Double*,
150                const Complex*,
151                Int*,
152                Int*,
153                Int*,
154                const Int*,
155                const Int*,
156                const Float*,
157                Int*,
158                Int*,
159                Complex*,
160                Float*,
161                Int*,
162                Int*,
163                Int *,
164                Int *,
165                Int*,
166                Int*,
167                Float*,
168                Int*,
169                Int*,
170                Double*);
171}
172void STGrid::call_ggridsd( Array<Double> &xypos,
173                           Array<Complex> &spectra,
174                           Int &nvispol,
175                           Int &nvischan,
176                           Array<Int> &flagtra,
177                           Array<Int> &flagrow,
178                           Array<Float> &weight,
179                           Int &nrow,
180                           Int &irow,
181                           Array<Complex> &gdata,
182                           Array<Float> &gwgt,
183                           Int &nx,
184                           Int &ny,
185                           Int &npol,
186                           Int &nchan,
187                           Int &support,
188                           Int &sampling,
189                           Vector<Float> &convFunc,
190                           Int *chanMap,
191                           Int *polMap )
192{
193  // parameters for gridding
194  Int idopsf = 0 ;
195  Int len = npol*nchan ;
196  Double *sumw_p = new Double[len] ;
197  {
198    Double *work_p = sumw_p ;
199    for ( Int i = 0 ; i < len ; i++ ) {
200      *work_p = 0.0 ;
201      work_p++ ;
202    }
203  }
204
205  // prepare pointer
206  Bool deletePos, deleteData, deleteWgt, deleteFlag, deleteFlagR, deleteConv, deleteDataG, deleteWgtG ;
207  Double *xy_p = xypos.getStorage( deletePos ) ;
208  const Complex *values_p = spectra.getStorage( deleteData ) ;
209  const Int *flag_p = flagtra.getStorage( deleteFlag ) ;
210  const Int *rflag_p = flagrow.getStorage( deleteFlagR ) ;
211  const Float *wgt_p = weight.getStorage( deleteWgt ) ;
212  Complex *grid_p = gdata.getStorage( deleteDataG ) ;
213  Float *wgrid_p = gwgt.getStorage( deleteWgtG ) ;
214  Float *conv_p = convFunc.getStorage( deleteConv ) ;
215
216  // pass copy of irow to ggridsd since it will be modified in theroutine
217  Int irowCopy = irow ;
218     
219  // call ggridsd
220  ggridsd( xy_p,
221           values_p,
222           &nvispol,
223           &nvischan,
224           &idopsf,
225           flag_p,
226           rflag_p,
227           wgt_p,
228           &nrow,
229           &irowCopy,
230           grid_p,
231           wgrid_p,
232           &nx,
233           &ny,
234           &npol,
235           &nchan,
236           &support,
237           &sampling,
238           conv_p,
239           chanMap,
240           polMap,
241           sumw_p ) ;
242
243  // finalization
244  xypos.putStorage( xy_p, deletePos ) ;
245  spectra.freeStorage( values_p, deleteData ) ;
246  flagtra.freeStorage( flag_p, deleteFlag ) ;
247  flagrow.freeStorage( rflag_p, deleteFlagR ) ;
248  weight.freeStorage( wgt_p, deleteWgt ) ;
249  gdata.putStorage( grid_p, deleteDataG ) ;
250  gwgt.putStorage( wgrid_p, deleteWgtG ) ;
251  convFunc.putStorage( conv_p, deleteConv ) ;
252  delete sumw_p ;
253}
254
255#define NEED_UNDERSCORES
256#if defined(NEED_UNDERSCORES)
257#define ggridsd2 ggridsd2_
258#endif
259extern "C" {
260   void ggridsd2(Double*,
261                 const Complex*,
262                 Int*,
263                 Int*,
264                 Int*,
265                 const Int*,
266                 const Int*,
267                 const Float*,
268                 Int*,
269                 Int*,
270                 Complex*,
271                 Float*,
272                 Int*,
273                 Complex*,
274                 Float*,
275                 Float*,
276                 Complex*,
277                 Float*,
278                 Float*,
279                 Int*,
280                 Int*,
281                 Int *,
282                 Int *,
283                 Int*,
284                 Int*,
285                 Float*,
286                 Int*,
287                 Int*,
288                 Double*);
289}
290void STGrid::call_ggridsd2( Array<Double> &xypos,
291                            Array<Complex> &spectra,
292                            Int &nvispol,
293                            Int &nvischan,
294                            Array<Int> &flagtra,
295                            Array<Int> &flagrow,
296                            Array<Float> &weight,
297                            Int &nrow,
298                            Int &irow,
299                            Array<Complex> &gdata,
300                            Array<Float> &gwgt,
301                            Array<Int> &npoints,
302                            Array<Complex> &clipmin,
303                            Array<Float> &clipwmin,
304                            Array<Float> &clipcmin,
305                            Array<Complex> &clipmax,
306                            Array<Float> &clipwmax,
307                            Array<Float> &clipcmax,
308                            Int &nx,
309                            Int &ny,
310                            Int &npol,
311                            Int &nchan,
312                            Int &support,
313                            Int &sampling,
314                            Vector<Float> &convFunc,
315                            Int *chanMap,
316                            Int *polMap )
317{
318  // parameters for gridding
319  Int idopsf = 0 ;
320  Int len = npol*nchan ;
321  Double *sumw_p = new Double[len] ;
322  {
323    Double *work_p = sumw_p ;
324    for ( Int i = 0 ; i < len ; i++ ) {
325      *work_p = 0.0 ;
326      work_p++ ;
327    }
328  }
329
330  // prepare pointer
331  Bool deletePos, deleteData, deleteWgt, deleteFlag, deleteFlagR, deleteConv, deleteDataG, deleteWgtG, deleteNpts, deleteCMin, deleteCWMin, deleteCCMin, deleteCMax, deleteCWMax, deleteCCMax ;
332  Double *xy_p = xypos.getStorage( deletePos ) ;
333  const Complex *values_p = spectra.getStorage( deleteData ) ;
334  const Int *flag_p = flagtra.getStorage( deleteFlag ) ;
335  const Int *rflag_p = flagrow.getStorage( deleteFlagR ) ;
336  const Float *wgt_p = weight.getStorage( deleteWgt ) ;
337  Complex *grid_p = gdata.getStorage( deleteDataG ) ;
338  Float *wgrid_p = gwgt.getStorage( deleteWgtG ) ;
339  Float *conv_p = convFunc.getStorage( deleteConv ) ;
340  Int *npts_p = npoints.getStorage( deleteNpts ) ;
341  Complex *cmin_p = clipmin.getStorage( deleteCMin ) ;
342  Float *cwmin_p = clipwmin.getStorage( deleteCWMin ) ;
343  Float *ccmin_p = clipcmin.getStorage( deleteCCMin ) ;
344  Complex *cmax_p = clipmax.getStorage( deleteCMax ) ;
345  Float *cwmax_p = clipwmax.getStorage( deleteCWMax ) ;
346  Float *ccmax_p = clipcmax.getStorage( deleteCCMax ) ;
347
348  // pass copy of irow to ggridsd since it will be modified in theroutine
349  Int irowCopy = irow ;
350     
351  // call ggridsd
352  ggridsd2( xy_p,
353            values_p,
354            &nvispol,
355            &nvischan,
356            &idopsf,
357            flag_p,
358            rflag_p,
359            wgt_p,
360            &nrow,
361            &irowCopy,
362            grid_p,
363            wgrid_p,
364            npts_p,
365            cmin_p,
366            cwmin_p,
367            ccmin_p,
368            cmax_p,
369            cwmax_p,
370            ccmax_p,
371            &nx,
372            &ny,
373            &npol,
374            &nchan,
375            &support,
376            &sampling,
377            conv_p,
378            chanMap,
379            polMap,
380            sumw_p ) ;
381
382  // finalization
383  xypos.putStorage( xy_p, deletePos ) ;
384  spectra.freeStorage( values_p, deleteData ) ;
385  flagtra.freeStorage( flag_p, deleteFlag ) ;
386  flagrow.freeStorage( rflag_p, deleteFlagR ) ;
387  weight.freeStorage( wgt_p, deleteWgt ) ;
388  gdata.putStorage( grid_p, deleteDataG ) ;
389  gwgt.putStorage( wgrid_p, deleteWgtG ) ;
390  convFunc.putStorage( conv_p, deleteConv ) ;
391  clipmin.putStorage( cmin_p, deleteCMin ) ;
392  clipwmin.putStorage( cwmin_p, deleteCWMin ) ;
393  clipcmin.putStorage( ccmin_p, deleteCCMin ) ;
394  clipmax.putStorage( cmax_p, deleteCMax ) ;
395  clipwmax.putStorage( cwmax_p, deleteCWMax ) ;
396  clipcmax.putStorage( ccmax_p, deleteCCMax ) ;
397  delete sumw_p ;
398}
399
400void STGrid::grid()
401{
402  LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
403  double t0,t1 ;
404
405  // data selection
406  t0 = mathutil::gettimeofday_sec() ;
407  selectData() ;
408  t1 = mathutil::gettimeofday_sec() ;
409  os << "selectData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
410
411  setupGrid() ;
412  setupArray() ;
413
414  if ( wtype_.compare("UNIFORM") != 0 &&
415       wtype_.compare("TINT") != 0 &&
416       wtype_.compare("TSYS") != 0 &&
417       wtype_.compare("TINTSYS") != 0 ) {
418    LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
419    os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
420    wtype_ = "UNIFORM" ;
421  }
422
423  // grid parameter
424  os << LogIO::DEBUGGING ;
425  os << "----------" << endl ;
426  os << "Data selection summary" << endl ;
427  os << "   ifno = " << ifno_ << endl ;
428  os << "   pollist = " << pollist_ << endl ;
429  os << "   scanlist = " << scanlist_ << endl ;
430  os << "----------" << endl ;
431  os << "Grid parameter summary" << endl ;
432  os << "   (nx,ny) = (" << nx_ << "," << ny_ << ")" << endl ;
433  os << "   (cellx,celly) = (" << cellx_ << "," << celly_ << ")" << endl ;
434  os << "   center = " << center_ << endl ;
435  os << "   weighting = " << wtype_ << endl ;
436  os << "   convfunc = " << convType_ << " with support " << convSupport_ << endl ;
437  os << "   doclip = " << (doclip_?"True":"False") << endl ;
438  os << "----------" << LogIO::POST ;
439  os << LogIO::NORMAL ;
440
441  Bool doAll = examine() ;
442
443  if ( doAll )
444    gridPerPol() ;
445  else if ( doclip_ )
446    gridPerRowWithClipping() ;
447  else
448    gridPerRow() ;
449
450}
451
452Bool STGrid::examine()
453{
454  // TODO: nchunk_ must be determined from nchan_, npol_, and (nx_,ny_)
455  //       by considering data size to be allocated for ggridsd input/output
456  nchunk_ = 400 ;
457  Bool b = nchunk_ >= nrow_ ;
458  nchunk_ = min( nchunk_, nrow_ ) ;
459  vshape_ = IPosition( 1, nchunk_ ) ;
460  wshape_ = IPosition( 2, nchan_, nchunk_ ) ;
461  dshape_ = IPosition( 2, 2, nchunk_ ) ;
462  return b ;
463}
464
465struct STGChunk {
466  Int nrow ;
467  Array<Complex> spectra;
468  Array<Int> flagtra;
469  Array<Int> rflag;
470  Array<Float> weight;
471  Array<Double> direction;
472  STGChunk(IPosition const &wshape, IPosition const &vshape,
473           IPosition const &dshape)
474    : spectra(wshape), flagtra(wshape), rflag(vshape), weight(wshape),
475      direction(dshape)
476  { }
477};
478
479struct STCommonData {
480  Int gnx;
481  Int gny;
482  Int *chanMap;
483  Vector<Float> convFunc ;
484  Array<Complex> gdataArrC;
485  Array<Float> gwgtArr;
486  STCommonData(IPosition const &gshape, Array<Float> const &data)
487    : gdataArrC(gshape, 0.0), gwgtArr(data) {}
488};
489
490struct STCommonDataWithClipping {
491  Int gnx;
492  Int gny;
493  Int *chanMap;
494  Vector<Float> convFunc ;
495  Array<Complex> gdataArrC;
496  Array<Float> gwgtArr;
497  Array<Int> npoints ;
498  Array<Complex> clipMin ;
499  Array<Float> clipWMin ;
500  Array<Float> clipCMin ;
501  Array<Complex> clipMax ;
502  Array<Float> clipWMax ;
503  Array<Float> clipCMax ; 
504  STCommonDataWithClipping(IPosition const &gshape,
505                           IPosition const &pshape,
506                           Array<Float> const &data)
507    : gdataArrC(gshape, 0.0),
508      gwgtArr(data),
509      npoints(pshape, 0),
510      clipMin(gshape, Complex(FLT_MAX,0.0)),
511      clipWMin(gshape, 0.0),
512      clipCMin(gshape, 0.0),
513      clipMax(gshape, Complex(-FLT_MAX,0.0)),
514      clipWMax(gshape, 0.0),
515      clipCMax(gshape, 0.0)
516  {}
517};
518
519#define DO_AHEAD 3
520
521struct STContext {
522  STCommonData &common;
523  FIFO<STGChunk *, DO_AHEAD> queue;
524  STGrid *const self;
525  const Int pol;
526  STContext(STGrid *obj, STCommonData &common, Int pol)
527    : self(obj), common(common), pol(pol) {}
528};
529
530struct STContextWithClipping {
531  STCommonDataWithClipping &common;
532  FIFO<STGChunk *, DO_AHEAD> queue;
533  STGrid *const self;
534  const Int pol;
535  STContextWithClipping(STGrid *obj, STCommonDataWithClipping &common, Int pol)
536    : self(obj), common(common), pol(pol) {}
537};
538
539
540bool STGrid::produceChunk(void *ctx) throw(PCException)
541{
542  STContext &context = *(STContext *)ctx;
543  if ( context.self->nprocessed_ >= context.self->nrow_ ) {
544    return false;
545  }
546  STGChunk *chunk = new STGChunk(context.self->wshape_,
547                                 context.self->vshape_,
548                                 context.self->dshape_);
549
550  double t0 = mathutil::gettimeofday_sec() ;
551  chunk->nrow = context.self->getDataChunk(
552        context.self->wshape_, context.self->vshape_, context.self->dshape_,
553        chunk->spectra, chunk->direction,
554        chunk->flagtra, chunk->rflag, chunk->weight);
555  double t1 = mathutil::gettimeofday_sec() ;
556  context.self->eGetData_ += t1-t0 ;
557
558  context.queue.lock();
559  context.queue.put(chunk);
560  context.queue.unlock();
561  return true;
562}
563
564void STGrid::consumeChunk(void *ctx) throw(PCException)
565{
566  STContext &context = *(STContext *)ctx;
567  STGChunk *chunk = NULL;
568  try {
569    context.queue.lock();
570    chunk = context.queue.get();
571    context.queue.unlock();
572  } catch (FullException &e) {
573    context.queue.unlock();
574    // TODO: log error
575    throw PCException();
576  }
577
578  double t0, t1 ;
579  // world -> pixel
580  Array<Double> xypos( context.self->dshape_ ) ;
581  t0 = mathutil::gettimeofday_sec() ;
582  context.self->toPixel( chunk->direction, xypos ) ;
583  t1 = mathutil::gettimeofday_sec() ;
584  context.self->eToPixel_ += t1-t0 ;
585   
586  // call ggridsd
587  Int nvispol = 1 ;
588  Int irow = -1 ;
589  t0 = mathutil::gettimeofday_sec() ;
590  context.self->call_ggridsd( xypos,
591                chunk->spectra,
592                nvispol,
593                context.self->nchan_,
594                chunk->flagtra,
595                chunk->rflag,
596                chunk->weight,
597                chunk->nrow,
598                irow,
599                context.common.gdataArrC,
600                context.common.gwgtArr,
601                context.common.gnx,
602                context.common.gny,
603                context.self->npol_,
604                context.self->nchan_,
605                context.self->convSupport_,
606                context.self->convSampling_,
607                context.common.convFunc,
608                context.common.chanMap,
609                (Int*)&context.pol ) ;
610  t1 = mathutil::gettimeofday_sec() ;
611  context.self->eGGridSD_ += t1-t0 ;
612 
613  delete chunk;
614}
615
616void STGrid::gridPerRow()
617{
618  LogIO os( LogOrigin("STGrid", "gridPerRow", WHERE) ) ;
619  double t0, t1 ;
620
621
622  // grid data
623  // Extend grid plane with convSupport_
624  //   Int gnx = nx_+convSupport_*2 ;
625  //   Int gny = ny_+convSupport_*2 ;
626  Int gnx = nx_;
627  Int gny = ny_;
628
629  IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
630  // 2011/12/20 TN
631  // data_ and gwgtArr share storage
632  data_.resize( gshape ) ;
633  data_ = 0.0 ;
634  STCommonData common = STCommonData(gshape, data_);
635  common.gnx = gnx ;
636  common.gny = gny ;
637
638  // parameters for gridding
639  Int *chanMap = new Int[nchan_] ;
640  for ( Int i = 0 ; i < nchan_ ; i++ ) {
641    chanMap[i] = i ;
642  }
643  common.chanMap = chanMap;
644
645  // convolution kernel
646  t0 = mathutil::gettimeofday_sec() ;
647  setConvFunc( common.convFunc ) ;
648  t1 = mathutil::gettimeofday_sec() ;
649  os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
650
651  // for performance check
652  eGetData_ = 0.0 ;
653  eToPixel_ = 0.0 ;
654  eGGridSD_ = 0.0 ;
655  double eInitPol = 0.0 ;
656
657  Broker broker = Broker(produceChunk, consumeChunk);
658  for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) {
659    initTable( ifile ) ;
660
661    os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ;   
662    for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
663      t0 = mathutil::gettimeofday_sec() ;
664      initPol( ipol ) ; // set ptab_ and attach()
665      t1 = mathutil::gettimeofday_sec() ;
666      eInitPol += t1-t0 ;
667     
668      STContext context(this, common, ipol);
669     
670      os << "start pol " << ipol << LogIO::POST ;
671     
672      nprocessed_ = 0 ;
673#if 1
674      broker.runProducerAsMasterThread(&context, DO_AHEAD);
675#else
676      for (;;) {
677        bool produced = produceChunk(&context);
678        if (! produced) {
679          break;
680        }
681        consumeChunk(&context);
682      }
683#endif
684
685      os << "end pol " << ipol << LogIO::POST ;
686
687    }
688    os << "end table " << ifile << LogIO::POST ;   
689  }
690  os << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ;
691  os << "getData: elapsed time is " << eGetData_-eToInt-eGetWeight << " sec." << LogIO::POST ;
692  os << "toPixel: elapsed time is " << eToPixel_ << " sec." << LogIO::POST ;
693  os << "ggridsd: elapsed time is " << eGGridSD_ << " sec." << LogIO::POST ;
694  os << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
695  os << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
696 
697  delete chanMap ;
698
699  // set data
700  setData( common.gdataArrC, common.gwgtArr ) ;
701
702}
703
704void STGrid::consumeChunkWithClipping(void *ctx) throw(PCException)
705{
706  STContextWithClipping &context = *(STContextWithClipping *)ctx;
707  STGChunk *chunk = NULL;
708  try {
709    context.queue.lock();
710    chunk = context.queue.get();
711    context.queue.unlock();
712  } catch (FullException &e) {
713    context.queue.unlock();
714    // TODO: log error
715    throw PCException();
716  }
717
718  double t0, t1 ;
719  // world -> pixel
720  Array<Double> xypos( context.self->dshape_ ) ;
721  t0 = mathutil::gettimeofday_sec() ;
722  context.self->toPixel( chunk->direction, xypos ) ;
723  t1 = mathutil::gettimeofday_sec() ;
724  context.self->eToPixel_ += t1-t0 ;
725   
726  // call ggridsd
727  Int nvispol = 1 ;
728  Int irow = -1 ;
729  t0 = mathutil::gettimeofday_sec() ;
730  context.self->call_ggridsd2( xypos,
731                chunk->spectra,
732                nvispol,
733                context.self->nchan_,
734                chunk->flagtra,
735                chunk->rflag,
736                chunk->weight,
737                chunk->nrow,
738                irow,
739                context.common.gdataArrC,
740                context.common.gwgtArr,
741                context.common.npoints,
742                context.common.clipMin,
743                context.common.clipWMin,
744                context.common.clipCMin,
745                context.common.clipMax,
746                context.common.clipWMax,
747                context.common.clipCMax,
748                context.common.gnx,
749                context.common.gny,
750                context.self->npol_,
751                context.self->nchan_,
752                context.self->convSupport_,
753                context.self->convSampling_,
754                context.common.convFunc,
755                context.common.chanMap,
756                (Int*)&context.pol ) ;
757  t1 = mathutil::gettimeofday_sec() ;
758  context.self->eGGridSD_ += t1-t0 ;
759 
760  delete chunk;
761}
762
763void STGrid::gridPerRowWithClipping()
764{
765  LogIO os( LogOrigin("STGrid", "gridPerRowWithClipping", WHERE) ) ;
766  double t0, t1 ;
767
768
769  // grid data
770  // Extend grid plane with convSupport_
771  //   Int gnx = nx_+convSupport_*2 ;
772  //   Int gny = ny_+convSupport_*2 ;
773  Int gnx = nx_;
774  Int gny = ny_;
775
776  IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
777  IPosition pshape( 3, gnx, gny, npol_ ) ;
778  // 2011/12/20 TN
779  // data_ and gwgtArr share storage
780  data_.resize( gshape ) ;
781  data_ = 0.0 ;
782  //STCommonData common = STCommonData(gshape, data_);
783  STCommonDataWithClipping common = STCommonDataWithClipping( gshape,
784                                                              pshape,
785                                                              data_ ) ;
786  common.gnx = gnx ;
787  common.gny = gny ;
788
789  // parameters for gridding
790  Int *chanMap = new Int[nchan_] ;
791  for ( Int i = 0 ; i < nchan_ ; i++ ) {
792    chanMap[i] = i ;
793  }
794  common.chanMap = chanMap;
795
796  // convolution kernel
797  t0 = mathutil::gettimeofday_sec() ;
798  setConvFunc( common.convFunc ) ;
799  t1 = mathutil::gettimeofday_sec() ;
800  os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
801
802  // for performance check
803  eGetData_ = 0.0 ;
804  eToPixel_ = 0.0 ;
805  eGGridSD_ = 0.0 ;
806  double eInitPol = 0.0 ;
807
808  //Broker broker = Broker(produceChunk, consumeChunk);
809  Broker broker = Broker(produceChunk, consumeChunkWithClipping);
810  for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) {
811    initTable( ifile ) ;
812
813    os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ;   
814    for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
815      t0 = mathutil::gettimeofday_sec() ;
816      initPol( ipol ) ; // set ptab_ and attach()
817      t1 = mathutil::gettimeofday_sec() ;
818      eInitPol += t1-t0 ;
819     
820      //STContext context(this, common, ipol);
821      STContextWithClipping context(this, common, ipol);
822     
823      os << "start pol " << ipol << LogIO::POST ;
824     
825      nprocessed_ = 0 ;
826#if 1
827      broker.runProducerAsMasterThread(&context, DO_AHEAD);
828#else
829      for (;;) {
830        bool produced = produceChunk(&context);
831        if (! produced) {
832          break;
833        }
834        //consumeChunk(&context);
835        consumeChunkWithClipping(&context);
836      }
837#endif
838
839      os << "end pol " << ipol << LogIO::POST ;
840
841    }
842    os << "end table " << ifile << LogIO::POST ;   
843  }
844  os << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ;
845  os << "getData: elapsed time is " << eGetData_-eToInt-eGetWeight << " sec." << LogIO::POST ;
846  os << "toPixel: elapsed time is " << eToPixel_ << " sec." << LogIO::POST ;
847  os << "ggridsd2: elapsed time is " << eGGridSD_ << " sec." << LogIO::POST ;
848  os << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
849  os << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
850 
851  delete chanMap ;
852
853  // clip min and max in each grid
854//   os << "BEFORE CLIPPING" << LogIO::POST ;
855//   os << "gdataArrC=" << common.gdataArrC << LogIO::POST ;
856//   os << "gwgtArr=" << common.gwgtArr << LogIO::POST ;
857  t0 = mathutil::gettimeofday_sec() ;
858  clipMinMax( common.gdataArrC,
859              common.gwgtArr,
860              common.npoints,
861              common.clipMin,
862              common.clipWMin,
863              common.clipCMin,
864              common.clipMax,
865              common.clipWMax,
866              common.clipCMax ) ;
867  t1 = mathutil::gettimeofday_sec() ;
868  os << "clipMinMax: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
869//   os << "AFTER CLIPPING" << LogIO::POST ;
870//   os << "gdataArrC=" << common.gdataArrC << LogIO::POST ;
871//   os << "gwgtArr=" << common.gwgtArr << LogIO::POST ;
872
873  // set data
874  setData( common.gdataArrC, common.gwgtArr ) ;
875
876}
877
878void STGrid::clipMinMax( Array<Complex> &grid,
879                         Array<Float> &weight,
880                         Array<Int> &npoints,
881                         Array<Complex> &clipmin,
882                         Array<Float> &clipwmin,
883                         Array<Float> &clipcmin,
884                         Array<Complex> &clipmax,
885                         Array<Float> &clipwmax,
886                         Array<Float> &clipcmax )
887{
888  //LogIO os( LogOrigin("STGrid","clipMinMax",WHERE) ) ;
889
890  // prepare pointers
891  Bool delG, delW, delNP, delCMin, delCWMin, delCCMin, delCMax, delCWMax, delCCMax ;
892  Complex *grid_p = grid.getStorage( delG ) ;
893  Float *wgt_p = weight.getStorage( delW ) ;
894  const Int *npts_p = npoints.getStorage( delNP ) ;
895  const Complex *cmin_p = clipmin.getStorage( delCMin ) ;
896  const Float *cwmin_p = clipwmin.getStorage( delCWMin ) ;
897  const Float *ccmin_p = clipcmin.getStorage( delCCMin ) ;
898  const Complex *cmax_p = clipmax.getStorage( delCMax ) ;
899  const Float *cwmax_p = clipwmax.getStorage( delCWMax ) ;
900  const Float *ccmax_p = clipcmax.getStorage( delCCMax ) ;
901
902  const IPosition &gshape = grid.shape() ;
903  long offset = gshape[0] * gshape[1] * gshape[2] ; // nx * ny * npol
904  Int nchan = gshape[3] ;
905  long origin = nchan * offset ;
906  for ( long i = 0 ; i < offset ; i++ ) {
907    if ( *npts_p > 2 ) {
908      for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
909        // clip minimum and maximum
910        *grid_p -= (*cmin_p)*(*cwmin_p)*(*ccmin_p)
911          + (*cmax_p)*(*cwmax_p)*(*ccmax_p) ;
912        *wgt_p -= (*cwmin_p)*(*ccmin_p)
913          + (*cwmax_p)*(*ccmax_p) ;
914       
915        grid_p += offset ;
916        wgt_p += offset ;
917        cmin_p += offset ;
918        cwmin_p += offset ;
919        ccmin_p += offset ;
920        cmax_p += offset ;
921        cwmax_p += offset ;
922        ccmax_p += offset ;
923      }
924      grid_p -= origin ;
925      wgt_p -= origin ;
926      cmin_p -= origin ;
927      cwmin_p -= origin ;
928      ccmin_p -= origin ;
929      cmax_p -= origin ;
930      cwmax_p -= origin ;
931      ccmax_p -= origin ;
932    }
933    grid_p++ ;
934    wgt_p++ ;
935    npts_p++ ;
936    cmin_p++ ;
937    cwmin_p++ ;
938    ccmin_p++ ;
939    cmax_p++ ;
940    cwmax_p++ ;
941    ccmax_p++ ;
942  }
943  grid_p -= offset ;
944  wgt_p -= offset ;
945  npts_p -= offset ;
946  cmin_p -= offset ;
947  cwmin_p -= offset ;
948  ccmin_p -= offset ;
949  cmax_p -= offset ;
950  cwmax_p -= offset ;
951  ccmax_p -= offset ; 
952
953  // finalization
954  grid.putStorage( grid_p, delG ) ;
955  weight.putStorage( wgt_p, delW ) ;
956  npoints.freeStorage( npts_p, delNP ) ;
957  clipmin.freeStorage( cmin_p, delCMin ) ;
958  clipwmin.freeStorage( cwmin_p, delCWMin ) ;
959  clipcmin.freeStorage( ccmin_p, delCCMin ) ;
960  clipmax.freeStorage( cmax_p, delCMax ) ;
961  clipwmax.freeStorage( cwmax_p, delCWMax ) ;
962  clipcmax.freeStorage( ccmax_p, delCCMax ) ;
963}
964
965void STGrid::gridPerPol()
966{
967  LogIO os( LogOrigin("STGrid", "gridPerPol", WHERE) ) ;
968  double t0, t1 ;
969
970  // convolution kernel
971  Vector<Float> convFunc ;
972  t0 = mathutil::gettimeofday_sec() ;
973  setConvFunc( convFunc ) ;
974  t1 = mathutil::gettimeofday_sec() ;
975  os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
976
977  // prepare grid data storage
978  Int gnx = nx_ ;
979  Int gny = ny_ ;
980//   // Extend grid plane with convSupport_
981//   Int gnx = nx_+convSupport_*2 ;
982//   Int gny = ny_+convSupport_*2 ;
983  IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
984  Array<Complex> gdataArrC( gshape, 0.0 ) ;
985  //Array<Float> gwgtArr( gshape, 0.0 ) ;
986  // 2011/12/20 TN
987  // data_ and weight array shares storage
988  data_.resize( gshape ) ;
989  data_ = 0.0 ;
990  Array<Float> gwgtArr( data_ ) ;
991
992  // maps
993  Int *chanMap = new Int[nchan_] ;
994  {
995    Int *work_p = chanMap ;
996    for ( Int i = 0 ; i < nchan_ ; i++ ) {
997      *work_p = i ;
998      work_p++ ;
999    }
1000  }
1001  Int polMap[1] ;
1002
1003  // some parameters for ggridsd
1004  Int nvispol = 1 ;
1005  Int irow = -1 ;
1006 
1007  // for performance check
1008  double eInitPol = 0.0 ;
1009  double eGetData = 0.0 ;
1010  double eToPixel = 0.0 ;
1011  double eGGridSD = 0.0 ;
1012  double eToInt = 0.0 ;
1013
1014  for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) {
1015    initTable( ifile ) ;
1016
1017    os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ;   
1018    // to read data from the table
1019    IPosition mshape( 2, nchan_, nrow_ ) ;
1020    IPosition dshape( 2, 2, nrow_ ) ;
1021    Array<Complex> spectra( mshape ) ;
1022    Array<Double> direction( dshape ) ;
1023    Array<Int> flagtra( mshape ) ;
1024    Array<Int> rflag( IPosition(1,nrow_) ) ;
1025    Array<Float> weight( mshape ) ;
1026    Array<Double> xypos( dshape ) ;
1027
1028    for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
1029      t0 = mathutil::gettimeofday_sec() ;
1030      initPol( ipol ) ;
1031      t1 = mathutil::gettimeofday_sec() ;
1032      eInitPol += t1-t0 ;
1033     
1034      os << "start pol " << ipol << LogIO::POST ;
1035     
1036      // retrieve data
1037      t0 = mathutil::gettimeofday_sec() ;
1038      getDataPerPol( spectra, direction, flagtra, rflag, weight ) ;
1039      t1 = mathutil::gettimeofday_sec() ;
1040      eGetData += t1-t0 ;
1041     
1042      // world -> pixel
1043      t0 = mathutil::gettimeofday_sec() ;
1044      toPixel( direction, xypos ) ; 
1045      t1 = mathutil::gettimeofday_sec() ;
1046      eToPixel += t1-t0 ;
1047     
1048      // call ggridsd
1049      polMap[0] = ipol ;
1050      t0 = mathutil::gettimeofday_sec() ;
1051      call_ggridsd( xypos,
1052                    spectra,
1053                    nvispol,
1054                    nchan_,
1055                    flagtra,
1056                    rflag,
1057                    weight,
1058                    nrow_,
1059                    irow,
1060                    gdataArrC,
1061                    gwgtArr,
1062                    gnx,
1063                    gny,
1064                    npol_,
1065                    nchan_,
1066                    convSupport_,
1067                    convSampling_,
1068                    convFunc,
1069                    chanMap,
1070                    polMap ) ;
1071      t1 = mathutil::gettimeofday_sec() ;
1072      eGGridSD += t1-t0 ;
1073
1074      os << "end pol " << ipol << LogIO::POST ;
1075
1076    }
1077    os << "end table " << ifile << LogIO::POST ;   
1078  }
1079  os << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ;
1080  os << "getData: elapsed time is " << eGetData-eToInt-eGetWeight << " sec." << LogIO::POST ;
1081  os << "toPixel: elapsed time is " << eToPixel << " sec." << LogIO::POST ;
1082  os << "ggridsd: elapsed time is " << eGGridSD << " sec." << LogIO::POST ;
1083  os << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
1084  os << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
1085
1086  // delete maps
1087  delete chanMap ;
1088
1089  setData( gdataArrC, gwgtArr ) ;
1090//   os << "gdataArr = " << gdataArr << LogIO::POST ;
1091//   os << "gwgtArr = " << gwgtArr << LogIO::POST ;
1092//   os << "data_ " << data_ << LogIO::POST ;
1093}
1094
1095void STGrid::initPol( Int ipol )
1096{
1097  LogIO os( LogOrigin("STGrid","initPol",WHERE) ) ;
1098  if ( npolOrg_ == 1 ) {
1099    os << "single polarization data." << LogIO::POST ;
1100    ptab_ = tab_ ;
1101  }
1102  else
1103    ptab_ = tab_( tab_.col("POLNO") == (uInt)ipol ) ;
1104
1105  attach( ptab_ ) ;
1106}
1107
1108void STGrid::initTable( uInt idx )
1109{
1110  tab_ = tableList_[idx] ;
1111  nrow_ = rows_[idx] ;
1112}
1113
1114void STGrid::setData( Array<Complex> &gdata,
1115                      Array<Float> &gwgt )
1116{
1117  // 2011/12/20 TN
1118  // gwgt and data_ share storage
1119  LogIO os( LogOrigin("STGrid","setData",WHERE) ) ;
1120  double t0, t1 ;
1121  t0 = mathutil::gettimeofday_sec() ;
1122  uInt len = data_.nelements() ;
1123  const Complex *w1_p ;
1124  Float *w2_p ;
1125  Bool b1, b2 ;
1126  const Complex *gdata_p = gdata.getStorage( b1 ) ;
1127  Float *gwgt_p = data_.getStorage( b2 ) ;
1128  w1_p = gdata_p ;
1129  w2_p = gwgt_p ;
1130  for ( uInt i = 0 ; i < len ; i++ ) {
1131    if ( *w2_p > 0.0 ) *w2_p = (*w1_p).real() / *w2_p ;
1132    w1_p++ ;
1133    w2_p++ ;
1134  }
1135  gdata.freeStorage( gdata_p, b1 ) ;
1136  data_.putStorage( gwgt_p, b2 ) ;
1137  t1 = mathutil::gettimeofday_sec() ;
1138  os << "setData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
1139}
1140
1141void STGrid::setupGrid()
1142{
1143  Double xmin,xmax,ymin,ymax ;
1144  mapExtent( xmin, xmax, ymin, ymax ) ;
1145 
1146  setupGrid( nxUI_, nyUI_, cellxUI_, cellyUI_,
1147             xmin, xmax, ymin, ymax, centerUI_ ) ;
1148}
1149
1150void STGrid::setupGrid( Int &nx,
1151                        Int &ny,
1152                        String &cellx,
1153                        String &celly,
1154                        Double &xmin,
1155                        Double &xmax,
1156                        Double &ymin,
1157                        Double &ymax,
1158                        String &center )
1159{
1160  LogIO os( LogOrigin("STGrid","setupGrid",WHERE) ) ;
1161  //cout << "nx=" << nx << ", ny=" << ny << endl ;
1162
1163  // center position
1164  if ( center.size() == 0 ) {
1165    center_(0) = 0.5 * ( xmin + xmax ) ;
1166    center_(1) = 0.5 * ( ymin + ymax ) ;
1167  }
1168  else {
1169    String::size_type pos0 = center.find( " " ) ;
1170    if ( pos0 == String::npos ) {
1171      throw AipsError( "bad string format in parameter center" ) ;
1172    }
1173    String::size_type pos1 = center.find( " ", pos0+1 ) ;
1174    String typestr, xstr, ystr ;
1175    if ( pos1 != String::npos ) {
1176      typestr = center.substr( 0, pos0 ) ;
1177      xstr = center.substr( pos0+1, pos1-pos0 ) ;
1178      ystr = center.substr( pos1+1 ) ;
1179      // todo: convert to J2000 (or direction ref for DIRECTION column)
1180    }
1181    else {
1182      typestr = "J2000" ;
1183      xstr = center.substr( 0, pos0 ) ;
1184      ystr = center.substr( pos0+1 ) ;
1185    }
1186    QuantumHolder qh ;
1187    String err ;
1188    qh.fromString( err, xstr ) ;
1189    Quantum<Double> xcen = qh.asQuantumDouble() ;
1190    qh.fromString( err, ystr ) ;
1191    Quantum<Double> ycen = qh.asQuantumDouble() ;
1192    center_(0) = xcen.getValue( "rad" ) ;
1193    center_(1) = ycen.getValue( "rad" ) ;
1194  }
1195
1196
1197  nx_ = nx ;
1198  ny_ = ny ;
1199  if ( nx < 0 && ny > 0 ) {
1200    nx_ = ny ;
1201    ny_ = ny ;
1202  }
1203  if ( ny < 0 && nx > 0 ) {
1204    nx_ = nx ;
1205    ny_ = nx ;
1206  }
1207
1208  //Double wx = xmax - xmin ;
1209  //Double wy = ymax - ymin ;
1210  Double wx = max( abs(xmax-center_(0)), abs(xmin-center_(0)) ) * 2 ;
1211  Double wy = max( abs(ymax-center_(1)), abs(ymin-center_(1)) ) * 2 ;
1212  // take 10% margin
1213  wx *= 1.10 ;
1214  wy *= 1.10 ;
1215
1216  Quantum<Double> qcellx ;
1217  Quantum<Double> qcelly ;
1218  //cout << "nx_ = " << nx_ << ",  ny_ = " << ny_ << endl ;
1219  if ( cellx.size() != 0 && celly.size() != 0 ) {
1220    readQuantity( qcellx, cellx ) ;
1221    readQuantity( qcelly, celly ) ;
1222  }
1223  else if ( celly.size() != 0 ) {
1224    os << "Using celly to x-axis..." << LogIO::POST ;
1225    readQuantity( qcelly, celly ) ;
1226    qcellx = qcelly ;
1227  }
1228  else if ( cellx.size() != 0 ) {
1229    os << "Using cellx to y-axis..." << LogIO::POST ;
1230    readQuantity( qcellx, cellx ) ;
1231    qcelly = qcellx ;
1232  }
1233  else {
1234    if ( nx_ < 0 ) {
1235      os << "No user preference in grid setting. Using default..." << LogIO::POST ;
1236      readQuantity( qcellx, "1.0arcmin" ) ;
1237      qcelly = qcellx ;
1238    }
1239    else {
1240      if ( wx == 0.0 ) {
1241        os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
1242        wx = 0.00290888 ;
1243      }
1244      if ( wy == 0.0 ) {
1245        os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
1246        wy = 0.00290888 ;
1247      }
1248      qcellx = Quantum<Double>( wx/nx_, "rad" ) ;
1249      qcelly = Quantum<Double>( wy/ny_, "rad" ) ;
1250    }
1251  }
1252  cellx_ = qcellx.getValue( "rad" ) ;
1253  celly_ = qcelly.getValue( "rad" ) ;
1254  if ( nx_ < 0 ) {
1255    if ( wx == 0.0 ) {
1256      os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
1257      wx = 0.00290888 ;
1258    }
1259    if ( wy == 0.0 ) {
1260      os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
1261      wy = 0.00290888 ;
1262    }
1263    nx_ = Int( ceil( wx/cellx_ ) ) ;
1264    ny_ = Int( ceil( wy/celly_ ) ) ;
1265  }
1266}
1267
1268void STGrid::mapExtent( Double &xmin, Double &xmax,
1269                        Double &ymin, Double &ymax )
1270{
1271  //LogIO os( LogOrigin("STGrid","mapExtent",WHERE) ) ;
1272  directionCol_.attach( tableList_[0], "DIRECTION" ) ;
1273  Matrix<Double> direction = directionCol_.getColumn() ;
1274  //os << "dirCol.nrow() = " << dirCol.nrow() << LogIO::POST ;
1275  minMax( xmin, xmax, direction.row( 0 ) ) ;
1276  minMax( ymin, ymax, direction.row( 1 ) ) ;
1277  Double amin, amax, bmin, bmax ;
1278  for ( uInt i = 1 ; i < nfile_ ; i++ ) {
1279    directionCol_.attach( tableList_[i], "DIRECTION" ) ;
1280    direction.assign( directionCol_.getColumn() ) ;
1281    //os << "dirCol.nrow() = " << dirCol.nrow() << LogIO::POST ;
1282    minMax( amin, amax, direction.row( 0 ) ) ;
1283    minMax( bmin, bmax, direction.row( 1 ) ) ;
1284    xmin = min( xmin, amin ) ;
1285    xmax = max( xmax, amax ) ;
1286    ymin = min( ymin, bmin ) ;
1287    ymax = max( ymax, bmax ) ;
1288  }
1289  //os << "(xmin,xmax)=(" << xmin << "," << xmax << ")" << LogIO::POST ;
1290  //os << "(ymin,ymax)=(" << ymin << "," << ymax << ")" << LogIO::POST ;
1291}
1292
1293void STGrid::selectData()
1294{
1295  LogIO os( LogOrigin("STGrid","selectData",WHERE) ) ;   
1296  Int ifno = ifno_ ;
1297  tableList_.resize( nfile_ ) ;
1298  if ( ifno == -1 ) {
1299    Table taborg( infileList_[0] ) ;
1300    ROScalarColumn<uInt> ifnoCol( taborg, "IFNO" ) ;
1301    ifno = ifnoCol( 0 ) ;
1302    os << LogIO::WARN
1303       << "IFNO is not given. Using default IFNO: " << ifno << LogIO::POST ;
1304  }
1305  for ( uInt i = 0 ; i < nfile_ ; i++ ) {
1306    Table taborg( infileList_[i] ) ;
1307    TableExprNode node ;
1308    if ( isMultiIF( taborg ) ) {
1309      os << "apply selection on IFNO" << LogIO::POST ;
1310      node = taborg.col("IFNO") == ifno ;
1311    }
1312    if ( scanlist_.size() > 0 ) {
1313      os << "apply selection on SCANNO" << LogIO::POST ;
1314      node = node && taborg.col("SCANNO").in( scanlist_ ) ;
1315    }
1316    if ( node.isNull() ) {
1317      tableList_[i] = taborg ;
1318    }
1319    else {
1320      tableList_[i] = taborg( node ) ;
1321    }
1322    os << "tableList_[" << i << "].nrow()=" << tableList_[i].nrow() << LogIO::POST ;
1323    if ( tableList_[i].nrow() == 0 ) {
1324      os << LogIO::SEVERE
1325         << "No corresponding rows for given selection: IFNO " << ifno ;
1326      if ( scanlist_.size() > 0 )
1327        os << " SCANNO " << scanlist_ ;
1328      os << LogIO::EXCEPTION ;
1329    }
1330  }
1331}
1332
1333Bool STGrid::isMultiIF( Table &tab )
1334{
1335  ROScalarColumn<uInt> ifnoCol( tab, "IFNO" ) ;
1336  Vector<uInt> ifnos = ifnoCol.getColumn() ;
1337  return anyNE( ifnos, ifnos[0] ) ;
1338}
1339
1340void STGrid::attach( Table &tab )
1341{
1342  // attach to table
1343  spectraCol_.attach( tab, "SPECTRA" ) ;
1344  flagtraCol_.attach( tab, "FLAGTRA" ) ;
1345  directionCol_.attach( tab, "DIRECTION" ) ;
1346  flagRowCol_.attach( tab, "FLAGROW" ) ;
1347  tsysCol_.attach( tab, "TSYS" ) ;
1348  intervalCol_.attach( tab, "INTERVAL" ) ;
1349}
1350
1351void STGrid::getDataPerPol( Array<Float> &spectra,
1352                            Array<Double> &direction,
1353                            Array<uChar> &flagtra,
1354                            Array<uInt> &rflag,
1355                            Array<Float> &weight )
1356{
1357  LogIO os( LogOrigin("STGrid","getDataPerPol",WHERE) ) ;
1358
1359  // 2011/12/22 TN
1360  // weight and tsys shares its storage
1361  Array<Float> tsys( weight ) ;
1362  Array<Double> tint( rflag.shape() ) ;
1363
1364  Vector<uInt> rflagVec( rflag ) ;
1365  Vector<Double> tintVec( tint ) ;
1366
1367  spectraCol_.getColumn( spectra ) ;
1368  flagtraCol_.getColumn( flagtra ) ;
1369  flagRowCol_.getColumn( rflagVec ) ;
1370  directionCol_.getColumn( direction ) ;
1371  intervalCol_.getColumn( tintVec ) ;
1372  Vector<Float> tsysTemp = tsysCol_( 0 ) ;
1373  if ( tsysTemp.nelements() == (uInt)nchan_ ) {
1374    tsysCol_.getColumn( tsys ) ;
1375  }
1376  else {
1377    tsys = tsysTemp[0] ;
1378  }
1379 
1380  double t0,t1 ;
1381  t0 = mathutil::gettimeofday_sec() ;
1382  getWeight( weight, tsys, tint ) ;
1383  t1 = mathutil::gettimeofday_sec() ;
1384  eGetWeight += t1-t0 ;
1385}
1386
1387void STGrid::getDataPerPol( Array<Complex> &spectra,
1388                            Array<Double> &direction,
1389                            Array<Int> &flagtra,
1390                            Array<Int> &rflag,
1391                            Array<Float> &weight )
1392{
1393  LogIO os( LogOrigin("STGrid","getDataPerPol",WHERE) ) ;
1394  double t0, t1 ;
1395
1396  Array<uChar> flagUC( flagtra.shape() ) ;
1397  Array<uInt> rflagUI( rflag.shape() ) ;
1398  Array<Float> spectraF( spectra.shape() ) ;
1399  getDataPerPol( spectraF, direction, flagUC, rflagUI, weight ) ;
1400
1401  convertArray( spectra, spectraF ) ;
1402  t0 = mathutil::gettimeofday_sec() ;
1403  toInt( flagUC, flagtra ) ;
1404  toInt( rflagUI, rflag ) ;
1405  t1 = mathutil::gettimeofday_sec() ;
1406  eToInt += t1-t0 ;
1407  //os << "toInt: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
1408}
1409
1410Int STGrid::getDataChunk(
1411                         IPosition const &wshape,
1412                         IPosition const &vshape,
1413                         IPosition const &dshape,
1414                         Array<Complex> &spectra,
1415                         Array<Double> &direction,
1416                         Array<Int> &flagtra,
1417                         Array<Int> &rflag,
1418                         Array<Float> &weight )
1419{
1420  LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
1421
1422  Array<Float> spectraF_(wshape);
1423  Array<uChar> flagtraUC_(wshape);
1424  Array<uInt> rflagUI_(vshape);
1425  Int nrow = getDataChunk( spectraF_, direction, flagtraUC_, rflagUI_, weight ) ;
1426  if ( nrow < nchunk_ ) {
1427    spectra.resize( spectraF_.shape() ) ;
1428    flagtra.resize( flagtraUC_.shape() ) ;
1429    rflag.resize( rflagUI_.shape() ) ;
1430  }
1431  double t0, t1 ;
1432  t0 = mathutil::gettimeofday_sec() ;
1433  convertArray( spectra, spectraF_ ) ;
1434  toInt( flagtraUC_, flagtra ) ;
1435  toInt( rflagUI_, rflag ) ;
1436  t1 = mathutil::gettimeofday_sec() ;
1437  eToInt = t1 - t0 ;
1438 
1439  return nrow ;
1440}
1441
1442#if 0
1443Int STGrid::getDataChunk( Array<Complex> &spectra,
1444                          Array<Double> &direction,
1445                          Array<Int> &flagtra,
1446                          Array<Int> &rflag,
1447                          Array<Float> &weight )
1448{
1449  LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
1450  Int nrow = getDataChunk( spectraF_, direction, flagtraUC_, rflagUI_, weight ) ;
1451  if ( nrow < nchunk_ ) {
1452    spectra.resize( spectraF_.shape() ) ;
1453    flagtra.resize( flagtraUC_.shape() ) ;
1454    rflag.resize( rflagUI_.shape() ) ;
1455  }
1456  double t0, t1 ;
1457  t0 = mathutil::gettimeofday_sec() ;
1458  convertArray( spectra, spectraF_ ) ;
1459  toInt( flagtraUC_, flagtra ) ;
1460  toInt( rflagUI_, rflag ) ;
1461  t1 = mathutil::gettimeofday_sec() ;
1462  eToInt = t1 - t0 ;
1463 
1464  return nrow ;
1465}
1466#endif
1467
1468Int STGrid::getDataChunk( Array<Float> &spectra,
1469                          Array<Double> &direction,
1470                          Array<uChar> &flagtra,
1471                          Array<uInt> &rflag,
1472                          Array<Float> &weight )
1473{
1474  LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
1475  Int nrow = spectra.shape()[1] ;
1476  Int remainingRow = nrow_ - nprocessed_ ;
1477  if ( remainingRow < nrow ) {
1478    nrow = remainingRow ;
1479    IPosition mshape( 2, nchan_, nrow ) ;
1480    IPosition vshape( 1, nrow ) ;
1481    spectra.resize( mshape ) ;
1482    flagtra.resize( mshape ) ;
1483    direction.resize( IPosition(2,2,nrow) ) ;
1484    rflag.resize( vshape ) ;
1485    weight.resize( mshape ) ;
1486  }
1487  // 2011/12/22 TN
1488  // tsys shares its storage with weight
1489  Array<Float> tsys( weight ) ;
1490  Array<Double> tint( rflag.shape() ) ;
1491
1492  Vector<uInt> rflagVec( rflag ) ;
1493  Vector<Double> tintVec( tint ) ;
1494
1495  RefRows rows( nprocessed_, nprocessed_+nrow-1, 1 ) ;
1496  //os<<LogIO::DEBUGGING<<"nprocessed_="<<nprocessed_<<": rows.nrows()="<<rows.nrows()<<LogIO::POST ;
1497  spectraCol_.getColumnCells( rows, spectra ) ;
1498  flagtraCol_.getColumnCells( rows, flagtra ) ;
1499  directionCol_.getColumnCells( rows, direction ) ;
1500  flagRowCol_.getColumnCells( rows, rflagVec ) ;
1501  intervalCol_.getColumnCells( rows, tintVec ) ;
1502  Vector<Float> tsysTemp = tsysCol_( nprocessed_ ) ;
1503  if ( tsysTemp.nelements() == (uInt)nchan_ )
1504    tsysCol_.getColumnCells( rows, tsys ) ;
1505  else
1506    tsys = tsysTemp[0] ;
1507
1508  double t0,t1 ;
1509  t0 = mathutil::gettimeofday_sec() ;
1510  getWeight( weight, tsys, tint ) ;
1511  t1 = mathutil::gettimeofday_sec() ;
1512  eGetWeight += t1-t0 ;
1513
1514  nprocessed_ += nrow ;
1515 
1516  return nrow ;
1517}
1518
1519void STGrid::setupArray()
1520{
1521  LogIO os( LogOrigin("STGrid","setupArray",WHERE) ) ;
1522  ROScalarColumn<uInt> polnoCol( tableList_[0], "POLNO" ) ;
1523  Vector<uInt> pols = polnoCol.getColumn() ;
1524  //os << pols << LogIO::POST ;
1525  Vector<uInt> pollistOrg ;
1526  npolOrg_ = 0 ;
1527  uInt polno ;
1528  for ( uInt i = 0 ; i < polnoCol.nrow() ; i++ ) {
1529    //polno = polnoCol( i ) ;
1530    polno = pols( i ) ;
1531    if ( allNE( pollistOrg, polno ) ) {
1532      pollistOrg.resize( npolOrg_+1, True ) ;
1533      pollistOrg[npolOrg_] = polno ;
1534      npolOrg_++ ;
1535    }
1536  }
1537  if ( pollist_.size() == 0 )
1538    pollist_ = pollistOrg ;
1539  else {
1540    Vector<uInt> newlist ;
1541    uInt newsize = 0 ;
1542    for ( uInt i = 0 ; i < pollist_.size() ; i++ ) {
1543      if ( anyEQ( pollistOrg, pollist_[i] ) ) {
1544        newlist.resize( newsize+1, True ) ;
1545        newlist[newsize] = pollist_[i] ;
1546        newsize++ ;
1547      }
1548    }
1549    pollist_.assign( newlist ) ;
1550  }
1551  npol_ = pollist_.size() ;
1552  if ( npol_ == 0 ) {
1553    os << LogIO::SEVERE << "Empty pollist" << LogIO::EXCEPTION ;
1554  }
1555  rows_.resize( nfile_ ) ;
1556  for ( uInt i = 0 ; i < nfile_ ; i++ ) {
1557    rows_[i] = tableList_[i].nrow() / npolOrg_ ;
1558    if ( nrow_ < rows_[i] )
1559      nrow_ = rows_[i] ;
1560  }
1561  flagtraCol_.attach( tableList_[0], "FLAGTRA" ) ;
1562  nchan_ = flagtraCol_( 0 ).nelements() ;
1563//   os << "npol_ = " << npol_ << "(" << pollist_ << ")" << endl
1564//      << "nchan_ = " << nchan_ << endl
1565//      << "nrow_ = " << nrow_ << LogIO::POST ;
1566}
1567
1568void STGrid::getWeight( Array<Float> &w,
1569                              Array<Float> &tsys,
1570                              Array<Double> &tint )
1571{
1572  LogIO os( LogOrigin("STGrid","getWeight",WHERE) ) ;
1573
1574  // 2011/12/22 TN
1575  // w (weight) and tsys share storage
1576  IPosition refShape = tsys.shape() ;
1577  Int nchan = refShape[0] ;
1578  Int nrow = refShape[1] ;
1579//   os << "nchan=" << nchan << ", nrow=" << nrow << LogIO::POST ;
1580//   os << "w.shape()=" << w.shape() << endl
1581//      << "tsys.shape()=" << tsys.shape() << endl
1582//      << "tint.shape()=" << tint.shape() << LogIO::POST ;
1583
1584  // set weight
1585  if ( wtype_.compare( "UNIFORM" ) == 0 ) {
1586    w = 1.0 ;
1587  }
1588  else if ( wtype_.compare( "TINT" ) == 0 ) {
1589    Bool b0, b1 ;
1590    Float *w_p = w.getStorage( b0 ) ;
1591    Float *w0_p = w_p ;
1592    const Double *ti_p = tint.getStorage( b1 ) ;
1593    const Double *w1_p = ti_p ;
1594    for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1595      for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1596        *w0_p = *w1_p ;
1597        w0_p++ ;
1598      }
1599      w1_p++ ;
1600    }
1601    w.putStorage( w_p, b0 ) ;
1602    tint.freeStorage( ti_p, b1 ) ;
1603  }
1604  else if ( wtype_.compare( "TSYS" ) == 0 ) {
1605    Bool b0 ;
1606    Float *w_p = w.getStorage( b0 ) ;
1607    Float *w0_p = w_p ;
1608    for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1609      for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1610        Float temp = *w0_p ;
1611        *w0_p = 1.0 / ( temp * temp ) ;
1612        w0_p++ ;
1613      }
1614    }
1615    w.putStorage( w_p, b0 ) ;
1616  }
1617  else if ( wtype_.compare( "TINTSYS" ) == 0 ) {
1618    Bool b0, b1 ;
1619    Float *w_p = w.getStorage( b0 ) ;
1620    Float *w0_p = w_p ;
1621    const Double *ti_p = tint.getStorage( b1 ) ;
1622    const Double *w1_p = ti_p ;
1623    for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1624      Float interval = *w1_p ;
1625      for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1626        Float temp = *w0_p ;
1627        *w0_p = interval / ( temp * temp ) ;
1628        w0_p++ ;
1629      }
1630      w1_p++ ;
1631    }
1632    w.putStorage( w_p, b0 ) ;
1633    tint.freeStorage( ti_p, b1 ) ;
1634  }
1635  else {
1636    //LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
1637    //os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
1638    w = 1.0 ;
1639  }
1640}
1641
1642void STGrid::toInt( Array<uChar> &u, Array<Int> &v )
1643{
1644  uInt len = u.nelements() ;
1645  Int *int_p = new Int[len] ;
1646  Bool deleteIt ;
1647  const uChar *data_p = u.getStorage( deleteIt ) ;
1648  Int *i_p = int_p ;
1649  const uChar *u_p = data_p ;
1650  for ( uInt i = 0 ; i < len ; i++ ) {
1651    *i_p = ( *u_p == 0 ) ? 0 : 1 ;
1652    i_p++ ;
1653    u_p++ ;
1654  }
1655  u.freeStorage( data_p, deleteIt ) ;
1656  v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
1657}
1658
1659void STGrid::toInt( Array<uInt> &u, Array<Int> &v )
1660{
1661  uInt len = u.nelements() ;
1662  Int *int_p = new Int[len] ;
1663  Bool deleteIt ;
1664  const uInt *data_p = u.getStorage( deleteIt ) ;
1665  Int *i_p = int_p ;
1666  const uInt *u_p = data_p ;
1667  for ( uInt i = 0 ; i < len ; i++ ) {
1668    *i_p = ( *u_p == 0 ) ? 0 : 1 ;
1669    i_p++ ;
1670    u_p++ ;
1671  }
1672  u.freeStorage( data_p, deleteIt ) ;
1673  v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
1674}
1675
1676void STGrid::toPixel( Array<Double> &world, Array<Double> &pixel )
1677{
1678  // gridding will be done on (nx_+2*convSupport_) x (ny_+2*convSupport_)
1679  // grid plane to avoid unexpected behavior on grid edge
1680  Block<Double> pixc( 2 ) ;
1681  pixc[0] = Double( nx_-1 ) * 0.5 ;
1682  pixc[1] = Double( ny_-1 ) * 0.5 ;
1683//   pixc[0] = Double( nx_+2*convSupport_-1 ) * 0.5 ;
1684//   pixc[1] = Double( ny_+2*convSupport_-1 ) * 0.5 ;
1685  uInt nrow = world.shape()[1] ;
1686  Bool bw, bp ;
1687  const Double *w_p = world.getStorage( bw ) ;
1688  Double *p_p = pixel.getStorage( bp ) ;
1689  const Double *ww_p = w_p ;
1690  Double *wp_p = p_p ;
1691  for ( uInt i = 0 ; i < nrow ; i++ ) {
1692    *wp_p = pixc[0] + ( *ww_p - center_[0] ) / cellx_ ;
1693    wp_p++ ;
1694    ww_p++ ;
1695    *wp_p = pixc[1] + ( *ww_p - center_[1] ) / celly_ ;
1696    wp_p++ ;
1697    ww_p++ ;
1698  }
1699  world.freeStorage( w_p, bw ) ;
1700  pixel.putStorage( p_p, bp ) ; 
1701//   String gridfile = "grid."+convType_+"."+String::toString(convSupport_)+".dat" ;
1702//   ofstream ofs( gridfile.c_str(), ios::out ) ;
1703//   ofs << "center " << center_(0) << " " << pixc(0)
1704//       << " " << center_(1) << " " << pixc(1) << endl ;
1705//   for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
1706//     ofs << irow ;
1707//     for ( uInt i = 0 ; i < 2 ; i++ ) {
1708//       ofs << " " << world(i, irow) << " " << pixel(i, irow) ;
1709//     }
1710//     ofs << endl ;
1711//   }
1712//   ofs.close() ;
1713}
1714
1715void STGrid::boxFunc( Vector<Float> &convFunc, Int &convSize )
1716{
1717  convFunc = 0.0 ;
1718  for ( Int i = 0 ; i < convSize/2 ; i++ )
1719    convFunc(i) = 1.0 ;
1720}
1721
1722#define NEED_UNDERSCORES
1723#if defined(NEED_UNDERSCORES)
1724#define grdsf grdsf_
1725#endif
1726extern "C" {
1727   void grdsf(Double*, Double*);
1728}
1729void STGrid::spheroidalFunc( Vector<Float> &convFunc )
1730{
1731  convFunc = 0.0 ;
1732  for ( Int i = 0 ; i < convSampling_*convSupport_ ; i++ ) {
1733    Double nu = Double(i) / Double(convSupport_*convSampling_) ;
1734    Double val ;
1735    grdsf( &nu, &val ) ;
1736    convFunc(i) = ( 1.0 - nu * nu ) * val ;
1737  }
1738}
1739
1740void STGrid::gaussFunc( Vector<Float> &convFunc )
1741{
1742  convFunc = 0.0 ;
1743  // HWHM of the Gaussian is convSupport_ / 4
1744  // To take into account Gaussian tail, kernel cutoff is set to 4 * HWHM
1745  Int len = convSampling_ * convSupport_ ;
1746  Double hwhm = len * 0.25 ;
1747  for ( Int i = 0 ; i < len ; i++ ) {
1748    Double val = Double(i) / hwhm ;
1749    convFunc(i) = exp( -log(2)*val*val ) ;
1750  }
1751}
1752
1753void STGrid::pbFunc( Vector<Float> &convFunc )
1754{
1755  convFunc = 0.0 ;
1756}
1757
1758void STGrid::setConvFunc( Vector<Float> &convFunc )
1759{
1760  convSupport_ = userSupport_ ;
1761  if ( convType_ == "BOX" ) {
1762    if ( convSupport_ < 0 )
1763      convSupport_ = 0 ;
1764    Int convSize = convSampling_ * ( 2 * convSupport_ + 2 )  ;
1765    convFunc.resize( convSize ) ;
1766    boxFunc( convFunc, convSize ) ;
1767  }
1768  else if ( convType_ == "SF" ) {
1769    if ( convSupport_ < 0 )
1770      convSupport_ = 3 ;
1771    Int convSize = convSampling_ * ( 2 * convSupport_ + 2 )  ;
1772    convFunc.resize( convSize ) ;
1773    spheroidalFunc( convFunc ) ;
1774  }
1775  else if ( convType_ == "GAUSS" ) {
1776    // to take into account Gaussian tail
1777    if ( convSupport_ < 0 )
1778      convSupport_ = 12 ; // 3 * 4
1779    else {
1780      convSupport_ = userSupport_ * 4 ;
1781    }
1782    Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
1783    convFunc.resize( convSize ) ;
1784    gaussFunc( convFunc ) ;
1785  }
1786  else if ( convType_ == "PB" ) {
1787    if ( convSupport_ < 0 )
1788      convSupport_ = 0 ;
1789    pbFunc( convFunc ) ;
1790  }
1791  else {
1792    throw AipsError( "Unsupported convolution function" ) ;
1793  }
1794}
1795
1796string STGrid::saveData( string outfile )
1797{
1798  LogIO os( LogOrigin("STGrid", "saveData", WHERE) ) ;
1799  double t0, t1 ;
1800  t0 = mathutil::gettimeofday_sec() ;
1801
1802  //Int polno = 0 ;
1803  String outfile_ ;
1804  if ( outfile.size() == 0 ) {
1805    if ( infileList_[0].lastchar() == '/' ) {
1806      outfile_ = infileList_[0].substr( 0, infileList_[0].size()-1 ) ;
1807    }
1808    else {
1809      outfile_ = infileList_[0] ;
1810    }
1811    outfile_ += ".grid" ;
1812  }
1813  else {
1814    outfile_ = outfile ;
1815  }
1816  Table tab ;
1817  prepareTable( tab, outfile_ ) ;
1818  IPosition dshape = data_.shape() ;
1819  Int nrow = nx_ * ny_ * npol_ ;
1820  tab.rwKeywordSet().define( "nPol", npol_ ) ;
1821  tab.addRow( nrow ) ;
1822  Vector<Double> cpix( 2 ) ;
1823  cpix(0) = Double( nx_ - 1 ) * 0.5 ;
1824  cpix(1) = Double( ny_ - 1 ) * 0.5 ;
1825  Vector<Double> dir( 2 ) ;
1826  ArrayColumn<Double> directionCol( tab, "DIRECTION" ) ;
1827  ArrayColumn<Float> spectraCol( tab, "SPECTRA" ) ;
1828  ScalarColumn<uInt> polnoCol( tab, "POLNO" ) ;
1829  Int irow = 0 ;
1830  Vector<Float> sp( nchan_ ) ;
1831  Bool bsp, bdata ;
1832  const Float *data_p = data_.getStorage( bdata ) ;
1833  Float *wsp_p, *sp_p ;
1834  const Float *wdata_p = data_p ;
1835  long step = nx_ * ny_ * npol_ ;
1836  long offset ;
1837  for ( Int iy = 0 ; iy < ny_ ; iy++ ) {
1838    dir(1) = center_(1) - ( cpix(1) - (Double)iy ) * celly_ ;
1839    for ( Int ix = 0 ; ix < nx_ ; ix++ ) {
1840      dir(0) = center_(0) - ( cpix(0) - (Double)ix ) * cellx_ ;
1841      for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
1842        offset = ix + iy * nx_ + ipol * nx_ * ny_ ;
1843        //os << "offset = " << offset << LogIO::POST ;
1844        sp_p = sp.getStorage( bsp ) ;
1845        wsp_p = sp_p ;
1846        wdata_p = data_p + offset ;
1847        for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
1848          *wsp_p = *wdata_p ;
1849          wsp_p++ ;
1850          wdata_p += step ;
1851        }
1852        sp.putStorage( sp_p, bsp ) ;
1853        spectraCol.put( irow, sp ) ;
1854        directionCol.put( irow, dir ) ;
1855        polnoCol.put( irow, pollist_[ipol] ) ;
1856        irow++ ;
1857      }
1858    }
1859  }
1860  data_.freeStorage( data_p, bdata ) ;
1861
1862  t1 = mathutil::gettimeofday_sec() ;
1863  os << "saveData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
1864
1865  return outfile_ ;
1866}
1867
1868void STGrid::prepareTable( Table &tab, String &name )
1869{
1870  Table t( infileList_[0], Table::Old ) ;
1871  t.deepCopy( name, Table::New, False, t.endianFormat(), True ) ;
1872  tab = Table( name, Table::Update ) ;
1873}
1874
1875}
Note: See TracBrowser for help on using the repository browser.