source: trunk/src/STGrid.cpp @ 2629

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

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: test_sdgrid

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Several fixes about R.A. distribution.


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