source: trunk/src/STGrid.cpp @ 2638

Last change on this file since 2638 was 2638, 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...

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