source: trunk/src/STGrid.cpp @ 2671

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

New Development: No

JIRA Issue: Yes CAS-4429

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

Implemented gauss and gjinc gridding. Added some arguments to
pass necessary parameters for those grid functions.

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