source: trunk/src/STGrid.cpp @ 3106

Last change on this file since 3106 was 3106, checked in by Takeshi Nakazato, 8 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes/No?

Interface Changes: Yes/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...


Check-in asap modifications from Jim regarding casacore namespace conversion.
File size: 63.6 KB
Line 
1//
2// C++ Implementation: STGrid
3//
4// Description:
5//
6//
7// Author: Takeshi Nakazato <takeshi.nakazato@nao.ac.jp>, (C) 2011
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12
13// ASAP
14// STGrid must be included first to avoid compiler warning related
15// with _XOPEN_SOURCE
16#include "STGrid.h"
17#include "MathUtils.h"
18#include <atnf/PKSIO/SrcType.h>
19
20#include <casa/BasicSL/String.h>
21#include <casa/Arrays/Vector.h>
22#include <casa/Arrays/ArrayMath.h>
23#include <casa/Quanta/Quantum.h>
24#include <casa/Quanta/QuantumHolder.h>
25#include <casa/Utilities/CountedPtr.h>
26#include <casa/Logging/LogIO.h>
27
28#include <coordinates/Coordinates/DirectionCoordinate.h>
29
30#include <tables/Tables/Table.h>
31#include <tables/Tables/TableRecord.h>
32#include <tables/Tables/TableRow.h>
33#include <tables/TaQL/ExprNode.h>
34#include <tables/Tables/ScalarColumn.h>
35#include <tables/Tables/ArrayColumn.h>
36#include <tables/Tables/TableCopy.h>
37
38#include <measures/Measures/MDirection.h>
39
40using namespace std ;
41using namespace concurrent ;
42using namespace casacore;
43using namespace asap ;
44
45namespace asap {
46
47// for performance check
48double eToInt = 0.0 ;
49double eGetWeight = 0.0 ;
50
51// constructor
52STGrid::STGrid()
53  : vshape_( 1 ), wshape_( 2 ), dshape_( 2 )
54{
55  init() ;
56}
57
58STGrid::STGrid( const string infile )
59  : vshape_( 1 ), wshape_( 2 ), dshape_( 2 )
60{
61  init() ;
62
63  setFileIn( infile ) ;
64}
65
66STGrid::STGrid( const vector<string> infile )
67{
68  init() ;
69
70  setFileList( infile ) ;
71}
72
73void  STGrid::init()
74{
75  ifno_ = -1 ;
76  nx_ = -1 ;
77  ny_ = -1 ;
78  npol_ = 0 ;
79  nchan_ = 0 ;
80  nrow_ = 0 ;
81  cellx_ = 0.0 ;
82  celly_ = 0.0 ;
83  center_ = Vector<Double> ( 2, 0.0 ) ;
84  convType_ = "BOX" ;
85  wtype_ = "UNIFORM" ;
86  convSupport_ = -1 ;
87  userSupport_ = -1 ;
88  truncate_ = "";
89  gwidth_ = "";
90  jwidth_ = "";
91  convSampling_ = 100 ;
92  nprocessed_ = 0 ;
93  nchunk_ = 0 ;
94
95  // initialize user input
96  nxUI_ = -1 ;
97  nyUI_ = -1 ;
98  cellxUI_ = "" ;
99  cellyUI_ = "" ;
100  centerUI_ = "" ;
101  doclip_ = False ;
102}
103
104void STGrid::setFileIn( const string infile )
105{
106  nfile_ = 1 ;
107  String name( infile ) ;
108  infileList_.resize( nfile_ ) ;
109  infileList_[0] = String(infile) ;
110}
111
112void STGrid::setFileList( const vector<string> infile )
113{
114  nfile_ = infile.size() ;
115  infileList_.resize( nfile_ ) ;
116  for ( uInt i = 0 ; i < nfile_ ; i++ ) {
117    infileList_[i] = infile[i] ;
118  }
119}
120
121void STGrid::setPolList( vector<unsigned int> pols )
122{
123  pollist_.assign( Vector<uInt>( pols ) ) ;
124}
125
126void STGrid::setScanList( vector<unsigned int> scans )
127{
128  scanlist_.assign( Vector<uInt>( scans ) ) ;
129}
130
131void STGrid::setWeight( const string wType )
132{
133  wtype_ = String( wType ) ;
134  wtype_.upcase() ;
135}
136
137void STGrid::defineImage( int nx,
138                          int ny,
139                          string scellx,
140                          string scelly,
141                          string scenter )
142{
143  nxUI_ = (Int)nx ;
144  nyUI_ = (Int)ny ;
145  cellxUI_ = String( scellx ) ;
146  cellyUI_ = String( scelly ) ;
147  centerUI_ = String( scenter ) ;
148}
149 
150void STGrid::setFunc( string convType,
151                      int convSupport,
152                      string truncate,
153                      string gwidth,
154                      string jwidth )
155{
156  convType_ = String( convType ) ;
157  convType_.upcase() ;
158  userSupport_ = (Int)convSupport ;
159  truncate_ = String( truncate );
160  gwidth_ = String( gwidth );
161  jwidth_ = String( jwidth );
162}
163
164#define NEED_UNDERSCORES
165#if defined(NEED_UNDERSCORES)
166#define ggridsd ggridsd_
167#endif
168extern "C" {
169   void ggridsd(Double*,
170                const Complex*,
171                Int*,
172                Int*,
173                Int*,
174                const Int*,
175                const Int*,
176                const Float*,
177                Int*,
178                Int*,
179                Complex*,
180                Float*,
181                Int*,
182                Int*,
183                Int *,
184                Int *,
185                Int*,
186                Int*,
187                Float*,
188                Int*,
189                Int*,
190                Double*);
191}
192void STGrid::call_ggridsd( Array<Double> &xypos,
193                           Array<Complex> &spectra,
194                           Int &nvispol,
195                           Int &nvischan,
196                           Array<Int> &flagtra,
197                           Array<Int> &flagrow,
198                           Array<Float> &weight,
199                           Int &nrow,
200                           Int &irow,
201                           Array<Complex> &gdata,
202                           Array<Float> &gwgt,
203                           Int &nx,
204                           Int &ny,
205                           Int &npol,
206                           Int &nchan,
207                           Int &support,
208                           Int &sampling,
209                           Vector<Float> &convFunc,
210                           Int *chanMap,
211                           Int *polMap )
212{
213  // parameters for gridding
214  Int idopsf = 0 ;
215  Int len = npol*nchan ;
216  Double *sumw_p = new Double[len] ;
217  {
218    Double *work_p = sumw_p ;
219    for ( Int i = 0 ; i < len ; i++ ) {
220      *work_p = 0.0 ;
221      work_p++ ;
222    }
223  }
224
225  // prepare pointer
226  Bool deletePos, deleteData, deleteWgt, deleteFlag, deleteFlagR, deleteConv, deleteDataG, deleteWgtG ;
227  Double *xy_p = xypos.getStorage( deletePos ) ;
228  const Complex *values_p = spectra.getStorage( deleteData ) ;
229  const Int *flag_p = flagtra.getStorage( deleteFlag ) ;
230  const Int *rflag_p = flagrow.getStorage( deleteFlagR ) ;
231  const Float *wgt_p = weight.getStorage( deleteWgt ) ;
232  Complex *grid_p = gdata.getStorage( deleteDataG ) ;
233  Float *wgrid_p = gwgt.getStorage( deleteWgtG ) ;
234  Float *conv_p = convFunc.getStorage( deleteConv ) ;
235
236  // pass copy of irow to ggridsd since it will be modified in theroutine
237  Int irowCopy = irow ;
238     
239  // call ggridsd
240  ggridsd( xy_p,
241           values_p,
242           &nvispol,
243           &nvischan,
244           &idopsf,
245           flag_p,
246           rflag_p,
247           wgt_p,
248           &nrow,
249           &irowCopy,
250           grid_p,
251           wgrid_p,
252           &nx,
253           &ny,
254           &npol,
255           &nchan,
256           &support,
257           &sampling,
258           conv_p,
259           chanMap,
260           polMap,
261           sumw_p ) ;
262
263  // finalization
264  xypos.putStorage( xy_p, deletePos ) ;
265  spectra.freeStorage( values_p, deleteData ) ;
266  flagtra.freeStorage( flag_p, deleteFlag ) ;
267  flagrow.freeStorage( rflag_p, deleteFlagR ) ;
268  weight.freeStorage( wgt_p, deleteWgt ) ;
269  gdata.putStorage( grid_p, deleteDataG ) ;
270  gwgt.putStorage( wgrid_p, deleteWgtG ) ;
271  convFunc.putStorage( conv_p, deleteConv ) ;
272  delete sumw_p ;
273}
274
275#define NEED_UNDERSCORES
276#if defined(NEED_UNDERSCORES)
277#define ggridsd2 ggridsd2_
278#endif
279extern "C" {
280   void ggridsd2(Double*,
281                 const Complex*,
282                 Int*,
283                 Int*,
284                 Int*,
285                 const Int*,
286                 const Int*,
287                 const Float*,
288                 Int*,
289                 Int*,
290                 Complex*,
291                 Float*,
292                 Int*,
293                 Complex*,
294                 Float*,
295                 Float*,
296                 Complex*,
297                 Float*,
298                 Float*,
299                 Int*,
300                 Int*,
301                 Int *,
302                 Int *,
303                 Int*,
304                 Int*,
305                 Float*,
306                 Int*,
307                 Int*,
308                 Double*);
309}
310void STGrid::call_ggridsd2( Array<Double> &xypos,
311                            Array<Complex> &spectra,
312                            Int &nvispol,
313                            Int &nvischan,
314                            Array<Int> &flagtra,
315                            Array<Int> &flagrow,
316                            Array<Float> &weight,
317                            Int &nrow,
318                            Int &irow,
319                            Array<Complex> &gdata,
320                            Array<Float> &gwgt,
321                            Array<Int> &npoints,
322                            Array<Complex> &clipmin,
323                            Array<Float> &clipwmin,
324                            Array<Float> &clipcmin,
325                            Array<Complex> &clipmax,
326                            Array<Float> &clipwmax,
327                            Array<Float> &clipcmax,
328                            Int &nx,
329                            Int &ny,
330                            Int &npol,
331                            Int &nchan,
332                            Int &support,
333                            Int &sampling,
334                            Vector<Float> &convFunc,
335                            Int *chanMap,
336                            Int *polMap )
337{
338  // parameters for gridding
339  Int idopsf = 0 ;
340  Int len = npol*nchan ;
341  Double *sumw_p = new Double[len] ;
342  {
343    Double *work_p = sumw_p ;
344    for ( Int i = 0 ; i < len ; i++ ) {
345      *work_p = 0.0 ;
346      work_p++ ;
347    }
348  }
349
350  // prepare pointer
351  Bool deletePos, deleteData, deleteWgt, deleteFlag, deleteFlagR, deleteConv, deleteDataG, deleteWgtG, deleteNpts, deleteCMin, deleteCWMin, deleteCCMin, deleteCMax, deleteCWMax, deleteCCMax ;
352  Double *xy_p = xypos.getStorage( deletePos ) ;
353  const Complex *values_p = spectra.getStorage( deleteData ) ;
354  const Int *flag_p = flagtra.getStorage( deleteFlag ) ;
355  const Int *rflag_p = flagrow.getStorage( deleteFlagR ) ;
356  const Float *wgt_p = weight.getStorage( deleteWgt ) ;
357  Complex *grid_p = gdata.getStorage( deleteDataG ) ;
358  Float *wgrid_p = gwgt.getStorage( deleteWgtG ) ;
359  Float *conv_p = convFunc.getStorage( deleteConv ) ;
360  Int *npts_p = npoints.getStorage( deleteNpts ) ;
361  Complex *cmin_p = clipmin.getStorage( deleteCMin ) ;
362  Float *cwmin_p = clipwmin.getStorage( deleteCWMin ) ;
363  Float *ccmin_p = clipcmin.getStorage( deleteCCMin ) ;
364  Complex *cmax_p = clipmax.getStorage( deleteCMax ) ;
365  Float *cwmax_p = clipwmax.getStorage( deleteCWMax ) ;
366  Float *ccmax_p = clipcmax.getStorage( deleteCCMax ) ;
367
368  // pass copy of irow to ggridsd since it will be modified in theroutine
369  Int irowCopy = irow ;
370     
371  // call ggridsd
372  ggridsd2( xy_p,
373            values_p,
374            &nvispol,
375            &nvischan,
376            &idopsf,
377            flag_p,
378            rflag_p,
379            wgt_p,
380            &nrow,
381            &irowCopy,
382            grid_p,
383            wgrid_p,
384            npts_p,
385            cmin_p,
386            cwmin_p,
387            ccmin_p,
388            cmax_p,
389            cwmax_p,
390            ccmax_p,
391            &nx,
392            &ny,
393            &npol,
394            &nchan,
395            &support,
396            &sampling,
397            conv_p,
398            chanMap,
399            polMap,
400            sumw_p ) ;
401
402  // finalization
403  xypos.putStorage( xy_p, deletePos ) ;
404  spectra.freeStorage( values_p, deleteData ) ;
405  flagtra.freeStorage( flag_p, deleteFlag ) ;
406  flagrow.freeStorage( rflag_p, deleteFlagR ) ;
407  weight.freeStorage( wgt_p, deleteWgt ) ;
408  gdata.putStorage( grid_p, deleteDataG ) ;
409  gwgt.putStorage( wgrid_p, deleteWgtG ) ;
410  convFunc.putStorage( conv_p, deleteConv ) ;
411  clipmin.putStorage( cmin_p, deleteCMin ) ;
412  clipwmin.putStorage( cwmin_p, deleteCWMin ) ;
413  clipcmin.putStorage( ccmin_p, deleteCCMin ) ;
414  clipmax.putStorage( cmax_p, deleteCMax ) ;
415  clipwmax.putStorage( cwmax_p, deleteCWMax ) ;
416  clipcmax.putStorage( ccmax_p, deleteCCMax ) ;
417  delete sumw_p ;
418}
419
420void STGrid::grid()
421{
422  LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
423  double t0,t1 ;
424
425  // data selection
426  t0 = mathutil::gettimeofday_sec() ;
427  selectData() ;
428  t1 = mathutil::gettimeofday_sec() ;
429  os << LogIO::DEBUGGING << "selectData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
430
431  setupGrid() ;
432  setupArray() ;
433
434  if ( wtype_.compare("UNIFORM") != 0 &&
435       wtype_.compare("TINT") != 0 &&
436       wtype_.compare("TSYS") != 0 &&
437       wtype_.compare("TINTSYS") != 0 ) {
438    LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
439    os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
440    wtype_ = "UNIFORM" ;
441  }
442
443  // Warn if gauss or gjinc gridding with non-square cells
444  if ((cellx_ != celly_) && (convType_=="GAUSS"||convType_=="GJINC")) {
445    os << LogIO::WARN
446       << "The " << convType_ << " gridding doesn't support non-square grid." << endl
447       << "Result may be wrong." << LogIO::POST;
448  }
449
450  // grid parameter
451  os << LogIO::DEBUGGING ;
452  os << "----------" << endl ;
453  os << "Data selection summary" << endl ;
454  os << "   ifno = " << ifno_ << endl ;
455  os << "   pollist = " << pollist_ << endl ;
456  os << "   scanlist = " << scanlist_ << endl ;
457  os << "----------" << endl ;
458  os << "Grid parameter summary" << endl ;
459  os << "   (nx,ny) = (" << nx_ << "," << ny_ << ")" << endl ;
460  os << "   (cellx,celly) = (" << cellx_ << "," << celly_ << ")" << endl ;
461  os << "   center = " << center_ << endl ;
462  os << "   weighting = " << wtype_ << endl ;
463  os << "   convfunc = " << convType_ << endl;
464  if (convType_ == "GAUSS") {
465    os << "      gwidth = " << gwidth_ << endl;
466    os << "      truncate = " << truncate_ << endl;
467  }
468  else if (convType_ == "GJINC") {
469    os << "      gwidth = " << gwidth_ << endl;
470    os << "      jwidth = " << jwidth_ << endl;
471    os << "      truncate = " << truncate_ << endl;
472  }
473  else {
474    os << "      support = " << userSupport_ << endl;
475  }
476  os << "   doclip = " << (doclip_?"True":"False") << endl ;
477  os << "----------" << LogIO::POST ;
478  os << LogIO::NORMAL ;
479
480  if ( doclip_ )
481    gridPerRowWithClipping() ;
482  else
483    gridPerRow() ;
484}
485
486void STGrid::updateChunkShape()
487{
488  // TODO: nchunk_ must be determined from nchan_, npol_, and (nx_,ny_)
489  //       by considering data size to be allocated for ggridsd input/output
490  nchunk_ = 400 ;
491  nchunk_ = min( nchunk_, nrow_ ) ;
492  vshape_ = IPosition( 1, nchunk_ ) ;
493  wshape_ = IPosition( 2, nchan_, nchunk_ ) ;
494  dshape_ = IPosition( 2, 2, nchunk_ ) ;
495}
496
497struct STGChunk {
498  Int nrow ;
499  Array<Complex> spectra;
500  Array<Int> flagtra;
501  Array<Int> rflag;
502  Array<Float> weight;
503  Array<Double> direction;
504  STGChunk(IPosition const &wshape, IPosition const &vshape,
505           IPosition const &dshape)
506    : spectra(wshape), flagtra(wshape), rflag(vshape), weight(wshape),
507      direction(dshape)
508  { }
509};
510
511struct STCommonData {
512  Int gnx;
513  Int gny;
514  Int *chanMap;
515  Vector<Float> convFunc ;
516  Array<Complex> gdataArrC;
517  Array<Float> gwgtArr;
518  STCommonData(IPosition const &gshape, Array<Float> const &data)
519    : gdataArrC(gshape, 0.0), gwgtArr(data) {}
520};
521
522struct STCommonDataWithClipping {
523  Int gnx;
524  Int gny;
525  Int *chanMap;
526  Vector<Float> convFunc ;
527  Array<Complex> gdataArrC;
528  Array<Float> gwgtArr;
529  Array<Int> npoints ;
530  Array<Complex> clipMin ;
531  Array<Float> clipWMin ;
532  Array<Float> clipCMin ;
533  Array<Complex> clipMax ;
534  Array<Float> clipWMax ;
535  Array<Float> clipCMax ; 
536  STCommonDataWithClipping(IPosition const &gshape,
537                           IPosition const &pshape,
538                           Array<Float> const &data)
539    : gdataArrC(gshape, 0.0),
540      gwgtArr(data),
541      npoints(pshape, 0),
542      clipMin(gshape, Complex(FLT_MAX,0.0)),
543      clipWMin(gshape, 0.0),
544      clipCMin(gshape, 0.0),
545      clipMax(gshape, Complex(-FLT_MAX,0.0)),
546      clipWMax(gshape, 0.0),
547      clipCMax(gshape, 0.0)
548  {}
549};
550
551#define DO_AHEAD 3
552
553struct STContext {
554  STCommonData &common;
555  FIFO<STGChunk *, DO_AHEAD> queue;
556  STGrid *const self;
557  const Int pol;
558  STContext(STGrid *obj, STCommonData &common, Int pol)
559    : common(common), self(obj), pol(pol) {}
560};
561
562struct STContextWithClipping {
563  STCommonDataWithClipping &common;
564  FIFO<STGChunk *, DO_AHEAD> queue;
565  STGrid *const self;
566  const Int pol;
567  STContextWithClipping(STGrid *obj, STCommonDataWithClipping &common, Int pol)
568    : common(common), self(obj), pol(pol) {}
569};
570
571
572bool STGrid::produceChunk(void *ctx) throw(PCException)
573{
574  STContext &context = *(STContext *)ctx;
575  if ( context.self->nprocessed_ >= context.self->nrow_ ) {
576    return false;
577  }
578  STGChunk *chunk = new STGChunk(context.self->wshape_,
579                                 context.self->vshape_,
580                                 context.self->dshape_);
581
582  double t0 = mathutil::gettimeofday_sec() ;
583  chunk->nrow = context.self->getDataChunk(
584        context.self->wshape_, context.self->vshape_, context.self->dshape_,
585        chunk->spectra, chunk->direction,
586        chunk->flagtra, chunk->rflag, chunk->weight);
587  double t1 = mathutil::gettimeofday_sec() ;
588  context.self->eGetData_ += t1-t0 ;
589
590  context.queue.lock();
591  context.queue.put(chunk);
592  context.queue.unlock();
593  return true;
594}
595
596void STGrid::consumeChunk(void *ctx) throw(PCException)
597{
598  STContext &context = *(STContext *)ctx;
599  STGChunk *chunk = NULL;
600  try {
601    context.queue.lock();
602    chunk = context.queue.get();
603    context.queue.unlock();
604  } catch (FullException &e) {
605    context.queue.unlock();
606    // TODO: log error
607    throw PCException();
608  }
609
610  double t0, t1 ;
611  // world -> pixel
612  Array<Double> xypos( context.self->dshape_ ) ;
613  t0 = mathutil::gettimeofday_sec() ;
614  context.self->toPixel( chunk->direction, xypos ) ;
615  t1 = mathutil::gettimeofday_sec() ;
616  context.self->eToPixel_ += t1-t0 ;
617   
618  // call ggridsd
619  Int nvispol = 1 ;
620  Int irow = -1 ;
621  t0 = mathutil::gettimeofday_sec() ;
622  context.self->call_ggridsd( xypos,
623                chunk->spectra,
624                nvispol,
625                context.self->nchan_,
626                chunk->flagtra,
627                chunk->rflag,
628                chunk->weight,
629                chunk->nrow,
630                irow,
631                context.common.gdataArrC,
632                context.common.gwgtArr,
633                context.common.gnx,
634                context.common.gny,
635                context.self->npol_,
636                context.self->nchan_,
637                context.self->convSupport_,
638                context.self->convSampling_,
639                context.common.convFunc,
640                context.common.chanMap,
641                (Int*)&context.pol ) ;
642  t1 = mathutil::gettimeofday_sec() ;
643  context.self->eGGridSD_ += t1-t0 ;
644 
645  delete chunk;
646}
647
648void STGrid::gridPerRow()
649{
650  LogIO os( LogOrigin("STGrid", "gridPerRow", WHERE) ) ;
651  double t0, t1 ;
652
653
654  // grid data
655  // Extend grid plane with convSupport_
656  //   Int gnx = nx_+convSupport_*2 ;
657  //   Int gny = ny_+convSupport_*2 ;
658  Int gnx = nx_;
659  Int gny = ny_;
660
661  IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
662  // 2011/12/20 TN
663  // data_ and gwgtArr share storage
664  data_.resize( gshape ) ;
665  data_ = 0.0 ;
666  flag_.resize( gshape ) ;
667  flag_ = (uChar)0;
668  STCommonData common = STCommonData(gshape, data_);
669  common.gnx = gnx ;
670  common.gny = gny ;
671
672  // parameters for gridding
673  Int *chanMap = new Int[nchan_] ;
674  for ( Int i = 0 ; i < nchan_ ; i++ ) {
675    chanMap[i] = i ;
676  }
677  common.chanMap = chanMap;
678
679  // convolution kernel
680  t0 = mathutil::gettimeofday_sec() ;
681  setConvFunc( common.convFunc ) ;
682  t1 = mathutil::gettimeofday_sec() ;
683  os << LogIO::DEBUGGING << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
684
685  // for performance check
686  eGetData_ = 0.0 ;
687  eToPixel_ = 0.0 ;
688  eGGridSD_ = 0.0 ;
689  double eInitPol = 0.0 ;
690
691  for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) {
692    initTable( ifile ) ;
693
694    os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ;   
695    Broker broker = Broker(produceChunk, consumeChunk);
696    for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
697      t0 = mathutil::gettimeofday_sec() ;
698      initPol( ipol ) ; // set ptab_ and attach()
699      t1 = mathutil::gettimeofday_sec() ;
700      eInitPol += t1-t0 ;
701     
702      STContext context(this, common, ipol);
703     
704      os << "start pol " << ipol << LogIO::POST ;
705     
706      nprocessed_ = 0 ;
707#if 1
708      broker.runProducerAsMasterThread(&context, DO_AHEAD);
709#else
710      for (;;) {
711        bool produced = produceChunk(&context);
712        if (! produced) {
713          break;
714        }
715        consumeChunk(&context);
716      }
717#endif
718
719      os << "end pol " << ipol << LogIO::POST ;
720
721    }
722    os << "end table " << ifile << LogIO::POST ;   
723  }
724  os << LogIO::DEBUGGING << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ;
725  os << LogIO::DEBUGGING << "getData: elapsed time is " << eGetData_-eToInt-eGetWeight << " sec." << LogIO::POST ;
726  os << LogIO::DEBUGGING << "toPixel: elapsed time is " << eToPixel_ << " sec." << LogIO::POST ;
727  os << LogIO::DEBUGGING << "ggridsd: elapsed time is " << eGGridSD_ << " sec." << LogIO::POST ;
728  os << LogIO::DEBUGGING << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
729  os << LogIO::DEBUGGING << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
730 
731  delete chanMap ;
732
733  // set data
734  setData( common.gdataArrC, common.gwgtArr ) ;
735
736}
737
738void STGrid::consumeChunkWithClipping(void *ctx) throw(PCException)
739{
740  STContextWithClipping &context = *(STContextWithClipping *)ctx;
741  STGChunk *chunk = NULL;
742  try {
743    context.queue.lock();
744    chunk = context.queue.get();
745    context.queue.unlock();
746  } catch (FullException &e) {
747    context.queue.unlock();
748    // TODO: log error
749    throw PCException();
750  }
751
752  double t0, t1 ;
753  // world -> pixel
754  Array<Double> xypos( context.self->dshape_ ) ;
755  t0 = mathutil::gettimeofday_sec() ;
756  context.self->toPixel( chunk->direction, xypos ) ;
757  t1 = mathutil::gettimeofday_sec() ;
758  context.self->eToPixel_ += t1-t0 ;
759   
760  // call ggridsd
761  Int nvispol = 1 ;
762  Int irow = -1 ;
763  t0 = mathutil::gettimeofday_sec() ;
764  context.self->call_ggridsd2( xypos,
765                chunk->spectra,
766                nvispol,
767                context.self->nchan_,
768                chunk->flagtra,
769                chunk->rflag,
770                chunk->weight,
771                chunk->nrow,
772                irow,
773                context.common.gdataArrC,
774                context.common.gwgtArr,
775                context.common.npoints,
776                context.common.clipMin,
777                context.common.clipWMin,
778                context.common.clipCMin,
779                context.common.clipMax,
780                context.common.clipWMax,
781                context.common.clipCMax,
782                context.common.gnx,
783                context.common.gny,
784                context.self->npol_,
785                context.self->nchan_,
786                context.self->convSupport_,
787                context.self->convSampling_,
788                context.common.convFunc,
789                context.common.chanMap,
790                (Int*)&context.pol ) ;
791  t1 = mathutil::gettimeofday_sec() ;
792  context.self->eGGridSD_ += t1-t0 ;
793 
794  delete chunk;
795}
796
797void STGrid::gridPerRowWithClipping()
798{
799  LogIO os( LogOrigin("STGrid", "gridPerRowWithClipping", WHERE) ) ;
800  double t0, t1 ;
801
802
803  // grid data
804  // Extend grid plane with convSupport_
805  //   Int gnx = nx_+convSupport_*2 ;
806  //   Int gny = ny_+convSupport_*2 ;
807  Int gnx = nx_;
808  Int gny = ny_;
809
810  IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
811  IPosition pshape( 3, gnx, gny, npol_ ) ;
812  // 2011/12/20 TN
813  // data_ and gwgtArr share storage
814  data_.resize( gshape ) ;
815  data_ = 0.0 ;
816  flag_.resize( gshape ) ;
817  flag_ = (uChar)0;
818  STCommonDataWithClipping common = STCommonDataWithClipping( gshape,
819                                                              pshape,
820                                                              data_ ) ;
821  common.gnx = gnx ;
822  common.gny = gny ;
823
824  // parameters for gridding
825  Int *chanMap = new Int[nchan_] ;
826  for ( Int i = 0 ; i < nchan_ ; i++ ) {
827    chanMap[i] = i ;
828  }
829  common.chanMap = chanMap;
830
831  // convolution kernel
832  t0 = mathutil::gettimeofday_sec() ;
833  setConvFunc( common.convFunc ) ;
834  t1 = mathutil::gettimeofday_sec() ;
835  os << LogIO::DEBUGGING << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
836
837  // for performance check
838  eGetData_ = 0.0 ;
839  eToPixel_ = 0.0 ;
840  eGGridSD_ = 0.0 ;
841  double eInitPol = 0.0 ;
842
843  for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) {
844    initTable( ifile ) ;
845
846    os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ;   
847    Broker broker = Broker(produceChunk, consumeChunkWithClipping);
848    for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
849      t0 = mathutil::gettimeofday_sec() ;
850      initPol( ipol ) ; // set ptab_ and attach()
851      t1 = mathutil::gettimeofday_sec() ;
852      eInitPol += t1-t0 ;
853     
854      STContextWithClipping context(this, common, ipol);
855     
856      os << "start pol " << ipol << LogIO::POST ;
857     
858      nprocessed_ = 0 ;
859#if 1
860      broker.runProducerAsMasterThread(&context, DO_AHEAD);
861#else
862      for (;;) {
863        bool produced = produceChunk(&context);
864        if (! produced) {
865          break;
866        }
867        consumeChunkWithClipping(&context);
868      }
869#endif
870
871      os << "end pol " << ipol << LogIO::POST ;
872
873    }
874    os << "end table " << ifile << LogIO::POST ;   
875  }
876  os << LogIO::DEBUGGING << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ;
877  os << LogIO::DEBUGGING << "getData: elapsed time is " << eGetData_-eToInt-eGetWeight << " sec." << LogIO::POST ;
878  os << LogIO::DEBUGGING << "toPixel: elapsed time is " << eToPixel_ << " sec." << LogIO::POST ;
879  os << LogIO::DEBUGGING << "ggridsd2: elapsed time is " << eGGridSD_ << " sec." << LogIO::POST ;
880  os << LogIO::DEBUGGING << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
881  os << LogIO::DEBUGGING << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
882 
883  delete chanMap ;
884
885  // clip min and max in each grid
886//   os << "BEFORE CLIPPING" << LogIO::POST ;
887//   os << "gdataArrC=" << common.gdataArrC << LogIO::POST ;
888//   os << "gwgtArr=" << common.gwgtArr << LogIO::POST ;
889  t0 = mathutil::gettimeofday_sec() ;
890  clipMinMax( common.gdataArrC,
891              common.gwgtArr,
892              common.npoints,
893              common.clipMin,
894              common.clipWMin,
895              common.clipCMin,
896              common.clipMax,
897              common.clipWMax,
898              common.clipCMax ) ;
899  t1 = mathutil::gettimeofday_sec() ;
900  os << LogIO::DEBUGGING << "clipMinMax: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
901//   os << "AFTER CLIPPING" << LogIO::POST ;
902//   os << "gdataArrC=" << common.gdataArrC << LogIO::POST ;
903//   os << "gwgtArr=" << common.gwgtArr << LogIO::POST ;
904
905  // set data
906  setData( common.gdataArrC, common.gwgtArr ) ;
907
908}
909
910void STGrid::clipMinMax( Array<Complex> &grid,
911                         Array<Float> &weight,
912                         Array<Int> &npoints,
913                         Array<Complex> &clipmin,
914                         Array<Float> &clipwmin,
915                         Array<Float> &clipcmin,
916                         Array<Complex> &clipmax,
917                         Array<Float> &clipwmax,
918                         Array<Float> &clipcmax )
919{
920  //LogIO os( LogOrigin("STGrid","clipMinMax",WHERE) ) ;
921
922  // prepare pointers
923  Bool delG, delW, delNP, delCMin, delCWMin, delCCMin, delCMax, delCWMax, delCCMax ;
924  Complex *grid_p = grid.getStorage( delG ) ;
925  Float *wgt_p = weight.getStorage( delW ) ;
926  const Int *npts_p = npoints.getStorage( delNP ) ;
927  const Complex *cmin_p = clipmin.getStorage( delCMin ) ;
928  const Float *cwmin_p = clipwmin.getStorage( delCWMin ) ;
929  const Float *ccmin_p = clipcmin.getStorage( delCCMin ) ;
930  const Complex *cmax_p = clipmax.getStorage( delCMax ) ;
931  const Float *cwmax_p = clipwmax.getStorage( delCWMax ) ;
932  const Float *ccmax_p = clipcmax.getStorage( delCCMax ) ;
933
934  const IPosition &gshape = grid.shape() ;
935  long offset = gshape[0] * gshape[1] * gshape[2] ; // nx * ny * npol
936  Int nchan = gshape[3] ;
937  long origin = nchan * offset ;
938  for ( long i = 0 ; i < offset ; i++ ) {
939    if ( *npts_p > 2 ) {
940      for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
941        // clip minimum and maximum
942        *grid_p -= (*cmin_p)*(*cwmin_p)*(*ccmin_p)
943          + (*cmax_p)*(*cwmax_p)*(*ccmax_p) ;
944        *wgt_p -= (*cwmin_p)*(*ccmin_p)
945          + (*cwmax_p)*(*ccmax_p) ;
946       
947        grid_p += offset ;
948        wgt_p += offset ;
949        cmin_p += offset ;
950        cwmin_p += offset ;
951        ccmin_p += offset ;
952        cmax_p += offset ;
953        cwmax_p += offset ;
954        ccmax_p += offset ;
955      }
956      grid_p -= origin ;
957      wgt_p -= origin ;
958      cmin_p -= origin ;
959      cwmin_p -= origin ;
960      ccmin_p -= origin ;
961      cmax_p -= origin ;
962      cwmax_p -= origin ;
963      ccmax_p -= origin ;
964    }
965    grid_p++ ;
966    wgt_p++ ;
967    npts_p++ ;
968    cmin_p++ ;
969    cwmin_p++ ;
970    ccmin_p++ ;
971    cmax_p++ ;
972    cwmax_p++ ;
973    ccmax_p++ ;
974  }
975  grid_p -= offset ;
976  wgt_p -= offset ;
977  npts_p -= offset ;
978  cmin_p -= offset ;
979  cwmin_p -= offset ;
980  ccmin_p -= offset ;
981  cmax_p -= offset ;
982  cwmax_p -= offset ;
983  ccmax_p -= offset ; 
984
985  // finalization
986  grid.putStorage( grid_p, delG ) ;
987  weight.putStorage( wgt_p, delW ) ;
988  npoints.freeStorage( npts_p, delNP ) ;
989  clipmin.freeStorage( cmin_p, delCMin ) ;
990  clipwmin.freeStorage( cwmin_p, delCWMin ) ;
991  clipcmin.freeStorage( ccmin_p, delCCMin ) ;
992  clipmax.freeStorage( cmax_p, delCMax ) ;
993  clipwmax.freeStorage( cwmax_p, delCWMax ) ;
994  clipcmax.freeStorage( ccmax_p, delCCMax ) ;
995}
996
997void STGrid::initPol( Int ipol )
998{
999  LogIO os( LogOrigin("STGrid","initPol",WHERE) ) ;
1000  if ( npolOrg_ == 1 ) {
1001    os << "single polarization data." << LogIO::POST ;
1002    ptab_ = tab_ ;
1003  }
1004  else
1005    ptab_ = tab_( tab_.col("POLNO") == pollist_[ipol] ) ;
1006
1007  attach( ptab_ ) ;
1008}
1009
1010void STGrid::initTable( uInt idx )
1011{
1012  tab_ = tableList_[idx] ;
1013  nrow_ = rows_[idx] ;
1014  updateChunkShape() ;
1015}
1016
1017void STGrid::setData( Array<Complex> &gdata,
1018                      Array<Float> &gwgt )
1019{
1020  // 2011/12/20 TN
1021  // gwgt and data_ share storage
1022  LogIO os( LogOrigin("STGrid","setData",WHERE) ) ;
1023  double t0, t1 ;
1024  t0 = mathutil::gettimeofday_sec() ;
1025  uInt len = data_.nelements() ;
1026  Bool b1, b2, b3 ;
1027  const Complex *gdata_p = gdata.getStorage( b1 ) ;
1028  Float *gwgt_p = gwgt.getStorage( b2 ) ; // storage shared with data_
1029  uChar *gflg_p = flag_.getStorage( b3 ) ;
1030  for ( uInt i = 0 ; i < len ; i++ ) {
1031    if (gwgt_p[i] > 0.0) {
1032      gwgt_p[i] = (gdata_p[i]).real() / gwgt_p[i];
1033      gflg_p[i] = (uChar)0;
1034    }
1035    else {
1036      gflg_p[i] = (uChar)1;
1037    }
1038  }
1039  gdata.freeStorage( gdata_p, b1 ) ;
1040  data_.putStorage( gwgt_p, b2 ) ;
1041  flag_.putStorage( gflg_p, b3 ) ;
1042  t1 = mathutil::gettimeofday_sec() ;
1043  os << LogIO::DEBUGGING << "setData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
1044}
1045
1046void STGrid::setupGrid()
1047{
1048  Double xmin,xmax,ymin,ymax ;
1049  mapExtent( xmin, xmax, ymin, ymax ) ;
1050 
1051  setupGrid( nxUI_, nyUI_, cellxUI_, cellyUI_,
1052             xmin, xmax, ymin, ymax, centerUI_ ) ;
1053}
1054
1055void STGrid::setupGrid( Int &nx,
1056                        Int &ny,
1057                        String &cellx,
1058                        String &celly,
1059                        Double &xmin,
1060                        Double &xmax,
1061                        Double &ymin,
1062                        Double &ymax,
1063                        String &center )
1064{
1065  LogIO os( LogOrigin("STGrid","setupGrid",WHERE) ) ;
1066  //cout << "nx=" << nx << ", ny=" << ny << endl ;
1067
1068  // center position
1069  if ( center.size() == 0 ) {
1070    center_(0) = 0.5 * ( xmin + xmax ) ;
1071    center_(1) = 0.5 * ( ymin + ymax ) ;
1072  }
1073  else {
1074    String::size_type pos0 = center.find( " " ) ;
1075    if ( pos0 == String::npos ) {
1076      throw AipsError( "bad string format in parameter center" ) ;
1077    }
1078    String::size_type pos1 = center.find( " ", pos0+1 ) ;
1079    String typestr, xstr, ystr ;
1080    if ( pos1 != String::npos ) {
1081      typestr = center.substr( 0, pos0 ) ;
1082      xstr = center.substr( pos0+1, pos1-pos0 ) ;
1083      ystr = center.substr( pos1+1 ) ;
1084      // todo: convert to J2000 (or direction ref for DIRECTION column)
1085    }
1086    else {
1087      typestr = "J2000" ;
1088      xstr = center.substr( 0, pos0 ) ;
1089      ystr = center.substr( pos0+1 ) ;
1090    }
1091    QuantumHolder qh ;
1092    String err ;
1093    qh.fromString( err, xstr ) ;
1094    Quantum<Double> xcen = qh.asQuantumDouble() ;
1095    qh.fromString( err, ystr ) ;
1096    Quantum<Double> ycen = qh.asQuantumDouble() ;
1097    center_(0) = xcen.getValue( "rad" ) ;
1098    center_(1) = ycen.getValue( "rad" ) ;
1099    double base = 0.5 * (xmin + xmax) ;
1100    int maxrotate = 1 ;
1101    int nelem = 2 * maxrotate + 1 ;
1102    double *sep = new double[nelem] ;
1103    for ( int i = 0 ; i < nelem ; i++ )
1104      sep[i] = abs(base - center_[0] - (i-maxrotate) * C::_2pi) ;
1105//     os << "sep[0]=" << sep[0] << endl 
1106//        << "sep[1]=" << sep[1] << endl
1107//        << "sep[2]=" << sep[2] << LogIO::POST ;
1108    int idx = 0 ;
1109    base = sep[0] ;
1110    int nrotate = 0 ;
1111    while ( idx < nelem ) {
1112      if ( base > sep[idx] ) {
1113        base = sep[idx] ;
1114        nrotate = idx ;
1115      }
1116      idx++ ;
1117    }
1118    delete sep ;
1119    nrotate -= maxrotate ;
1120//     os << "nrotate = " << nrotate << LogIO::POST ;
1121    center_[0] += nrotate * C::_2pi ;
1122  }
1123//   os << "xmin=" << xmin << LogIO::POST ;
1124//   os << "center_=" << center_ << LogIO::POST ;
1125
1126  nx_ = nx ;
1127  ny_ = ny ;
1128  if ( nx < 0 && ny > 0 ) {
1129    nx_ = ny ;
1130    ny_ = ny ;
1131  }
1132  if ( ny < 0 && nx > 0 ) {
1133    nx_ = nx ;
1134    ny_ = nx ;
1135  }
1136
1137  //Double wx = xmax - xmin ;
1138  //Double wy = ymax - ymin ;
1139  Double wx = max( abs(xmax-center_(0)), abs(xmin-center_(0)) ) * 2 ;
1140  Double wy = max( abs(ymax-center_(1)), abs(ymin-center_(1)) ) * 2 ;
1141
1142  Quantum<Double> qcellx ;
1143  Quantum<Double> qcelly ;
1144  //cout << "nx_ = " << nx_ << ",  ny_ = " << ny_ << endl ;
1145  if ( cellx.size() != 0 && celly.size() != 0 ) {
1146    readQuantity( qcellx, cellx ) ;
1147    readQuantity( qcelly, celly ) ;
1148  }
1149  else if ( celly.size() != 0 ) {
1150    os << "Using celly to x-axis..." << LogIO::POST ;
1151    readQuantity( qcelly, celly ) ;
1152    qcellx = qcelly ;
1153  }
1154  else if ( cellx.size() != 0 ) {
1155    os << "Using cellx to y-axis..." << LogIO::POST ;
1156    readQuantity( qcellx, cellx ) ;
1157    qcelly = qcellx ;
1158  }
1159  else {
1160    if ( nx_ < 0 ) {
1161      os << "No user preference in grid setting. Using default..." << LogIO::POST ;
1162      readQuantity( qcellx, "1.0arcmin" ) ;
1163      qcelly = qcellx ;
1164    }
1165    else {
1166      if ( wx == 0.0 ) {
1167        os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
1168        wx = 0.00290888 ;
1169      }
1170      if ( wy == 0.0 ) {
1171        os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
1172        wy = 0.00290888 ;
1173      }
1174      if (nx_ > 1) {
1175        qcellx = Quantum<Double>(wx / (nx_ - 1) * cos(center_[1]), "rad");
1176      }
1177      else {
1178        qcellx = Quantum<Double>( 1.1f * wx / nx_ * cos(center_[1]), "rad" );
1179      }
1180      if (ny_ > 1) {
1181        qcelly = Quantum<Double>(wy / (ny_ - 1), "rad");
1182      }
1183      else {
1184        qcelly = Quantum<Double>( 1.1f * wy / ny_, "rad" );
1185      }
1186    }
1187  }
1188  cellx_ = qcellx.getValue( "rad" ) ;
1189  celly_ = qcelly.getValue( "rad" ) ;
1190  //os << "cellx_=" << cellx_ << ", celly_=" << celly_ << ", cos("<<center_(1)<<")=" << cos(center_(1)) << LogIO::POST ;
1191  if ( nx_ < 0 ) {
1192    if ( wx == 0.0 ) {
1193      os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
1194      wx = 0.00290888 ;
1195    }
1196    if ( wy == 0.0 ) {
1197      os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
1198      wy = 0.00290888 ;
1199    }
1200    nx_ = Int( ceil( wx/(cellx_/cos(center_[1])) ) ) ;
1201    ny_ = Int( ceil( wy/celly_ ) ) ;
1202  }
1203
1204  // create DirectionCoordinate
1205  Matrix<Double> xform(2,2) ;
1206  xform = 0.0 ;
1207  xform.diagonal() = 1.0 ;
1208  dircoord_ = new DirectionCoordinate(MDirection::J2000,
1209                                      Projection( Projection::SIN ),
1210                                      center_[0], center_[1],
1211                                      -cellx_, celly_,
1212                                      xform,
1213                                      0.5*Double(nx_-1),
1214                                      0.5*Double(ny_-1)) ;
1215}
1216
1217void STGrid::mapExtent( Double &xmin, Double &xmax,
1218                        Double &ymin, Double &ymax )
1219{
1220  //LogIO os( LogOrigin("STGrid","mapExtent",WHERE) ) ;
1221  directionCol_.attach( tableList_[0], "DIRECTION" ) ;
1222  Matrix<Double> direction = directionCol_.getColumn() ;
1223  //os << "dirCol.nrow() = " << dirCol.nrow() << LogIO::POST ;
1224  Vector<Double> ra( direction.row(0) ) ;
1225  mathutil::rotateRA( ra ) ;
1226  minMax( xmin, xmax, direction.row( 0 ) ) ;
1227  minMax( ymin, ymax, direction.row( 1 ) ) ;
1228  Double amin, amax, bmin, bmax ;
1229  for ( uInt i = 1 ; i < nfile_ ; i++ ) {
1230    directionCol_.attach( tableList_[i], "DIRECTION" ) ;
1231    direction.assign( directionCol_.getColumn() ) ;
1232    //os << "dirCol.nrow() = " << dirCol.nrow() << LogIO::POST ;
1233    // to make contiguous RA distribution (no 2pi jump)
1234    Vector<Double> ra( direction.row(0) ) ;
1235    mathutil::rotateRA( ra ) ;
1236    minMax( amin, amax, direction.row( 0 ) ) ;
1237    minMax( bmin, bmax, direction.row( 1 ) ) ;
1238    xmin = min( xmin, amin ) ;
1239    xmax = max( xmax, amax ) ;
1240    ymin = min( ymin, bmin ) ;
1241    ymax = max( ymax, bmax ) ;
1242  }
1243  //os << "(xmin,xmax)=(" << xmin << "," << xmax << ")" << LogIO::POST ;
1244  //os << "(ymin,ymax)=(" << ymin << "," << ymax << ")" << LogIO::POST ;
1245}
1246
1247void STGrid::table( Table &tab, uInt i )
1248{
1249  if ( i < nfile_ )
1250    tab = Table( infileList_[i] ) ;
1251}
1252
1253void STGrid::selectData()
1254{
1255  LogIO os( LogOrigin("STGrid","selectData",WHERE) ) ;   
1256  Int ifno = ifno_ ;
1257  tableList_.resize( nfile_ ) ;
1258  if ( ifno_ == -1 ) {
1259    //Table taborg( infileList_[0] ) ;
1260    Table taborg ;
1261    table( taborg, 0 ) ;
1262    ROScalarColumn<uInt> ifnoCol( taborg, "IFNO" ) ;
1263    ifno_ = ifnoCol( 0 ) ;
1264    os << LogIO::WARN
1265       << "IFNO is not given. Using default IFNO: " << ifno_ << LogIO::POST ;
1266  }
1267  for ( uInt i = 0 ; i < nfile_ ; i++ ) {
1268    //Table taborg( infileList_[i] ) ;
1269    Table taborg ;
1270    table( taborg, i ) ;
1271    TableExprNode node ;
1272    if ( ifno != -1 || isMultiIF( taborg ) ) {
1273      os << "apply selection on IFNO" << LogIO::POST ;
1274      node = taborg.col("IFNO") == ifno_ ;
1275    }
1276    if ( scanlist_.size() > 0 ) {
1277      os << "apply selection on SCANNO" << LogIO::POST ;
1278      node = node && taborg.col("SCANNO").in( scanlist_ ) ;
1279    }
1280    if ( node.isNull() ) {
1281      tableList_[i] = taborg ;
1282    }
1283    else {
1284      tableList_[i] = taborg( node ) ;
1285    }
1286    os << LogIO::DEBUGGING << "tableList_[" << i << "].nrow()=" << tableList_[i].nrow() << LogIO::POST ;
1287    if ( tableList_[i].nrow() == 0 ) {
1288      os << LogIO::SEVERE
1289         << "No corresponding rows for given selection: IFNO " << ifno_ ;
1290      if ( scanlist_.size() > 0 )
1291        os << " SCANNO " << scanlist_ ;
1292      os << LogIO::EXCEPTION ;
1293    }
1294  }
1295}
1296
1297Bool STGrid::isMultiIF( Table &tab )
1298{
1299  ROScalarColumn<uInt> ifnoCol( tab, "IFNO" ) ;
1300  Vector<uInt> ifnos = ifnoCol.getColumn() ;
1301  return anyNE( ifnos, ifnos[0] ) ;
1302}
1303
1304void STGrid::attach( Table &tab )
1305{
1306  // attach to table
1307  spectraCol_.attach( tab, "SPECTRA" ) ;
1308  flagtraCol_.attach( tab, "FLAGTRA" ) ;
1309  directionCol_.attach( tab, "DIRECTION" ) ;
1310  flagRowCol_.attach( tab, "FLAGROW" ) ;
1311  tsysCol_.attach( tab, "TSYS" ) ;
1312  intervalCol_.attach( tab, "INTERVAL" ) ;
1313}
1314
1315Int STGrid::getDataChunk(
1316                         IPosition const &wshape,
1317                         IPosition const &vshape,
1318                         IPosition const &/* dshape */,
1319                         Array<Complex> &spectra,
1320                         Array<Double> &direction,
1321                         Array<Int> &flagtra,
1322                         Array<Int> &rflag,
1323                         Array<Float> &weight )
1324{
1325  LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
1326
1327  Array<Float> spectraF_(wshape);
1328  Array<uChar> flagtraUC_(wshape);
1329  Array<uInt> rflagUI_(vshape);
1330  Int nrow = getDataChunk( spectraF_, direction, flagtraUC_, rflagUI_, weight ) ;
1331  if ( nrow < nchunk_ ) {
1332    spectra.resize( spectraF_.shape() ) ;
1333    flagtra.resize( flagtraUC_.shape() ) ;
1334    rflag.resize( rflagUI_.shape() ) ;
1335  }
1336  double t0, t1 ;
1337  t0 = mathutil::gettimeofday_sec() ;
1338  convertArray( spectra, spectraF_ ) ;
1339  toInt( flagtraUC_, flagtra ) ;
1340  toInt( rflagUI_, rflag ) ;
1341  t1 = mathutil::gettimeofday_sec() ;
1342  eToInt = t1 - t0 ;
1343 
1344  return nrow ;
1345}
1346
1347#if 0
1348Int STGrid::getDataChunk( Array<Complex> &spectra,
1349                          Array<Double> &direction,
1350                          Array<Int> &flagtra,
1351                          Array<Int> &rflag,
1352                          Array<Float> &weight )
1353{
1354  LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
1355  Int nrow = getDataChunk( spectraF_, direction, flagtraUC_, rflagUI_, weight ) ;
1356  if ( nrow < nchunk_ ) {
1357    spectra.resize( spectraF_.shape() ) ;
1358    flagtra.resize( flagtraUC_.shape() ) ;
1359    rflag.resize( rflagUI_.shape() ) ;
1360  }
1361  double t0, t1 ;
1362  t0 = mathutil::gettimeofday_sec() ;
1363  convertArray( spectra, spectraF_ ) ;
1364  toInt( flagtraUC_, flagtra ) ;
1365  toInt( rflagUI_, rflag ) ;
1366  t1 = mathutil::gettimeofday_sec() ;
1367  eToInt = t1 - t0 ;
1368 
1369  return nrow ;
1370}
1371#endif
1372
1373Int STGrid::getDataChunk( Array<Float> &spectra,
1374                          Array<Double> &direction,
1375                          Array<uChar> &flagtra,
1376                          Array<uInt> &rflag,
1377                          Array<Float> &weight )
1378{
1379  LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
1380  Int nrow = spectra.shape()[1] ;
1381  Int remainingRow = nrow_ - nprocessed_ ;
1382  if ( remainingRow < nrow ) {
1383    nrow = remainingRow ;
1384    IPosition mshape( 2, nchan_, nrow ) ;
1385    IPosition vshape( 1, nrow ) ;
1386    spectra.resize( mshape ) ;
1387    flagtra.resize( mshape ) ;
1388    direction.resize( IPosition(2,2,nrow) ) ;
1389    rflag.resize( vshape ) ;
1390    weight.resize( mshape ) ;
1391  }
1392  // 2011/12/22 TN
1393  // tsys shares its storage with weight
1394  Array<Float> tsys( weight ) ;
1395  Array<Double> tint( rflag.shape() ) ;
1396
1397  Vector<uInt> rflagVec( rflag ) ;
1398  Vector<Double> tintVec( tint ) ;
1399
1400  RefRows rows( nprocessed_, nprocessed_+nrow-1, 1 ) ;
1401  //os<<LogIO::DEBUGGING<<"nprocessed_="<<nprocessed_<<": rows.nrows()="<<rows.nrows()<<LogIO::POST ;
1402  spectraCol_.getColumnCells( rows, spectra ) ;
1403  flagtraCol_.getColumnCells( rows, flagtra ) ;
1404  directionCol_.getColumnCells( rows, direction ) ;
1405  // to make contiguous RA distribution (no 2pi jump)
1406  Vector<Double> v( Matrix<Double>(direction).row(0) ) ;
1407  mathutil::rotateRA( v ) ;
1408  flagRowCol_.getColumnCells( rows, rflagVec ) ;
1409  intervalCol_.getColumnCells( rows, tintVec ) ;
1410  Vector<Float> tsysTemp = tsysCol_( nprocessed_ ) ;
1411  if ( tsysTemp.nelements() == (uInt)nchan_ )
1412    tsysCol_.getColumnCells( rows, tsys ) ;
1413  else
1414    tsys = tsysTemp[0] ;
1415
1416  double t0,t1 ;
1417  t0 = mathutil::gettimeofday_sec() ;
1418  getWeight( weight, tsys, tint ) ;
1419  t1 = mathutil::gettimeofday_sec() ;
1420  eGetWeight += t1-t0 ;
1421
1422  nprocessed_ += nrow ;
1423 
1424  return nrow ;
1425}
1426
1427void STGrid::setupArray()
1428{
1429  LogIO os( LogOrigin("STGrid","setupArray",WHERE) ) ;
1430  ROScalarColumn<uInt> polnoCol( tableList_[0], "POLNO" ) ;
1431  Vector<uInt> pols = polnoCol.getColumn() ;
1432  //os << pols << LogIO::POST ;
1433  Vector<uInt> pollistOrg ;
1434  npolOrg_ = 0 ;
1435  uInt polno ;
1436  for ( uInt i = 0 ; i < polnoCol.nrow() ; i++ ) {
1437    //polno = polnoCol( i ) ;
1438    polno = pols( i ) ;
1439    if ( allNE( pollistOrg, polno ) ) {
1440      pollistOrg.resize( npolOrg_+1, True ) ;
1441      pollistOrg[npolOrg_] = polno ;
1442      npolOrg_++ ;
1443    }
1444  }
1445  if ( pollist_.size() == 0 )
1446    pollist_ = pollistOrg ;
1447  else {
1448    Vector<uInt> newlist ;
1449    uInt newsize = 0 ;
1450    for ( uInt i = 0 ; i < pollist_.size() ; i++ ) {
1451      if ( anyEQ( pollistOrg, pollist_[i] ) ) {
1452        newlist.resize( newsize+1, True ) ;
1453        newlist[newsize] = pollist_[i] ;
1454        newsize++ ;
1455      }
1456    }
1457    pollist_.assign( newlist ) ;
1458  }
1459  npol_ = pollist_.size() ;
1460  if ( npol_ == 0 ) {
1461    os << LogIO::SEVERE << "Empty pollist" << LogIO::EXCEPTION ;
1462  }
1463  rows_.resize( nfile_ ) ;
1464  for ( uInt i = 0 ; i < nfile_ ; i++ ) {
1465    rows_[i] = tableList_[i].nrow() / npolOrg_ ;
1466    //if ( nrow_ < rows_[i] )
1467    //  nrow_ = rows_[i] ;
1468  }
1469  flagtraCol_.attach( tableList_[0], "FLAGTRA" ) ;
1470  nchan_ = flagtraCol_( 0 ).nelements() ;
1471//   os << "npol_ = " << npol_ << "(" << pollist_ << ")" << endl
1472//      << "nchan_ = " << nchan_ << endl
1473//      << "nrow_ = " << nrow_ << LogIO::POST ;
1474}
1475
1476void STGrid::getWeight( Array<Float> &w,
1477                              Array<Float> &tsys,
1478                              Array<Double> &tint )
1479{
1480  LogIO os( LogOrigin("STGrid","getWeight",WHERE) ) ;
1481
1482  // 2011/12/22 TN
1483  // w (weight) and tsys share storage
1484  IPosition refShape = tsys.shape() ;
1485  Int nchan = refShape[0] ;
1486  Int nrow = refShape[1] ;
1487//   os << "nchan=" << nchan << ", nrow=" << nrow << LogIO::POST ;
1488//   os << "w.shape()=" << w.shape() << endl
1489//      << "tsys.shape()=" << tsys.shape() << endl
1490//      << "tint.shape()=" << tint.shape() << LogIO::POST ;
1491
1492  // set weight
1493  if ( wtype_.compare( "UNIFORM" ) == 0 ) {
1494    w = 1.0 ;
1495  }
1496  else if ( wtype_.compare( "TINT" ) == 0 ) {
1497    Bool b0, b1 ;
1498    Float *w_p = w.getStorage( b0 ) ;
1499    Float *w0_p = w_p ;
1500    const Double *ti_p = tint.getStorage( b1 ) ;
1501    const Double *w1_p = ti_p ;
1502    for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1503      for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1504        *w0_p = *w1_p ;
1505        w0_p++ ;
1506      }
1507      w1_p++ ;
1508    }
1509    w.putStorage( w_p, b0 ) ;
1510    tint.freeStorage( ti_p, b1 ) ;
1511  }
1512  else if ( wtype_.compare( "TSYS" ) == 0 ) {
1513    Bool b0 ;
1514    Float *w_p = w.getStorage( b0 ) ;
1515    Float *w0_p = w_p ;
1516    for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1517      for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1518        Float temp = *w0_p ;
1519        *w0_p = 1.0 / ( temp * temp ) ;
1520        w0_p++ ;
1521      }
1522    }
1523    w.putStorage( w_p, b0 ) ;
1524  }
1525  else if ( wtype_.compare( "TINTSYS" ) == 0 ) {
1526    Bool b0, b1 ;
1527    Float *w_p = w.getStorage( b0 ) ;
1528    Float *w0_p = w_p ;
1529    const Double *ti_p = tint.getStorage( b1 ) ;
1530    const Double *w1_p = ti_p ;
1531    for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1532      Float interval = *w1_p ;
1533      for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1534        Float temp = *w0_p ;
1535        *w0_p = interval / ( temp * temp ) ;
1536        w0_p++ ;
1537      }
1538      w1_p++ ;
1539    }
1540    w.putStorage( w_p, b0 ) ;
1541    tint.freeStorage( ti_p, b1 ) ;
1542  }
1543  else {
1544    //LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
1545    //os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
1546    w = 1.0 ;
1547  }
1548}
1549
1550void STGrid::toInt( Array<uChar> &u, Array<Int> &v )
1551{
1552  uInt len = u.nelements() ;
1553  Int *int_p = new Int[len] ;
1554  Bool deleteIt ;
1555  const uChar *data_p = u.getStorage( deleteIt ) ;
1556  Int *i_p = int_p ;
1557  const uChar *u_p = data_p ;
1558  for ( uInt i = 0 ; i < len ; i++ ) {
1559    *i_p = ( *u_p == 0 ) ? 0 : 1 ;
1560    i_p++ ;
1561    u_p++ ;
1562  }
1563  u.freeStorage( data_p, deleteIt ) ;
1564  v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
1565}
1566
1567void STGrid::toInt( Array<uInt> &u, Array<Int> &v )
1568{
1569  uInt len = u.nelements() ;
1570  Int *int_p = new Int[len] ;
1571  Bool deleteIt ;
1572  const uInt *data_p = u.getStorage( deleteIt ) ;
1573  Int *i_p = int_p ;
1574  const uInt *u_p = data_p ;
1575  for ( uInt i = 0 ; i < len ; i++ ) {
1576    *i_p = ( *u_p == 0 ) ? 0 : 1 ;
1577    i_p++ ;
1578    u_p++ ;
1579  }
1580  u.freeStorage( data_p, deleteIt ) ;
1581  v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
1582}
1583
1584void STGrid::toPixel( Array<Double> &world, Array<Double> &pixel )
1585{
1586  uInt nrow = world.shape()[1] ;
1587  Bool bw, bp ;
1588  Double *w_p = world.getStorage( bw ) ;
1589  Double *p_p = pixel.getStorage( bp ) ;
1590  Double *ww_p = w_p ;
1591  Double *wp_p = p_p ;
1592  IPosition vshape( 1, 2 ) ;
1593  Vector<Double> _world, _pixel ;
1594  for ( uInt i = 0 ; i < nrow ; i++ ) {
1595    _world.takeStorage( vshape, ww_p, SHARE ) ;
1596    _pixel.takeStorage( vshape, wp_p, SHARE ) ;
1597    dircoord_->toPixel( _pixel, _world ) ;
1598    ww_p += 2 ;
1599    wp_p += 2 ;
1600  }
1601  world.putStorage( w_p, bw ) ;
1602  pixel.putStorage( p_p, bp ) ;
1603}
1604
1605void STGrid::boxFunc( Vector<Float> &convFunc, Int &convSize )
1606{
1607  convFunc = 0.0 ;
1608  for ( Int i = 0 ; i < convSize/2 ; i++ )
1609    convFunc(i) = 1.0 ;
1610}
1611
1612#define NEED_UNDERSCORES
1613#if defined(NEED_UNDERSCORES)
1614#define grdsf grdsf_
1615#define grdgauss grdgauss_
1616#define grdjinc1 grdjinc1_
1617#endif
1618#if defined(USE_CASAPY)
1619extern "C" {
1620  void grdsf(Double*, Double*);
1621  void grdgauss(Double*, Double*, Double*);
1622  void grdjinc1(Double*, Double*, Int*, Double*);
1623}
1624#else
1625extern "C" {
1626  void grdsf(Double*, Double*);
1627}
1628void grdgauss(Double *hwhm, Double *val, Double *out)
1629{
1630  *out = exp(-log(2.0) * (*val / *hwhm) * (*val / *hwhm));
1631}
1632void grdjinc1(Double *c, Double *val, Int *normalize, Double *out)
1633{
1634  // Calculate J_1(x) using approximate formula
1635  Double x = C::pi * *val / *c;
1636  Double ax = fabs(x);
1637  Double ans;
1638  if ( ax < 8.0 ) {
1639    Double y = x * x;
1640    Double ans1 = x * (72362614232.0 + y * (-7895059235.0
1641                       + y * (242396853.1 + y * (-2972611.439
1642                       + y * (15704.48260 + y * (-30.16036606))))));
1643    Double ans2 = 144725228442.0 + y * (2300535178.0
1644                       + y * (18583304.74 + y * (99447.43394
1645                       + y * (376.9991397 + y * 1.0))));
1646    ans = ans1 / ans2;
1647  }
1648  else {
1649    Double z = 8.0 / ax;
1650    Double y = z * z;
1651    Double xx = ax - 2.356194491;
1652    Double ans1 = 1.0 + y * (0.183105e-2 + y * (-0.3516396496e-4
1653                      + y * (0.2457520174e-5 + y * (-0.240337019e-6))));
1654    Double ans2 = 0.04687499995 + y * (-0.2002690873e-3
1655                      + y * (0.8449199096e-5 + y * (-0.88228987e-6
1656                      + y * (0.105787412e-6))));
1657    ans = sqrt(0.636619772 / ax) * (cos(xx) * ans1
1658                                    - z * sin(xx) * ans2);
1659    if (x < 0.0)
1660      ans = -ans;
1661  }
1662 
1663  // Then, calculate Jinc
1664  if (x == 0.0) {
1665    *out = 0.5;
1666  }
1667  else {
1668    *out = ans / x;
1669  }
1670
1671  if (*normalize == 1)
1672    *out = *out / 0.5; 
1673}
1674#endif
1675void STGrid::spheroidalFunc( Vector<Float> &convFunc )
1676{
1677  convFunc = 0.0 ;
1678  for ( Int i = 0 ; i < convSampling_*convSupport_ ; i++ ) {
1679    Double nu = Double(i) / Double(convSupport_*convSampling_) ;
1680    Double val ;
1681    grdsf( &nu, &val ) ;
1682    convFunc(i) = ( 1.0 - nu * nu ) * val ;
1683  }
1684}
1685
1686void STGrid::gaussFunc( Vector<Float> &convFunc, Double hwhm, Double truncate )
1687{
1688  convFunc = 0.0 ;
1689  Int len = (Int)(truncate*Double(convSampling_)+0.5);
1690  Double out, val;
1691  for ( Int i = 0 ; i < len ; i++ ) {
1692    val = Double(i) / Double(convSampling_) ;
1693    grdgauss(&hwhm, &val, &out);
1694    convFunc(i) = out;
1695  }
1696}
1697
1698void STGrid::gjincFunc( Vector<Float> &convFunc, Double hwhm, Double c, Double truncate )
1699{
1700  convFunc = 0.0;
1701  Double out1, out2, val;
1702  Int normalize = 1;
1703  if  (truncate >= 0.0) {
1704    Int len = (Int)(truncate*Double(convSampling_)+0.5);
1705    for (Int i = 0 ; i < len ; i++) {
1706      val = Double(i) / Double(convSampling_);
1707      grdgauss(&hwhm, &val, &out1);
1708      grdjinc1(&c, &val, &normalize, &out2);
1709      convFunc(i) = out1 * out2;
1710    }
1711  }
1712  else {
1713    Int len = convFunc.nelements();
1714    for (Int i = 0 ; i < len ; i++) {
1715      val = Double(i) / Double(convSampling_);
1716      grdjinc1(&c, &val, &normalize, &out2);
1717      if (out2 <= 0.0) {
1718        LogIO os(LogOrigin("STGrid","gjincFunc",WHERE));
1719        os << LogIO::DEBUG1 << "convFunc is automatically truncated at radius " << val << LogIO::POST;
1720        break;
1721      }
1722      grdgauss(&hwhm, &val, &out1);
1723      convFunc(i) = out1 * out2;
1724    }
1725  }
1726}
1727
1728void STGrid::pbFunc( Vector<Float> &convFunc )
1729{
1730  convFunc = 0.0 ;
1731}
1732
1733vector<float> STGrid::getConvFunc()
1734{
1735  LogIO os(LogOrigin("STGrid","getConvFunc",WHERE));
1736  Vector<Float> convFunc;
1737  vector<float> out;
1738
1739  if (cellx_ <= 0.0 || celly_ <= 0.0) {
1740    selectData();
1741    setupGrid();
1742  }
1743
1744  if (convType_ == "BOX" || convType_ == "SF") {
1745    setConvFunc(convFunc);
1746  }
1747  else if (convType_ == "GAUSS") {
1748    Quantum<Double> q1,q2;
1749    readQuantity(q1,gwidth_);
1750    readQuantity(q2,truncate_);
1751//     if (celly_ <= 0.0
1752//         && ((!q1.getUnit().empty()&&q1.getUnit()!="pixel") ||
1753//             (!q2.getUnit().empty()&&q2.getUnit()!="pixel"))) {
1754//       throw AipsError("You have to call defineImage to get correct convFunc");
1755//     }
1756    setConvFunc(convFunc);
1757  }
1758  else if (convType_ == "GJINC") {
1759    Quantum<Double> q1,q2,q3;
1760    readQuantity(q1,gwidth_);
1761    readQuantity(q2,truncate_);
1762    readQuantity(q3,jwidth_);
1763//     if (celly_ <= 0.0
1764//         && ((!q1.getUnit().empty()&&q1.getUnit()!="pixel") ||
1765//             (!q2.getUnit().empty()&&q2.getUnit()!="pixel") ||
1766//             (!q3.getUnit().empty()&&q3.getUnit()!="pixel"))) {
1767//       throw AipsError("You have to call defineImage to get correct convFunc");
1768//     }
1769    setConvFunc(convFunc);
1770  }
1771  else if (convType_ == "PB") {
1772    throw AipsError("Grid function PB is not available");
1773  }
1774  else {
1775    throw AipsError("Unknown grid function: "+convType_);
1776  }
1777
1778  convFunc.tovector(out);
1779  return out;
1780}
1781
1782void STGrid::setConvFunc( Vector<Float> &convFunc )
1783{
1784  LogIO os(LogOrigin("STGrid","setConvFunc",WHERE));
1785  convSupport_ = userSupport_ ;
1786  if ( convType_ == "BOX" ) {
1787    if ( convSupport_ < 0 )
1788      convSupport_ = 0 ;
1789    Int convSize = convSampling_ * ( 2 * convSupport_ + 2 )  ;
1790    convFunc.resize( convSize ) ;
1791    boxFunc( convFunc, convSize ) ;
1792    os << LogIO::DEBUGGING
1793       << "convType_ = " << convType_ << endl
1794       << "convSupport_ = " << convSupport_ << LogIO::POST;
1795  }
1796  else if ( convType_ == "SF" ) {
1797    if ( convSupport_ < 0 )
1798      convSupport_ = 3 ;
1799    Int convSize = convSampling_ * ( 2 * convSupport_ + 2 )  ;
1800    convFunc.resize( convSize ) ;
1801    spheroidalFunc( convFunc ) ;
1802    os << LogIO::DEBUGGING
1803       << "convType_ = " << convType_ << endl
1804       << "convSupport_ = " << convSupport_ << LogIO::POST;
1805  }
1806  else if ( convType_ == "GAUSS" ) {
1807    // determine pixel gwidth
1808    // default is HWHM corresponding to b = 1.0 (Mangum et al. 2007)
1809    Double pixelGW = -1.0;
1810    Quantum<Double> q ;
1811    if (!gwidth_.empty()) {
1812      readQuantity( q, gwidth_ );
1813      if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
1814        pixelGW = q.getValue();
1815      }
1816      else {
1817        pixelGW = q.getValue("rad")/celly_;
1818      }
1819    }
1820    pixelGW = (pixelGW >= 0.0) ? pixelGW : sqrt(log(2.0));
1821    if (pixelGW < 0.0) {
1822      os << LogIO::SEVERE
1823         << "Negative width is specified for gaussian" << LogIO::EXCEPTION;
1824    }
1825    // determine truncation radius
1826    // default is 3 * HWHM
1827    Double truncate = -1.0;
1828    if (!truncate_.empty()) {
1829      readQuantity( q, truncate_ );
1830      if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
1831        truncate = q.getValue();
1832      }
1833      else {
1834        truncate = q.getValue("rad")/celly_;
1835      }
1836    }     
1837    //convSupport_ = (Int)(truncate+0.5);
1838    truncate = (truncate >= 0.0) ? truncate : 3.0 * pixelGW;
1839    convSupport_ = Int(truncate);
1840    convSupport_ += (((truncate-(Double)convSupport_) > 0.0) ? 1 : 0);
1841    Int convSize = convSampling_ * ( 2*convSupport_ + 2 ) ;
1842    convFunc.resize( convSize ) ;
1843    gaussFunc( convFunc, pixelGW, truncate ) ;
1844    os << LogIO::DEBUGGING
1845       << "convType_ = " << convType_ << endl
1846       << "convSupport_ = " << convSupport_ << endl
1847       << "truncate_ = " << truncate << "pixel" << endl
1848       << "gwidth_ = " << pixelGW << "pixel" << LogIO::POST;
1849  }
1850  else if ( convType_ == "GJINC" ) {
1851    // determine pixel gwidth
1852    // default is HWHM corresponding to b = 2.52 (Mangum et al. 2007)
1853    Double pixelGW = -1.0;
1854    Quantum<Double> q ;
1855    if (!gwidth_.empty()) {
1856      readQuantity( q, gwidth_ );
1857      if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
1858        pixelGW = q.getValue();
1859      }
1860      else {
1861        pixelGW = q.getValue("rad")/celly_;
1862      }
1863    }
1864    pixelGW = (pixelGW >= 0.0) ? pixelGW : sqrt(log(2.0)) * 2.52;
1865    if (pixelGW < 0.0) {
1866      os << LogIO::SEVERE
1867         << "Negative width is specified for gaussian" << LogIO::EXCEPTION;
1868    }
1869    // determine pixel c
1870    // default is c = 1.55 (Mangum et al. 2007)
1871    Double pixelJW = -1.0;
1872    if (!jwidth_.empty()) {
1873      readQuantity( q, jwidth_ );
1874      if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
1875        pixelJW = q.getValue();
1876      }
1877      else {
1878        pixelJW = q.getValue("rad")/celly_;
1879      }
1880    }
1881    pixelJW = (pixelJW >= 0.0) ? pixelJW : 1.55;
1882    if (pixelJW < 0.0) {
1883      os << LogIO::SEVERE
1884         << "Negative width is specified for jinc" << LogIO::EXCEPTION;
1885    }
1886    // determine truncation radius
1887    // default is -1.0 (truncate at first null)
1888    Double truncate = -1.0;
1889    if (!truncate_.empty()) {
1890      readQuantity( q, truncate_ );
1891      if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
1892        truncate = q.getValue();
1893      }
1894      else {
1895        truncate = q.getValue("rad")/celly_;
1896      }
1897    }     
1898    //convSupport_ = (truncate >= 0.0) ? (Int)(truncate+0.5) : (Int)(2*pixelJW+0.5);
1899    Double convSupportF = (truncate >= 0.0) ? truncate : (2*pixelJW);
1900    convSupport_ = (Int)convSupportF;
1901    convSupport_ += (((convSupportF-(Double)convSupport_) > 0.0) ? 1 : 0);
1902    Int convSize = convSampling_ * ( 2*convSupport_ + 2 ) ;
1903    convFunc.resize( convSize ) ;
1904    gjincFunc( convFunc, pixelGW, pixelJW, truncate ) ;
1905    os << LogIO::DEBUGGING
1906       << "convType_ = " << convType_ << endl
1907       << "convSupport_ = " << convSupport_ << endl
1908       << "truncate_ = " << truncate << "pixel" << endl
1909       << "gwidth_ = " << pixelGW << "pixel" << endl
1910       << "jwidth_ = " << pixelJW << "pixel" << LogIO::POST;
1911  }
1912  else if ( convType_ == "PB" ) {
1913    if ( convSupport_ < 0 )
1914      convSupport_ = 0 ;
1915    pbFunc( convFunc ) ;
1916  }
1917  else {
1918    throw AipsError( "Unsupported convolution function" ) ;
1919  }
1920}
1921
1922string STGrid::saveData( string outfile )
1923{
1924  LogIO os( LogOrigin("STGrid", "saveData", WHERE) ) ;
1925  double t0, t1 ;
1926  t0 = mathutil::gettimeofday_sec() ;
1927
1928  //Int polno = 0 ;
1929  String outfile_ ;
1930  if ( outfile.size() == 0 ) {
1931    if ( infileList_[0].lastchar() == '/' ) {
1932      outfile_ = infileList_[0].substr( 0, infileList_[0].size()-1 ) ;
1933    }
1934    else {
1935      outfile_ = infileList_[0] ;
1936    }
1937    outfile_ += ".grid" ;
1938  }
1939  else {
1940    outfile_ = outfile ;
1941  }
1942  Table tab ;
1943  prepareTable( tab, outfile_ ) ;
1944  fillTable( tab ) ;
1945
1946  t1 = mathutil::gettimeofday_sec() ;
1947  os << LogIO::DEBUGGING << "saveData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
1948
1949  return outfile_ ;
1950}
1951
1952void STGrid::prepareTable( Table &tab, String &name )
1953{
1954  Table t( infileList_[0], Table::Old ) ;
1955  t.deepCopy( name, Table::New, False, t.endianFormat(), True ) ;
1956  tab = Table( name, Table::Update ) ;
1957  // 2012/02/13 TN
1958  // explicitly copy subtables since no rows including subtables are
1959  // copied by Table::deepCopy with noRows=True
1960  //TableCopy::copySubTables( tab, t ) ;
1961  const TableRecord &inrec = t.keywordSet();
1962  TableRecord &outrec = tab.rwKeywordSet();
1963  for (uInt i = 0 ; i < inrec.nfields() ; i++) {
1964    if (inrec.type(i) == TpTable) {
1965      String name = inrec.name(i);
1966      Table intable = inrec.asTable(name);
1967      Table outtable = outrec.asTable(name);
1968      TableCopy::copyRows(outtable, intable);
1969    }
1970  }
1971}
1972
1973void STGrid::fillTable( Table &tab )
1974{
1975  //IPosition dshape = data_.shape() ;
1976  Int nrow = nx_ * ny_ * npol_ ;
1977  tab.rwKeywordSet().define( "nPol", npol_ ) ;
1978  tab.addRow( nrow ) ;
1979  Vector<Double> cpix( 2 ) ;
1980  cpix(0) = Double( nx_ - 1 ) * 0.5 ;
1981  cpix(1) = Double( ny_ - 1 ) * 0.5 ;
1982  Vector<Double> dir( 2 ) ;
1983  Vector<Double> pix( 2 );
1984  ArrayColumn<Double> directionCol( tab, "DIRECTION" ) ;
1985  ArrayColumn<Float> spectraCol( tab, "SPECTRA" ) ;
1986  ArrayColumn<uChar> flagtraCol( tab, "FLAGTRA" ) ;
1987  ScalarColumn<uInt> flagRowCol( tab, "FLAGROW" );
1988  ScalarColumn<uInt> polnoCol( tab, "POLNO" ) ;
1989  ScalarColumn<uInt> scannoCol( tab, "SCANNO" ) ;
1990  Int irow = 0 ;
1991  Vector<Float> sp( nchan_ ) ;
1992  Vector<uChar> flag( nchan_ ) ;
1993  Bool bsp, bdata, bflag ;
1994  const Float *data_p = data_.getStorage( bdata ) ;
1995  const uChar *flag_p = flag_.getStorage( bflag ) ;
1996  Float *wsp_p, *sp_p ;
1997  const Float *wdata_p = data_p ;
1998  const uChar *wflag_p = flag_p ;
1999  long step = nx_ * ny_ * npol_ ;
2000  long offset ;
2001  uInt scanno = 0 ;
2002  uChar rflag;
2003  for ( Int iy = 0 ; iy < ny_ ; iy++ ) {
2004    pix(1) = (Double)(iy);
2005    for ( Int ix = 0 ; ix < nx_ ; ix++ ) {
2006      pix(0) = (Double)(ix);
2007      dircoord_->toWorld(dir,pix);
2008      //os << "dir[" << ix << "," << iy << "]=" << dir << LogIO::POST;
2009      for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
2010        offset = ix + nx_ * (iy + ipol * ny_) ;
2011        //os << "offset = " << offset << LogIO::POST ;
2012        sp_p = sp.getStorage( bsp ) ;
2013        wsp_p = sp_p ;
2014        wdata_p = data_p + offset ;
2015        wflag_p = flag_p + offset ;
2016        rflag = ~0 ; //11111111
2017        for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
2018          *wsp_p = *wdata_p ;
2019          wsp_p++ ;
2020          wdata_p += step ;
2021          flag[ichan] = *wflag_p ;
2022          rflag &= flag[ichan] ;
2023          wflag_p += step ;
2024        }
2025        sp.putStorage( sp_p, bsp ) ;
2026        spectraCol.put( irow, sp ) ;
2027        flagtraCol.put( irow, flag ) ;
2028        flagRowCol.put( irow, ((rflag>0) ? 1 : 0) ) ;
2029        directionCol.put( irow, dir ) ;
2030        polnoCol.put( irow, pollist_[ipol] ) ;
2031        scannoCol.put( irow, scanno ) ;
2032        irow++ ;
2033      }
2034      scanno++ ;
2035    }
2036  }
2037  data_.freeStorage( data_p, bdata ) ;
2038  flag_.freeStorage( flag_p, bflag ) ;
2039
2040  fillMainColumns( tab ) ;
2041}
2042
2043void STGrid::fillMainColumns( Table &tab )
2044{
2045  // values for fill
2046  //Table t( infileList_[0], Table::Old ) ;
2047  Table t ;
2048  table( t, 0 ) ;
2049  Table tsel = t( t.col( "IFNO" ) == (uInt)ifno_, 1 ) ;
2050  ROTableRow row( tsel ) ;
2051  row.get( 0 ) ;
2052  const TableRecord &rec = row.record() ;
2053  uInt freqId = rec.asuInt( "FREQ_ID" ) ;
2054  uInt molId = rec.asuInt( "MOLECULE_ID" ) ;
2055  uInt tcalId = rec.asuInt( "TCAL_ID" ) ;
2056  uInt focusId = rec.asuInt( "FOCUS_ID" ) ;
2057  uInt weatherId = rec.asuInt( "WEATHER_ID" ) ;
2058  String srcname = rec.asString( "SRCNAME" ) ;
2059  String fieldname = rec.asString( "FIELDNAME" ) ;
2060  Vector<Float> defaultTsys( 1, 1.0 ) ;
2061  // @todo how to set flagtra for gridded spectra?
2062  Vector<uChar> flagtra = rec.asArrayuChar( "FLAGTRA" ) ;
2063  flagtra = (uChar)0 ;
2064  Float opacity = rec.asFloat( "OPACITY" ) ;
2065  Double srcvel = rec.asDouble( "SRCVELOCITY" ) ;
2066  Vector<Double> srcpm = rec.asArrayDouble( "SRCPROPERMOTION" ) ;
2067  Vector<Double> srcdir = rec.asArrayDouble( "SRCDIRECTION" ) ;
2068  Vector<Double> scanrate = rec.asArrayDouble( "SCANRATE" ) ;
2069  Double time = rec.asDouble( "TIME" ) ;
2070  Double interval = rec.asDouble( "INTERVAL" ) ;
2071
2072  // fill columns
2073  Int nrow = tab.nrow() ;
2074  ScalarColumn<uInt> ifnoCol( tab, "IFNO" ) ;
2075  ScalarColumn<uInt> beamnoCol(tab, "BEAMNO");
2076  ScalarColumn<uInt> freqIdCol( tab, "FREQ_ID" ) ;
2077  ScalarColumn<uInt> molIdCol( tab, "MOLECULE_ID" ) ;
2078  ScalarColumn<uInt> tcalidCol( tab, "TCAL_ID" ) ;
2079  ScalarColumn<Int> fitidCol( tab, "FIT_ID" ) ;
2080  ScalarColumn<uInt> focusidCol( tab, "FOCUS_ID" ) ;
2081  ScalarColumn<uInt> weatheridCol( tab, "WEATHER_ID" ) ;
2082  ArrayColumn<uChar> flagtraCol( tab, "FLAGTRA" ) ;
2083  ScalarColumn<uInt> rflagCol( tab, "FLAGROW" ) ;
2084  ArrayColumn<Float> tsysCol( tab, "TSYS" ) ;
2085  ScalarColumn<String> srcnameCol( tab, "SRCNAME" ) ;
2086  ScalarColumn<String> fieldnameCol( tab, "FIELDNAME" ) ;
2087  ScalarColumn<Int> srctypeCol( tab, "SRCTYPE" ) ;
2088  ScalarColumn<Float> opacityCol( tab, "OPACITY" ) ;
2089  ScalarColumn<Double> srcvelCol( tab, "SRCVELOCITY" ) ;
2090  ArrayColumn<Double> srcpmCol( tab, "SRCPROPERMOTION" ) ;
2091  ArrayColumn<Double> srcdirCol( tab, "SRCDIRECTION" ) ;
2092  ArrayColumn<Double> scanrateCol( tab, "SCANRATE" ) ;
2093  ScalarColumn<Double> timeCol( tab, "TIME" ) ;
2094  ScalarColumn<Double> intervalCol( tab, "INTERVAL" ) ;
2095  for ( Int i = 0 ; i < nrow ; i++ ) {
2096    ifnoCol.put( i, (uInt)ifno_ ) ;
2097    beamnoCol.put(i, 0);
2098    freqIdCol.put( i, freqId ) ;
2099    molIdCol.put( i, molId ) ;
2100    tcalidCol.put( i, tcalId ) ;
2101    fitidCol.put( i, -1 ) ;
2102    focusidCol.put( i, focusId ) ;
2103    weatheridCol.put( i, weatherId ) ;
2104    //flagtraCol.put( i, flagtra ) ;
2105    //rflagCol.put( i, 0 ) ;
2106    tsysCol.put( i, defaultTsys ) ;
2107    srcnameCol.put( i, srcname ) ;
2108    fieldnameCol.put( i, fieldname ) ;
2109    srctypeCol.put( i, (Int)SrcType::PSON ) ;
2110    opacityCol.put( i, opacity ) ;
2111    srcvelCol.put( i, srcvel ) ;
2112    srcpmCol.put( i, srcpm ) ;
2113    srcdirCol.put( i, srcdir ) ;
2114    scanrateCol.put( i, scanrate ) ;
2115    timeCol.put( i, time ) ;
2116    intervalCol.put( i, interval ) ;
2117    if ( (i + 1) % npol_ == 0 ) {
2118      time += interval / 86400.0;
2119    }
2120  }
2121}
2122
2123vector<int> STGrid::getResultantMapSize()
2124{
2125  vector<int> r(2);
2126  r[0] = nx_;
2127  r[1] = ny_;
2128  return r;
2129}
2130
2131vector<double> STGrid::getResultantCellSize()
2132{
2133  vector<double> r(2);
2134  r[0] = cellx_;
2135  r[1] = celly_;
2136  return r;
2137}
2138
2139// STGrid2
2140STGrid2::STGrid2()
2141  : STGrid()
2142{
2143}
2144
2145STGrid2::STGrid2( const ScantableWrapper &s )
2146  : STGrid()
2147{
2148  setScantable( s ) ;
2149}
2150
2151STGrid2::STGrid2( const vector<ScantableWrapper> &v )
2152  : STGrid()
2153{
2154  setScantableList( v ) ;
2155}
2156
2157void STGrid2::setScantable( const ScantableWrapper &s )
2158{
2159  nfile_ = 1 ;
2160  dataList_.resize( nfile_ ) ;
2161  dataList_[0] = s ;
2162  infileList_.resize( nfile_ ) ;
2163  infileList_[0] = s.getCP()->table().tableName() ;
2164}
2165
2166void STGrid2::setScantableList( const vector<ScantableWrapper> &v )
2167{
2168  nfile_ = v.size() ;
2169  dataList_.resize( nfile_ ) ;
2170  infileList_.resize( nfile_ ) ;
2171  for ( uInt i = 0 ; i < nfile_ ; i++ ) {
2172    dataList_[i] = v[i] ;
2173    infileList_[i] = v[i].getCP()->table().tableName() ;
2174  }
2175}
2176
2177ScantableWrapper STGrid2::getResultAsScantable( int tp )
2178{
2179  ScantableWrapper sw( tp ) ;
2180  CountedPtr<Scantable> s = sw.getCP() ;
2181  s->setHeader( dataList_[0].getCP()->getHeader() ) ;
2182  Table tout, tin ;
2183  String subt[] = { "FREQUENCIES", "FOCUS", "WEATHER",
2184                    "TCAL", "MOLECULES", "HISTORY", "FIT" } ;
2185  for ( uInt i = 0 ; i < 7 ; i++ ) {
2186    tout = s->table().rwKeywordSet().asTable(subt[i]) ;
2187    tin = dataList_[0].getCP()->table().rwKeywordSet().asTable(subt[i]) ;
2188    TableCopy::copyRows( tout, tin ) ;
2189    tout.rwKeywordSet() = tin.rwKeywordSet();
2190  }
2191  fillTable( s->table() ) ;
2192  return sw ;
2193}
2194
2195void STGrid2::table( Table &tab, uInt i )
2196{
2197  if ( i < nfile_ )
2198    tab = dataList_[i].getCP()->table() ;
2199}
2200
2201}
Note: See TracBrowser for help on using the repository browser.