source: trunk/src/STGrid.cpp @ 3015

Last change on this file since 3015 was 3015, checked in by Takeshi Nakazato, 9 years ago

New Development: No

JIRA Issue: Yes CAS-4362

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...

Fix for the issue on sdgrid that resulting map extent expands or shrinks when cell size is calculated by the task.
A cause of the issue is that the task always takes 10% margin when it calculates cell size automatically. In this fix, margin is set to 1 pixel.
File size: 63.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 = " << userSupport_ << 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    : common(common), self(obj), 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    : common(common), self(obj), 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  flag_.resize( gshape ) ;
664  flag_ = (uChar)0;
665  STCommonData common = STCommonData(gshape, data_);
666  common.gnx = gnx ;
667  common.gny = gny ;
668
669  // parameters for gridding
670  Int *chanMap = new Int[nchan_] ;
671  for ( Int i = 0 ; i < nchan_ ; i++ ) {
672    chanMap[i] = i ;
673  }
674  common.chanMap = chanMap;
675
676  // convolution kernel
677  t0 = mathutil::gettimeofday_sec() ;
678  setConvFunc( common.convFunc ) ;
679  t1 = mathutil::gettimeofday_sec() ;
680  os << LogIO::DEBUGGING << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
681
682  // for performance check
683  eGetData_ = 0.0 ;
684  eToPixel_ = 0.0 ;
685  eGGridSD_ = 0.0 ;
686  double eInitPol = 0.0 ;
687
688  for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) {
689    initTable( ifile ) ;
690
691    os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ;   
692    Broker broker = Broker(produceChunk, consumeChunk);
693    for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
694      t0 = mathutil::gettimeofday_sec() ;
695      initPol( ipol ) ; // set ptab_ and attach()
696      t1 = mathutil::gettimeofday_sec() ;
697      eInitPol += t1-t0 ;
698     
699      STContext context(this, common, ipol);
700     
701      os << "start pol " << ipol << LogIO::POST ;
702     
703      nprocessed_ = 0 ;
704#if 1
705      broker.runProducerAsMasterThread(&context, DO_AHEAD);
706#else
707      for (;;) {
708        bool produced = produceChunk(&context);
709        if (! produced) {
710          break;
711        }
712        consumeChunk(&context);
713      }
714#endif
715
716      os << "end pol " << ipol << LogIO::POST ;
717
718    }
719    os << "end table " << ifile << LogIO::POST ;   
720  }
721  os << LogIO::DEBUGGING << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ;
722  os << LogIO::DEBUGGING << "getData: elapsed time is " << eGetData_-eToInt-eGetWeight << " sec." << LogIO::POST ;
723  os << LogIO::DEBUGGING << "toPixel: elapsed time is " << eToPixel_ << " sec." << LogIO::POST ;
724  os << LogIO::DEBUGGING << "ggridsd: elapsed time is " << eGGridSD_ << " sec." << LogIO::POST ;
725  os << LogIO::DEBUGGING << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
726  os << LogIO::DEBUGGING << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
727 
728  delete chanMap ;
729
730  // set data
731  setData( common.gdataArrC, common.gwgtArr ) ;
732
733}
734
735void STGrid::consumeChunkWithClipping(void *ctx) throw(PCException)
736{
737  STContextWithClipping &context = *(STContextWithClipping *)ctx;
738  STGChunk *chunk = NULL;
739  try {
740    context.queue.lock();
741    chunk = context.queue.get();
742    context.queue.unlock();
743  } catch (FullException &e) {
744    context.queue.unlock();
745    // TODO: log error
746    throw PCException();
747  }
748
749  double t0, t1 ;
750  // world -> pixel
751  Array<Double> xypos( context.self->dshape_ ) ;
752  t0 = mathutil::gettimeofday_sec() ;
753  context.self->toPixel( chunk->direction, xypos ) ;
754  t1 = mathutil::gettimeofday_sec() ;
755  context.self->eToPixel_ += t1-t0 ;
756   
757  // call ggridsd
758  Int nvispol = 1 ;
759  Int irow = -1 ;
760  t0 = mathutil::gettimeofday_sec() ;
761  context.self->call_ggridsd2( xypos,
762                chunk->spectra,
763                nvispol,
764                context.self->nchan_,
765                chunk->flagtra,
766                chunk->rflag,
767                chunk->weight,
768                chunk->nrow,
769                irow,
770                context.common.gdataArrC,
771                context.common.gwgtArr,
772                context.common.npoints,
773                context.common.clipMin,
774                context.common.clipWMin,
775                context.common.clipCMin,
776                context.common.clipMax,
777                context.common.clipWMax,
778                context.common.clipCMax,
779                context.common.gnx,
780                context.common.gny,
781                context.self->npol_,
782                context.self->nchan_,
783                context.self->convSupport_,
784                context.self->convSampling_,
785                context.common.convFunc,
786                context.common.chanMap,
787                (Int*)&context.pol ) ;
788  t1 = mathutil::gettimeofday_sec() ;
789  context.self->eGGridSD_ += t1-t0 ;
790 
791  delete chunk;
792}
793
794void STGrid::gridPerRowWithClipping()
795{
796  LogIO os( LogOrigin("STGrid", "gridPerRowWithClipping", WHERE) ) ;
797  double t0, t1 ;
798
799
800  // grid data
801  // Extend grid plane with convSupport_
802  //   Int gnx = nx_+convSupport_*2 ;
803  //   Int gny = ny_+convSupport_*2 ;
804  Int gnx = nx_;
805  Int gny = ny_;
806
807  IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
808  IPosition pshape( 3, gnx, gny, npol_ ) ;
809  // 2011/12/20 TN
810  // data_ and gwgtArr share storage
811  data_.resize( gshape ) ;
812  data_ = 0.0 ;
813  flag_.resize( gshape ) ;
814  flag_ = (uChar)0;
815  STCommonDataWithClipping common = STCommonDataWithClipping( gshape,
816                                                              pshape,
817                                                              data_ ) ;
818  common.gnx = gnx ;
819  common.gny = gny ;
820
821  // parameters for gridding
822  Int *chanMap = new Int[nchan_] ;
823  for ( Int i = 0 ; i < nchan_ ; i++ ) {
824    chanMap[i] = i ;
825  }
826  common.chanMap = chanMap;
827
828  // convolution kernel
829  t0 = mathutil::gettimeofday_sec() ;
830  setConvFunc( common.convFunc ) ;
831  t1 = mathutil::gettimeofday_sec() ;
832  os << LogIO::DEBUGGING << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
833
834  // for performance check
835  eGetData_ = 0.0 ;
836  eToPixel_ = 0.0 ;
837  eGGridSD_ = 0.0 ;
838  double eInitPol = 0.0 ;
839
840  for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) {
841    initTable( ifile ) ;
842
843    os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ;   
844    Broker broker = Broker(produceChunk, consumeChunkWithClipping);
845    for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
846      t0 = mathutil::gettimeofday_sec() ;
847      initPol( ipol ) ; // set ptab_ and attach()
848      t1 = mathutil::gettimeofday_sec() ;
849      eInitPol += t1-t0 ;
850     
851      STContextWithClipping context(this, common, ipol);
852     
853      os << "start pol " << ipol << LogIO::POST ;
854     
855      nprocessed_ = 0 ;
856#if 1
857      broker.runProducerAsMasterThread(&context, DO_AHEAD);
858#else
859      for (;;) {
860        bool produced = produceChunk(&context);
861        if (! produced) {
862          break;
863        }
864        consumeChunkWithClipping(&context);
865      }
866#endif
867
868      os << "end pol " << ipol << LogIO::POST ;
869
870    }
871    os << "end table " << ifile << LogIO::POST ;   
872  }
873  os << LogIO::DEBUGGING << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ;
874  os << LogIO::DEBUGGING << "getData: elapsed time is " << eGetData_-eToInt-eGetWeight << " sec." << LogIO::POST ;
875  os << LogIO::DEBUGGING << "toPixel: elapsed time is " << eToPixel_ << " sec." << LogIO::POST ;
876  os << LogIO::DEBUGGING << "ggridsd2: elapsed time is " << eGGridSD_ << " sec." << LogIO::POST ;
877  os << LogIO::DEBUGGING << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
878  os << LogIO::DEBUGGING << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
879 
880  delete chanMap ;
881
882  // clip min and max in each grid
883//   os << "BEFORE CLIPPING" << LogIO::POST ;
884//   os << "gdataArrC=" << common.gdataArrC << LogIO::POST ;
885//   os << "gwgtArr=" << common.gwgtArr << LogIO::POST ;
886  t0 = mathutil::gettimeofday_sec() ;
887  clipMinMax( common.gdataArrC,
888              common.gwgtArr,
889              common.npoints,
890              common.clipMin,
891              common.clipWMin,
892              common.clipCMin,
893              common.clipMax,
894              common.clipWMax,
895              common.clipCMax ) ;
896  t1 = mathutil::gettimeofday_sec() ;
897  os << LogIO::DEBUGGING << "clipMinMax: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
898//   os << "AFTER CLIPPING" << LogIO::POST ;
899//   os << "gdataArrC=" << common.gdataArrC << LogIO::POST ;
900//   os << "gwgtArr=" << common.gwgtArr << LogIO::POST ;
901
902  // set data
903  setData( common.gdataArrC, common.gwgtArr ) ;
904
905}
906
907void STGrid::clipMinMax( Array<Complex> &grid,
908                         Array<Float> &weight,
909                         Array<Int> &npoints,
910                         Array<Complex> &clipmin,
911                         Array<Float> &clipwmin,
912                         Array<Float> &clipcmin,
913                         Array<Complex> &clipmax,
914                         Array<Float> &clipwmax,
915                         Array<Float> &clipcmax )
916{
917  //LogIO os( LogOrigin("STGrid","clipMinMax",WHERE) ) ;
918
919  // prepare pointers
920  Bool delG, delW, delNP, delCMin, delCWMin, delCCMin, delCMax, delCWMax, delCCMax ;
921  Complex *grid_p = grid.getStorage( delG ) ;
922  Float *wgt_p = weight.getStorage( delW ) ;
923  const Int *npts_p = npoints.getStorage( delNP ) ;
924  const Complex *cmin_p = clipmin.getStorage( delCMin ) ;
925  const Float *cwmin_p = clipwmin.getStorage( delCWMin ) ;
926  const Float *ccmin_p = clipcmin.getStorage( delCCMin ) ;
927  const Complex *cmax_p = clipmax.getStorage( delCMax ) ;
928  const Float *cwmax_p = clipwmax.getStorage( delCWMax ) ;
929  const Float *ccmax_p = clipcmax.getStorage( delCCMax ) ;
930
931  const IPosition &gshape = grid.shape() ;
932  long offset = gshape[0] * gshape[1] * gshape[2] ; // nx * ny * npol
933  Int nchan = gshape[3] ;
934  long origin = nchan * offset ;
935  for ( long i = 0 ; i < offset ; i++ ) {
936    if ( *npts_p > 2 ) {
937      for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
938        // clip minimum and maximum
939        *grid_p -= (*cmin_p)*(*cwmin_p)*(*ccmin_p)
940          + (*cmax_p)*(*cwmax_p)*(*ccmax_p) ;
941        *wgt_p -= (*cwmin_p)*(*ccmin_p)
942          + (*cwmax_p)*(*ccmax_p) ;
943       
944        grid_p += offset ;
945        wgt_p += offset ;
946        cmin_p += offset ;
947        cwmin_p += offset ;
948        ccmin_p += offset ;
949        cmax_p += offset ;
950        cwmax_p += offset ;
951        ccmax_p += offset ;
952      }
953      grid_p -= origin ;
954      wgt_p -= origin ;
955      cmin_p -= origin ;
956      cwmin_p -= origin ;
957      ccmin_p -= origin ;
958      cmax_p -= origin ;
959      cwmax_p -= origin ;
960      ccmax_p -= origin ;
961    }
962    grid_p++ ;
963    wgt_p++ ;
964    npts_p++ ;
965    cmin_p++ ;
966    cwmin_p++ ;
967    ccmin_p++ ;
968    cmax_p++ ;
969    cwmax_p++ ;
970    ccmax_p++ ;
971  }
972  grid_p -= offset ;
973  wgt_p -= offset ;
974  npts_p -= offset ;
975  cmin_p -= offset ;
976  cwmin_p -= offset ;
977  ccmin_p -= offset ;
978  cmax_p -= offset ;
979  cwmax_p -= offset ;
980  ccmax_p -= offset ; 
981
982  // finalization
983  grid.putStorage( grid_p, delG ) ;
984  weight.putStorage( wgt_p, delW ) ;
985  npoints.freeStorage( npts_p, delNP ) ;
986  clipmin.freeStorage( cmin_p, delCMin ) ;
987  clipwmin.freeStorage( cwmin_p, delCWMin ) ;
988  clipcmin.freeStorage( ccmin_p, delCCMin ) ;
989  clipmax.freeStorage( cmax_p, delCMax ) ;
990  clipwmax.freeStorage( cwmax_p, delCWMax ) ;
991  clipcmax.freeStorage( ccmax_p, delCCMax ) ;
992}
993
994void STGrid::initPol( Int ipol )
995{
996  LogIO os( LogOrigin("STGrid","initPol",WHERE) ) ;
997  if ( npolOrg_ == 1 ) {
998    os << "single polarization data." << LogIO::POST ;
999    ptab_ = tab_ ;
1000  }
1001  else
1002    ptab_ = tab_( tab_.col("POLNO") == pollist_[ipol] ) ;
1003
1004  attach( ptab_ ) ;
1005}
1006
1007void STGrid::initTable( uInt idx )
1008{
1009  tab_ = tableList_[idx] ;
1010  nrow_ = rows_[idx] ;
1011  updateChunkShape() ;
1012}
1013
1014void STGrid::setData( Array<Complex> &gdata,
1015                      Array<Float> &gwgt )
1016{
1017  // 2011/12/20 TN
1018  // gwgt and data_ share storage
1019  LogIO os( LogOrigin("STGrid","setData",WHERE) ) ;
1020  double t0, t1 ;
1021  t0 = mathutil::gettimeofday_sec() ;
1022  uInt len = data_.nelements() ;
1023  Bool b1, b2, b3 ;
1024  const Complex *gdata_p = gdata.getStorage( b1 ) ;
1025  Float *gwgt_p = gwgt.getStorage( b2 ) ; // storage shared with data_
1026  uChar *gflg_p = flag_.getStorage( b3 ) ;
1027  for ( uInt i = 0 ; i < len ; i++ ) {
1028    if (gwgt_p[i] > 0.0) {
1029      gwgt_p[i] = (gdata_p[i]).real() / gwgt_p[i];
1030      gflg_p[i] = (uChar)0;
1031    }
1032    else {
1033      gflg_p[i] = (uChar)1;
1034    }
1035  }
1036  gdata.freeStorage( gdata_p, b1 ) ;
1037  data_.putStorage( gwgt_p, b2 ) ;
1038  flag_.putStorage( gflg_p, b3 ) ;
1039  t1 = mathutil::gettimeofday_sec() ;
1040  os << LogIO::DEBUGGING << "setData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
1041}
1042
1043void STGrid::setupGrid()
1044{
1045  Double xmin,xmax,ymin,ymax ;
1046  mapExtent( xmin, xmax, ymin, ymax ) ;
1047 
1048  setupGrid( nxUI_, nyUI_, cellxUI_, cellyUI_,
1049             xmin, xmax, ymin, ymax, centerUI_ ) ;
1050}
1051
1052void STGrid::setupGrid( Int &nx,
1053                        Int &ny,
1054                        String &cellx,
1055                        String &celly,
1056                        Double &xmin,
1057                        Double &xmax,
1058                        Double &ymin,
1059                        Double &ymax,
1060                        String &center )
1061{
1062  LogIO os( LogOrigin("STGrid","setupGrid",WHERE) ) ;
1063  //cout << "nx=" << nx << ", ny=" << ny << endl ;
1064
1065  // center position
1066  if ( center.size() == 0 ) {
1067    center_(0) = 0.5 * ( xmin + xmax ) ;
1068    center_(1) = 0.5 * ( ymin + ymax ) ;
1069  }
1070  else {
1071    String::size_type pos0 = center.find( " " ) ;
1072    if ( pos0 == String::npos ) {
1073      throw AipsError( "bad string format in parameter center" ) ;
1074    }
1075    String::size_type pos1 = center.find( " ", pos0+1 ) ;
1076    String typestr, xstr, ystr ;
1077    if ( pos1 != String::npos ) {
1078      typestr = center.substr( 0, pos0 ) ;
1079      xstr = center.substr( pos0+1, pos1-pos0 ) ;
1080      ystr = center.substr( pos1+1 ) ;
1081      // todo: convert to J2000 (or direction ref for DIRECTION column)
1082    }
1083    else {
1084      typestr = "J2000" ;
1085      xstr = center.substr( 0, pos0 ) ;
1086      ystr = center.substr( pos0+1 ) ;
1087    }
1088    QuantumHolder qh ;
1089    String err ;
1090    qh.fromString( err, xstr ) ;
1091    Quantum<Double> xcen = qh.asQuantumDouble() ;
1092    qh.fromString( err, ystr ) ;
1093    Quantum<Double> ycen = qh.asQuantumDouble() ;
1094    center_(0) = xcen.getValue( "rad" ) ;
1095    center_(1) = ycen.getValue( "rad" ) ;
1096    double base = 0.5 * (xmin + xmax) ;
1097    int maxrotate = 1 ;
1098    int nelem = 2 * maxrotate + 1 ;
1099    double *sep = new double[nelem] ;
1100    for ( int i = 0 ; i < nelem ; i++ )
1101      sep[i] = abs(base - center_[0] - (i-maxrotate) * C::_2pi) ;
1102//     os << "sep[0]=" << sep[0] << endl 
1103//        << "sep[1]=" << sep[1] << endl
1104//        << "sep[2]=" << sep[2] << LogIO::POST ;
1105    int idx = 0 ;
1106    base = sep[0] ;
1107    int nrotate = 0 ;
1108    while ( idx < nelem ) {
1109      if ( base > sep[idx] ) {
1110        base = sep[idx] ;
1111        nrotate = idx ;
1112      }
1113      idx++ ;
1114    }
1115    delete sep ;
1116    nrotate -= maxrotate ;
1117//     os << "nrotate = " << nrotate << LogIO::POST ;
1118    center_[0] += nrotate * C::_2pi ;
1119  }
1120//   os << "xmin=" << xmin << LogIO::POST ;
1121//   os << "center_=" << center_ << LogIO::POST ;
1122
1123  nx_ = nx ;
1124  ny_ = ny ;
1125  if ( nx < 0 && ny > 0 ) {
1126    nx_ = ny ;
1127    ny_ = ny ;
1128  }
1129  if ( ny < 0 && nx > 0 ) {
1130    nx_ = nx ;
1131    ny_ = nx ;
1132  }
1133
1134  //Double wx = xmax - xmin ;
1135  //Double wy = ymax - ymin ;
1136  Double wx = max( abs(xmax-center_(0)), abs(xmin-center_(0)) ) * 2 ;
1137  Double wy = max( abs(ymax-center_(1)), abs(ymin-center_(1)) ) * 2 ;
1138
1139  Quantum<Double> qcellx ;
1140  Quantum<Double> qcelly ;
1141  //cout << "nx_ = " << nx_ << ",  ny_ = " << ny_ << endl ;
1142  if ( cellx.size() != 0 && celly.size() != 0 ) {
1143    readQuantity( qcellx, cellx ) ;
1144    readQuantity( qcelly, celly ) ;
1145  }
1146  else if ( celly.size() != 0 ) {
1147    os << "Using celly to x-axis..." << LogIO::POST ;
1148    readQuantity( qcelly, celly ) ;
1149    qcellx = qcelly ;
1150  }
1151  else if ( cellx.size() != 0 ) {
1152    os << "Using cellx to y-axis..." << LogIO::POST ;
1153    readQuantity( qcellx, cellx ) ;
1154    qcelly = qcellx ;
1155  }
1156  else {
1157    if ( nx_ < 0 ) {
1158      os << "No user preference in grid setting. Using default..." << LogIO::POST ;
1159      readQuantity( qcellx, "1.0arcmin" ) ;
1160      qcelly = qcellx ;
1161    }
1162    else {
1163      if ( wx == 0.0 ) {
1164        os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
1165        wx = 0.00290888 ;
1166      }
1167      if ( wy == 0.0 ) {
1168        os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
1169        wy = 0.00290888 ;
1170      }
1171      if (nx_ > 1) {
1172        qcellx = Quantum<Double>(wx / (nx_ - 1) * cos(center_[1]), "rad");
1173      }
1174      else {
1175        qcellx = Quantum<Double>( 1.1f * wx / nx_ * cos(center_[1]), "rad" );
1176      }
1177      if (ny_ > 1) {
1178        qcelly = Quantum<Double>(wy / (ny_ - 1), "rad");
1179      }
1180      else {
1181        qcelly = Quantum<Double>( 1.1f * wy / ny_, "rad" );
1182      }
1183    }
1184  }
1185  cellx_ = qcellx.getValue( "rad" ) ;
1186  celly_ = qcelly.getValue( "rad" ) ;
1187  //os << "cellx_=" << cellx_ << ", celly_=" << celly_ << ", cos("<<center_(1)<<")=" << cos(center_(1)) << LogIO::POST ;
1188  if ( nx_ < 0 ) {
1189    if ( wx == 0.0 ) {
1190      os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
1191      wx = 0.00290888 ;
1192    }
1193    if ( wy == 0.0 ) {
1194      os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
1195      wy = 0.00290888 ;
1196    }
1197    nx_ = Int( ceil( wx/(cellx_/cos(center_[1])) ) ) ;
1198    ny_ = Int( ceil( wy/celly_ ) ) ;
1199  }
1200
1201  // create DirectionCoordinate
1202  Matrix<Double> xform(2,2) ;
1203  xform = 0.0 ;
1204  xform.diagonal() = 1.0 ;
1205  dircoord_ = new DirectionCoordinate(MDirection::J2000,
1206                                      Projection( Projection::SIN ),
1207                                      center_[0], center_[1],
1208                                      -cellx_, celly_,
1209                                      xform,
1210                                      0.5*Double(nx_-1),
1211                                      0.5*Double(ny_-1)) ;
1212}
1213
1214void STGrid::mapExtent( Double &xmin, Double &xmax,
1215                        Double &ymin, Double &ymax )
1216{
1217  //LogIO os( LogOrigin("STGrid","mapExtent",WHERE) ) ;
1218  directionCol_.attach( tableList_[0], "DIRECTION" ) ;
1219  Matrix<Double> direction = directionCol_.getColumn() ;
1220  //os << "dirCol.nrow() = " << dirCol.nrow() << LogIO::POST ;
1221  Vector<Double> ra( direction.row(0) ) ;
1222  mathutil::rotateRA( ra ) ;
1223  minMax( xmin, xmax, direction.row( 0 ) ) ;
1224  minMax( ymin, ymax, direction.row( 1 ) ) ;
1225  Double amin, amax, bmin, bmax ;
1226  for ( uInt i = 1 ; i < nfile_ ; i++ ) {
1227    directionCol_.attach( tableList_[i], "DIRECTION" ) ;
1228    direction.assign( directionCol_.getColumn() ) ;
1229    //os << "dirCol.nrow() = " << dirCol.nrow() << LogIO::POST ;
1230    // to make contiguous RA distribution (no 2pi jump)
1231    Vector<Double> ra( direction.row(0) ) ;
1232    mathutil::rotateRA( ra ) ;
1233    minMax( amin, amax, direction.row( 0 ) ) ;
1234    minMax( bmin, bmax, direction.row( 1 ) ) ;
1235    xmin = min( xmin, amin ) ;
1236    xmax = max( xmax, amax ) ;
1237    ymin = min( ymin, bmin ) ;
1238    ymax = max( ymax, bmax ) ;
1239  }
1240  //os << "(xmin,xmax)=(" << xmin << "," << xmax << ")" << LogIO::POST ;
1241  //os << "(ymin,ymax)=(" << ymin << "," << ymax << ")" << LogIO::POST ;
1242}
1243
1244void STGrid::table( Table &tab, uInt i )
1245{
1246  if ( i < nfile_ )
1247    tab = Table( infileList_[i] ) ;
1248}
1249
1250void STGrid::selectData()
1251{
1252  LogIO os( LogOrigin("STGrid","selectData",WHERE) ) ;   
1253  Int ifno = ifno_ ;
1254  tableList_.resize( nfile_ ) ;
1255  if ( ifno_ == -1 ) {
1256    //Table taborg( infileList_[0] ) ;
1257    Table taborg ;
1258    table( taborg, 0 ) ;
1259    ROScalarColumn<uInt> ifnoCol( taborg, "IFNO" ) ;
1260    ifno_ = ifnoCol( 0 ) ;
1261    os << LogIO::WARN
1262       << "IFNO is not given. Using default IFNO: " << ifno_ << LogIO::POST ;
1263  }
1264  for ( uInt i = 0 ; i < nfile_ ; i++ ) {
1265    //Table taborg( infileList_[i] ) ;
1266    Table taborg ;
1267    table( taborg, i ) ;
1268    TableExprNode node ;
1269    if ( ifno != -1 || isMultiIF( taborg ) ) {
1270      os << "apply selection on IFNO" << LogIO::POST ;
1271      node = taborg.col("IFNO") == ifno_ ;
1272    }
1273    if ( scanlist_.size() > 0 ) {
1274      os << "apply selection on SCANNO" << LogIO::POST ;
1275      node = node && taborg.col("SCANNO").in( scanlist_ ) ;
1276    }
1277    if ( node.isNull() ) {
1278      tableList_[i] = taborg ;
1279    }
1280    else {
1281      tableList_[i] = taborg( node ) ;
1282    }
1283    os << LogIO::DEBUGGING << "tableList_[" << i << "].nrow()=" << tableList_[i].nrow() << LogIO::POST ;
1284    if ( tableList_[i].nrow() == 0 ) {
1285      os << LogIO::SEVERE
1286         << "No corresponding rows for given selection: IFNO " << ifno_ ;
1287      if ( scanlist_.size() > 0 )
1288        os << " SCANNO " << scanlist_ ;
1289      os << LogIO::EXCEPTION ;
1290    }
1291  }
1292}
1293
1294Bool STGrid::isMultiIF( Table &tab )
1295{
1296  ROScalarColumn<uInt> ifnoCol( tab, "IFNO" ) ;
1297  Vector<uInt> ifnos = ifnoCol.getColumn() ;
1298  return anyNE( ifnos, ifnos[0] ) ;
1299}
1300
1301void STGrid::attach( Table &tab )
1302{
1303  // attach to table
1304  spectraCol_.attach( tab, "SPECTRA" ) ;
1305  flagtraCol_.attach( tab, "FLAGTRA" ) ;
1306  directionCol_.attach( tab, "DIRECTION" ) ;
1307  flagRowCol_.attach( tab, "FLAGROW" ) ;
1308  tsysCol_.attach( tab, "TSYS" ) ;
1309  intervalCol_.attach( tab, "INTERVAL" ) ;
1310}
1311
1312Int STGrid::getDataChunk(
1313                         IPosition const &wshape,
1314                         IPosition const &vshape,
1315                         IPosition const &/* dshape */,
1316                         Array<Complex> &spectra,
1317                         Array<Double> &direction,
1318                         Array<Int> &flagtra,
1319                         Array<Int> &rflag,
1320                         Array<Float> &weight )
1321{
1322  LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
1323
1324  Array<Float> spectraF_(wshape);
1325  Array<uChar> flagtraUC_(wshape);
1326  Array<uInt> rflagUI_(vshape);
1327  Int nrow = getDataChunk( spectraF_, direction, flagtraUC_, rflagUI_, weight ) ;
1328  if ( nrow < nchunk_ ) {
1329    spectra.resize( spectraF_.shape() ) ;
1330    flagtra.resize( flagtraUC_.shape() ) ;
1331    rflag.resize( rflagUI_.shape() ) ;
1332  }
1333  double t0, t1 ;
1334  t0 = mathutil::gettimeofday_sec() ;
1335  convertArray( spectra, spectraF_ ) ;
1336  toInt( flagtraUC_, flagtra ) ;
1337  toInt( rflagUI_, rflag ) ;
1338  t1 = mathutil::gettimeofday_sec() ;
1339  eToInt = t1 - t0 ;
1340 
1341  return nrow ;
1342}
1343
1344#if 0
1345Int STGrid::getDataChunk( Array<Complex> &spectra,
1346                          Array<Double> &direction,
1347                          Array<Int> &flagtra,
1348                          Array<Int> &rflag,
1349                          Array<Float> &weight )
1350{
1351  LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
1352  Int nrow = getDataChunk( spectraF_, direction, flagtraUC_, rflagUI_, weight ) ;
1353  if ( nrow < nchunk_ ) {
1354    spectra.resize( spectraF_.shape() ) ;
1355    flagtra.resize( flagtraUC_.shape() ) ;
1356    rflag.resize( rflagUI_.shape() ) ;
1357  }
1358  double t0, t1 ;
1359  t0 = mathutil::gettimeofday_sec() ;
1360  convertArray( spectra, spectraF_ ) ;
1361  toInt( flagtraUC_, flagtra ) ;
1362  toInt( rflagUI_, rflag ) ;
1363  t1 = mathutil::gettimeofday_sec() ;
1364  eToInt = t1 - t0 ;
1365 
1366  return nrow ;
1367}
1368#endif
1369
1370Int STGrid::getDataChunk( Array<Float> &spectra,
1371                          Array<Double> &direction,
1372                          Array<uChar> &flagtra,
1373                          Array<uInt> &rflag,
1374                          Array<Float> &weight )
1375{
1376  LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
1377  Int nrow = spectra.shape()[1] ;
1378  Int remainingRow = nrow_ - nprocessed_ ;
1379  if ( remainingRow < nrow ) {
1380    nrow = remainingRow ;
1381    IPosition mshape( 2, nchan_, nrow ) ;
1382    IPosition vshape( 1, nrow ) ;
1383    spectra.resize( mshape ) ;
1384    flagtra.resize( mshape ) ;
1385    direction.resize( IPosition(2,2,nrow) ) ;
1386    rflag.resize( vshape ) ;
1387    weight.resize( mshape ) ;
1388  }
1389  // 2011/12/22 TN
1390  // tsys shares its storage with weight
1391  Array<Float> tsys( weight ) ;
1392  Array<Double> tint( rflag.shape() ) ;
1393
1394  Vector<uInt> rflagVec( rflag ) ;
1395  Vector<Double> tintVec( tint ) ;
1396
1397  RefRows rows( nprocessed_, nprocessed_+nrow-1, 1 ) ;
1398  //os<<LogIO::DEBUGGING<<"nprocessed_="<<nprocessed_<<": rows.nrows()="<<rows.nrows()<<LogIO::POST ;
1399  spectraCol_.getColumnCells( rows, spectra ) ;
1400  flagtraCol_.getColumnCells( rows, flagtra ) ;
1401  directionCol_.getColumnCells( rows, direction ) ;
1402  // to make contiguous RA distribution (no 2pi jump)
1403  Vector<Double> v( Matrix<Double>(direction).row(0) ) ;
1404  mathutil::rotateRA( v ) ;
1405  flagRowCol_.getColumnCells( rows, rflagVec ) ;
1406  intervalCol_.getColumnCells( rows, tintVec ) ;
1407  Vector<Float> tsysTemp = tsysCol_( nprocessed_ ) ;
1408  if ( tsysTemp.nelements() == (uInt)nchan_ )
1409    tsysCol_.getColumnCells( rows, tsys ) ;
1410  else
1411    tsys = tsysTemp[0] ;
1412
1413  double t0,t1 ;
1414  t0 = mathutil::gettimeofday_sec() ;
1415  getWeight( weight, tsys, tint ) ;
1416  t1 = mathutil::gettimeofday_sec() ;
1417  eGetWeight += t1-t0 ;
1418
1419  nprocessed_ += nrow ;
1420 
1421  return nrow ;
1422}
1423
1424void STGrid::setupArray()
1425{
1426  LogIO os( LogOrigin("STGrid","setupArray",WHERE) ) ;
1427  ROScalarColumn<uInt> polnoCol( tableList_[0], "POLNO" ) ;
1428  Vector<uInt> pols = polnoCol.getColumn() ;
1429  //os << pols << LogIO::POST ;
1430  Vector<uInt> pollistOrg ;
1431  npolOrg_ = 0 ;
1432  uInt polno ;
1433  for ( uInt i = 0 ; i < polnoCol.nrow() ; i++ ) {
1434    //polno = polnoCol( i ) ;
1435    polno = pols( i ) ;
1436    if ( allNE( pollistOrg, polno ) ) {
1437      pollistOrg.resize( npolOrg_+1, True ) ;
1438      pollistOrg[npolOrg_] = polno ;
1439      npolOrg_++ ;
1440    }
1441  }
1442  if ( pollist_.size() == 0 )
1443    pollist_ = pollistOrg ;
1444  else {
1445    Vector<uInt> newlist ;
1446    uInt newsize = 0 ;
1447    for ( uInt i = 0 ; i < pollist_.size() ; i++ ) {
1448      if ( anyEQ( pollistOrg, pollist_[i] ) ) {
1449        newlist.resize( newsize+1, True ) ;
1450        newlist[newsize] = pollist_[i] ;
1451        newsize++ ;
1452      }
1453    }
1454    pollist_.assign( newlist ) ;
1455  }
1456  npol_ = pollist_.size() ;
1457  if ( npol_ == 0 ) {
1458    os << LogIO::SEVERE << "Empty pollist" << LogIO::EXCEPTION ;
1459  }
1460  rows_.resize( nfile_ ) ;
1461  for ( uInt i = 0 ; i < nfile_ ; i++ ) {
1462    rows_[i] = tableList_[i].nrow() / npolOrg_ ;
1463    //if ( nrow_ < rows_[i] )
1464    //  nrow_ = rows_[i] ;
1465  }
1466  flagtraCol_.attach( tableList_[0], "FLAGTRA" ) ;
1467  nchan_ = flagtraCol_( 0 ).nelements() ;
1468//   os << "npol_ = " << npol_ << "(" << pollist_ << ")" << endl
1469//      << "nchan_ = " << nchan_ << endl
1470//      << "nrow_ = " << nrow_ << LogIO::POST ;
1471}
1472
1473void STGrid::getWeight( Array<Float> &w,
1474                              Array<Float> &tsys,
1475                              Array<Double> &tint )
1476{
1477  LogIO os( LogOrigin("STGrid","getWeight",WHERE) ) ;
1478
1479  // 2011/12/22 TN
1480  // w (weight) and tsys share storage
1481  IPosition refShape = tsys.shape() ;
1482  Int nchan = refShape[0] ;
1483  Int nrow = refShape[1] ;
1484//   os << "nchan=" << nchan << ", nrow=" << nrow << LogIO::POST ;
1485//   os << "w.shape()=" << w.shape() << endl
1486//      << "tsys.shape()=" << tsys.shape() << endl
1487//      << "tint.shape()=" << tint.shape() << LogIO::POST ;
1488
1489  // set weight
1490  if ( wtype_.compare( "UNIFORM" ) == 0 ) {
1491    w = 1.0 ;
1492  }
1493  else if ( wtype_.compare( "TINT" ) == 0 ) {
1494    Bool b0, b1 ;
1495    Float *w_p = w.getStorage( b0 ) ;
1496    Float *w0_p = w_p ;
1497    const Double *ti_p = tint.getStorage( b1 ) ;
1498    const Double *w1_p = ti_p ;
1499    for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1500      for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1501        *w0_p = *w1_p ;
1502        w0_p++ ;
1503      }
1504      w1_p++ ;
1505    }
1506    w.putStorage( w_p, b0 ) ;
1507    tint.freeStorage( ti_p, b1 ) ;
1508  }
1509  else if ( wtype_.compare( "TSYS" ) == 0 ) {
1510    Bool b0 ;
1511    Float *w_p = w.getStorage( b0 ) ;
1512    Float *w0_p = w_p ;
1513    for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1514      for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1515        Float temp = *w0_p ;
1516        *w0_p = 1.0 / ( temp * temp ) ;
1517        w0_p++ ;
1518      }
1519    }
1520    w.putStorage( w_p, b0 ) ;
1521  }
1522  else if ( wtype_.compare( "TINTSYS" ) == 0 ) {
1523    Bool b0, b1 ;
1524    Float *w_p = w.getStorage( b0 ) ;
1525    Float *w0_p = w_p ;
1526    const Double *ti_p = tint.getStorage( b1 ) ;
1527    const Double *w1_p = ti_p ;
1528    for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1529      Float interval = *w1_p ;
1530      for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1531        Float temp = *w0_p ;
1532        *w0_p = interval / ( temp * temp ) ;
1533        w0_p++ ;
1534      }
1535      w1_p++ ;
1536    }
1537    w.putStorage( w_p, b0 ) ;
1538    tint.freeStorage( ti_p, b1 ) ;
1539  }
1540  else {
1541    //LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
1542    //os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
1543    w = 1.0 ;
1544  }
1545}
1546
1547void STGrid::toInt( Array<uChar> &u, Array<Int> &v )
1548{
1549  uInt len = u.nelements() ;
1550  Int *int_p = new Int[len] ;
1551  Bool deleteIt ;
1552  const uChar *data_p = u.getStorage( deleteIt ) ;
1553  Int *i_p = int_p ;
1554  const uChar *u_p = data_p ;
1555  for ( uInt i = 0 ; i < len ; i++ ) {
1556    *i_p = ( *u_p == 0 ) ? 0 : 1 ;
1557    i_p++ ;
1558    u_p++ ;
1559  }
1560  u.freeStorage( data_p, deleteIt ) ;
1561  v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
1562}
1563
1564void STGrid::toInt( Array<uInt> &u, Array<Int> &v )
1565{
1566  uInt len = u.nelements() ;
1567  Int *int_p = new Int[len] ;
1568  Bool deleteIt ;
1569  const uInt *data_p = u.getStorage( deleteIt ) ;
1570  Int *i_p = int_p ;
1571  const uInt *u_p = data_p ;
1572  for ( uInt i = 0 ; i < len ; i++ ) {
1573    *i_p = ( *u_p == 0 ) ? 0 : 1 ;
1574    i_p++ ;
1575    u_p++ ;
1576  }
1577  u.freeStorage( data_p, deleteIt ) ;
1578  v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
1579}
1580
1581void STGrid::toPixel( Array<Double> &world, Array<Double> &pixel )
1582{
1583  uInt nrow = world.shape()[1] ;
1584  Bool bw, bp ;
1585  Double *w_p = world.getStorage( bw ) ;
1586  Double *p_p = pixel.getStorage( bp ) ;
1587  Double *ww_p = w_p ;
1588  Double *wp_p = p_p ;
1589  IPosition vshape( 1, 2 ) ;
1590  Vector<Double> _world, _pixel ;
1591  for ( uInt i = 0 ; i < nrow ; i++ ) {
1592    _world.takeStorage( vshape, ww_p, SHARE ) ;
1593    _pixel.takeStorage( vshape, wp_p, SHARE ) ;
1594    dircoord_->toPixel( _pixel, _world ) ;
1595    ww_p += 2 ;
1596    wp_p += 2 ;
1597  }
1598  world.putStorage( w_p, bw ) ;
1599  pixel.putStorage( p_p, bp ) ;
1600}
1601
1602void STGrid::boxFunc( Vector<Float> &convFunc, Int &convSize )
1603{
1604  convFunc = 0.0 ;
1605  for ( Int i = 0 ; i < convSize/2 ; i++ )
1606    convFunc(i) = 1.0 ;
1607}
1608
1609#define NEED_UNDERSCORES
1610#if defined(NEED_UNDERSCORES)
1611#define grdsf grdsf_
1612#define grdgauss grdgauss_
1613#define grdjinc1 grdjinc1_
1614#endif
1615#if defined(USE_CASAPY)
1616extern "C" {
1617  void grdsf(Double*, Double*);
1618  void grdgauss(Double*, Double*, Double*);
1619  void grdjinc1(Double*, Double*, Int*, Double*);
1620}
1621#else
1622extern "C" {
1623  void grdsf(Double*, Double*);
1624}
1625void grdgauss(Double *hwhm, Double *val, Double *out)
1626{
1627  *out = exp(-log(2.0) * (*val / *hwhm) * (*val / *hwhm));
1628}
1629void grdjinc1(Double *c, Double *val, Int *normalize, Double *out)
1630{
1631  // Calculate J_1(x) using approximate formula
1632  Double x = C::pi * *val / *c;
1633  Double ax = fabs(x);
1634  Double ans;
1635  if ( ax < 8.0 ) {
1636    Double y = x * x;
1637    Double ans1 = x * (72362614232.0 + y * (-7895059235.0
1638                       + y * (242396853.1 + y * (-2972611.439
1639                       + y * (15704.48260 + y * (-30.16036606))))));
1640    Double ans2 = 144725228442.0 + y * (2300535178.0
1641                       + y * (18583304.74 + y * (99447.43394
1642                       + y * (376.9991397 + y * 1.0))));
1643    ans = ans1 / ans2;
1644  }
1645  else {
1646    Double z = 8.0 / ax;
1647    Double y = z * z;
1648    Double xx = ax - 2.356194491;
1649    Double ans1 = 1.0 + y * (0.183105e-2 + y * (-0.3516396496e-4
1650                      + y * (0.2457520174e-5 + y * (-0.240337019e-6))));
1651    Double ans2 = 0.04687499995 + y * (-0.2002690873e-3
1652                      + y * (0.8449199096e-5 + y * (-0.88228987e-6
1653                      + y * (0.105787412e-6))));
1654    ans = sqrt(0.636619772 / ax) * (cos(xx) * ans1
1655                                    - z * sin(xx) * ans2);
1656    if (x < 0.0)
1657      ans = -ans;
1658  }
1659 
1660  // Then, calculate Jinc
1661  if (x == 0.0) {
1662    *out = 0.5;
1663  }
1664  else {
1665    *out = ans / x;
1666  }
1667
1668  if (*normalize == 1)
1669    *out = *out / 0.5; 
1670}
1671#endif
1672void STGrid::spheroidalFunc( Vector<Float> &convFunc )
1673{
1674  convFunc = 0.0 ;
1675  for ( Int i = 0 ; i < convSampling_*convSupport_ ; i++ ) {
1676    Double nu = Double(i) / Double(convSupport_*convSampling_) ;
1677    Double val ;
1678    grdsf( &nu, &val ) ;
1679    convFunc(i) = ( 1.0 - nu * nu ) * val ;
1680  }
1681}
1682
1683void STGrid::gaussFunc( Vector<Float> &convFunc, Double hwhm, Double truncate )
1684{
1685  convFunc = 0.0 ;
1686  Int len = (Int)(truncate*Double(convSampling_)+0.5);
1687  Double out, val;
1688  for ( Int i = 0 ; i < len ; i++ ) {
1689    val = Double(i) / Double(convSampling_) ;
1690    grdgauss(&hwhm, &val, &out);
1691    convFunc(i) = out;
1692  }
1693}
1694
1695void STGrid::gjincFunc( Vector<Float> &convFunc, Double hwhm, Double c, Double truncate )
1696{
1697  convFunc = 0.0;
1698  Double out1, out2, val;
1699  Int normalize = 1;
1700  if  (truncate >= 0.0) {
1701    Int len = (Int)(truncate*Double(convSampling_)+0.5);
1702    for (Int i = 0 ; i < len ; i++) {
1703      val = Double(i) / Double(convSampling_);
1704      grdgauss(&hwhm, &val, &out1);
1705      grdjinc1(&c, &val, &normalize, &out2);
1706      convFunc(i) = out1 * out2;
1707    }
1708  }
1709  else {
1710    Int len = convFunc.nelements();
1711    for (Int i = 0 ; i < len ; i++) {
1712      val = Double(i) / Double(convSampling_);
1713      grdjinc1(&c, &val, &normalize, &out2);
1714      if (out2 <= 0.0) {
1715        LogIO os(LogOrigin("STGrid","gjincFunc",WHERE));
1716        os << LogIO::DEBUG1 << "convFunc is automatically truncated at radius " << val << LogIO::POST;
1717        break;
1718      }
1719      grdgauss(&hwhm, &val, &out1);
1720      convFunc(i) = out1 * out2;
1721    }
1722  }
1723}
1724
1725void STGrid::pbFunc( Vector<Float> &convFunc )
1726{
1727  convFunc = 0.0 ;
1728}
1729
1730vector<float> STGrid::getConvFunc()
1731{
1732  LogIO os(LogOrigin("STGrid","getConvFunc",WHERE));
1733  Vector<Float> convFunc;
1734  vector<float> out;
1735
1736  if (cellx_ <= 0.0 || celly_ <= 0.0) {
1737    selectData();
1738    setupGrid();
1739  }
1740
1741  if (convType_ == "BOX" || convType_ == "SF") {
1742    setConvFunc(convFunc);
1743  }
1744  else if (convType_ == "GAUSS") {
1745    Quantum<Double> q1,q2;
1746    readQuantity(q1,gwidth_);
1747    readQuantity(q2,truncate_);
1748//     if (celly_ <= 0.0
1749//         && ((!q1.getUnit().empty()&&q1.getUnit()!="pixel") ||
1750//             (!q2.getUnit().empty()&&q2.getUnit()!="pixel"))) {
1751//       throw AipsError("You have to call defineImage to get correct convFunc");
1752//     }
1753    setConvFunc(convFunc);
1754  }
1755  else if (convType_ == "GJINC") {
1756    Quantum<Double> q1,q2,q3;
1757    readQuantity(q1,gwidth_);
1758    readQuantity(q2,truncate_);
1759    readQuantity(q3,jwidth_);
1760//     if (celly_ <= 0.0
1761//         && ((!q1.getUnit().empty()&&q1.getUnit()!="pixel") ||
1762//             (!q2.getUnit().empty()&&q2.getUnit()!="pixel") ||
1763//             (!q3.getUnit().empty()&&q3.getUnit()!="pixel"))) {
1764//       throw AipsError("You have to call defineImage to get correct convFunc");
1765//     }
1766    setConvFunc(convFunc);
1767  }
1768  else if (convType_ == "PB") {
1769    throw AipsError("Grid function PB is not available");
1770  }
1771  else {
1772    throw AipsError("Unknown grid function: "+convType_);
1773  }
1774
1775  convFunc.tovector(out);
1776  return out;
1777}
1778
1779void STGrid::setConvFunc( Vector<Float> &convFunc )
1780{
1781  LogIO os(LogOrigin("STGrid","setConvFunc",WHERE));
1782  convSupport_ = userSupport_ ;
1783  if ( convType_ == "BOX" ) {
1784    if ( convSupport_ < 0 )
1785      convSupport_ = 0 ;
1786    Int convSize = convSampling_ * ( 2 * convSupport_ + 2 )  ;
1787    convFunc.resize( convSize ) ;
1788    boxFunc( convFunc, convSize ) ;
1789    os << LogIO::DEBUGGING
1790       << "convType_ = " << convType_ << endl
1791       << "convSupport_ = " << convSupport_ << LogIO::POST;
1792  }
1793  else if ( convType_ == "SF" ) {
1794    if ( convSupport_ < 0 )
1795      convSupport_ = 3 ;
1796    Int convSize = convSampling_ * ( 2 * convSupport_ + 2 )  ;
1797    convFunc.resize( convSize ) ;
1798    spheroidalFunc( convFunc ) ;
1799    os << LogIO::DEBUGGING
1800       << "convType_ = " << convType_ << endl
1801       << "convSupport_ = " << convSupport_ << LogIO::POST;
1802  }
1803  else if ( convType_ == "GAUSS" ) {
1804    // determine pixel gwidth
1805    // default is HWHM corresponding to b = 1.0 (Mangum et al. 2007)
1806    Double pixelGW = -1.0;
1807    Quantum<Double> q ;
1808    if (!gwidth_.empty()) {
1809      readQuantity( q, gwidth_ );
1810      if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
1811        pixelGW = q.getValue();
1812      }
1813      else {
1814        pixelGW = q.getValue("rad")/celly_;
1815      }
1816    }
1817    pixelGW = (pixelGW >= 0.0) ? pixelGW : sqrt(log(2.0));
1818    if (pixelGW < 0.0) {
1819      os << LogIO::SEVERE
1820         << "Negative width is specified for gaussian" << LogIO::EXCEPTION;
1821    }
1822    // determine truncation radius
1823    // default is 3 * HWHM
1824    Double truncate = -1.0;
1825    if (!truncate_.empty()) {
1826      readQuantity( q, truncate_ );
1827      if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
1828        truncate = q.getValue();
1829      }
1830      else {
1831        truncate = q.getValue("rad")/celly_;
1832      }
1833    }     
1834    //convSupport_ = (Int)(truncate+0.5);
1835    truncate = (truncate >= 0.0) ? truncate : 3.0 * pixelGW;
1836    convSupport_ = Int(truncate);
1837    convSupport_ += (((truncate-(Double)convSupport_) > 0.0) ? 1 : 0);
1838    Int convSize = convSampling_ * ( 2*convSupport_ + 2 ) ;
1839    convFunc.resize( convSize ) ;
1840    gaussFunc( convFunc, pixelGW, truncate ) ;
1841    os << LogIO::DEBUGGING
1842       << "convType_ = " << convType_ << endl
1843       << "convSupport_ = " << convSupport_ << endl
1844       << "truncate_ = " << truncate << "pixel" << endl
1845       << "gwidth_ = " << pixelGW << "pixel" << LogIO::POST;
1846  }
1847  else if ( convType_ == "GJINC" ) {
1848    // determine pixel gwidth
1849    // default is HWHM corresponding to b = 2.52 (Mangum et al. 2007)
1850    Double pixelGW = -1.0;
1851    Quantum<Double> q ;
1852    if (!gwidth_.empty()) {
1853      readQuantity( q, gwidth_ );
1854      if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
1855        pixelGW = q.getValue();
1856      }
1857      else {
1858        pixelGW = q.getValue("rad")/celly_;
1859      }
1860    }
1861    pixelGW = (pixelGW >= 0.0) ? pixelGW : sqrt(log(2.0)) * 2.52;
1862    if (pixelGW < 0.0) {
1863      os << LogIO::SEVERE
1864         << "Negative width is specified for gaussian" << LogIO::EXCEPTION;
1865    }
1866    // determine pixel c
1867    // default is c = 1.55 (Mangum et al. 2007)
1868    Double pixelJW = -1.0;
1869    if (!jwidth_.empty()) {
1870      readQuantity( q, jwidth_ );
1871      if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
1872        pixelJW = q.getValue();
1873      }
1874      else {
1875        pixelJW = q.getValue("rad")/celly_;
1876      }
1877    }
1878    pixelJW = (pixelJW >= 0.0) ? pixelJW : 1.55;
1879    if (pixelJW < 0.0) {
1880      os << LogIO::SEVERE
1881         << "Negative width is specified for jinc" << LogIO::EXCEPTION;
1882    }
1883    // determine truncation radius
1884    // default is -1.0 (truncate at first null)
1885    Double truncate = -1.0;
1886    if (!truncate_.empty()) {
1887      readQuantity( q, truncate_ );
1888      if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
1889        truncate = q.getValue();
1890      }
1891      else {
1892        truncate = q.getValue("rad")/celly_;
1893      }
1894    }     
1895    //convSupport_ = (truncate >= 0.0) ? (Int)(truncate+0.5) : (Int)(2*pixelJW+0.5);
1896    Double convSupportF = (truncate >= 0.0) ? truncate : (2*pixelJW);
1897    convSupport_ = (Int)convSupportF;
1898    convSupport_ += (((convSupportF-(Double)convSupport_) > 0.0) ? 1 : 0);
1899    Int convSize = convSampling_ * ( 2*convSupport_ + 2 ) ;
1900    convFunc.resize( convSize ) ;
1901    gjincFunc( convFunc, pixelGW, pixelJW, truncate ) ;
1902    os << LogIO::DEBUGGING
1903       << "convType_ = " << convType_ << endl
1904       << "convSupport_ = " << convSupport_ << endl
1905       << "truncate_ = " << truncate << "pixel" << endl
1906       << "gwidth_ = " << pixelGW << "pixel" << endl
1907       << "jwidth_ = " << pixelJW << "pixel" << LogIO::POST;
1908  }
1909  else if ( convType_ == "PB" ) {
1910    if ( convSupport_ < 0 )
1911      convSupport_ = 0 ;
1912    pbFunc( convFunc ) ;
1913  }
1914  else {
1915    throw AipsError( "Unsupported convolution function" ) ;
1916  }
1917}
1918
1919string STGrid::saveData( string outfile )
1920{
1921  LogIO os( LogOrigin("STGrid", "saveData", WHERE) ) ;
1922  double t0, t1 ;
1923  t0 = mathutil::gettimeofday_sec() ;
1924
1925  //Int polno = 0 ;
1926  String outfile_ ;
1927  if ( outfile.size() == 0 ) {
1928    if ( infileList_[0].lastchar() == '/' ) {
1929      outfile_ = infileList_[0].substr( 0, infileList_[0].size()-1 ) ;
1930    }
1931    else {
1932      outfile_ = infileList_[0] ;
1933    }
1934    outfile_ += ".grid" ;
1935  }
1936  else {
1937    outfile_ = outfile ;
1938  }
1939  Table tab ;
1940  prepareTable( tab, outfile_ ) ;
1941  fillTable( tab ) ;
1942
1943  t1 = mathutil::gettimeofday_sec() ;
1944  os << LogIO::DEBUGGING << "saveData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
1945
1946  return outfile_ ;
1947}
1948
1949void STGrid::prepareTable( Table &tab, String &name )
1950{
1951  Table t( infileList_[0], Table::Old ) ;
1952  t.deepCopy( name, Table::New, False, t.endianFormat(), True ) ;
1953  tab = Table( name, Table::Update ) ;
1954  // 2012/02/13 TN
1955  // explicitly copy subtables since no rows including subtables are
1956  // copied by Table::deepCopy with noRows=True
1957  //TableCopy::copySubTables( tab, t ) ;
1958  const TableRecord &inrec = t.keywordSet();
1959  TableRecord &outrec = tab.rwKeywordSet();
1960  for (uInt i = 0 ; i < inrec.nfields() ; i++) {
1961    if (inrec.type(i) == TpTable) {
1962      String name = inrec.name(i);
1963      Table intable = inrec.asTable(name);
1964      Table outtable = outrec.asTable(name);
1965      TableCopy::copyRows(outtable, intable);
1966    }
1967  }
1968}
1969
1970void STGrid::fillTable( Table &tab )
1971{
1972  //IPosition dshape = data_.shape() ;
1973  Int nrow = nx_ * ny_ * npol_ ;
1974  tab.rwKeywordSet().define( "nPol", npol_ ) ;
1975  tab.addRow( nrow ) ;
1976  Vector<Double> cpix( 2 ) ;
1977  cpix(0) = Double( nx_ - 1 ) * 0.5 ;
1978  cpix(1) = Double( ny_ - 1 ) * 0.5 ;
1979  Vector<Double> dir( 2 ) ;
1980  Vector<Double> pix( 2 );
1981  ArrayColumn<Double> directionCol( tab, "DIRECTION" ) ;
1982  ArrayColumn<Float> spectraCol( tab, "SPECTRA" ) ;
1983  ArrayColumn<uChar> flagtraCol( tab, "FLAGTRA" ) ;
1984  ScalarColumn<uInt> flagRowCol( tab, "FLAGROW" );
1985  ScalarColumn<uInt> polnoCol( tab, "POLNO" ) ;
1986  ScalarColumn<uInt> scannoCol( tab, "SCANNO" ) ;
1987  Int irow = 0 ;
1988  Vector<Float> sp( nchan_ ) ;
1989  Vector<uChar> flag( nchan_ ) ;
1990  Bool bsp, bdata, bflag ;
1991  const Float *data_p = data_.getStorage( bdata ) ;
1992  const uChar *flag_p = flag_.getStorage( bflag ) ;
1993  Float *wsp_p, *sp_p ;
1994  const Float *wdata_p = data_p ;
1995  const uChar *wflag_p = flag_p ;
1996  long step = nx_ * ny_ * npol_ ;
1997  long offset ;
1998  uInt scanno = 0 ;
1999  uChar rflag;
2000  for ( Int iy = 0 ; iy < ny_ ; iy++ ) {
2001    pix(1) = (Double)(iy);
2002    for ( Int ix = 0 ; ix < nx_ ; ix++ ) {
2003      pix(0) = (Double)(ix);
2004      dircoord_->toWorld(dir,pix);
2005      //os << "dir[" << ix << "," << iy << "]=" << dir << LogIO::POST;
2006      for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
2007        offset = ix + nx_ * (iy + ipol * ny_) ;
2008        //os << "offset = " << offset << LogIO::POST ;
2009        sp_p = sp.getStorage( bsp ) ;
2010        wsp_p = sp_p ;
2011        wdata_p = data_p + offset ;
2012        wflag_p = flag_p + offset ;
2013        rflag = ~0 ; //11111111
2014        for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
2015          *wsp_p = *wdata_p ;
2016          wsp_p++ ;
2017          wdata_p += step ;
2018          flag[ichan] = *wflag_p ;
2019          rflag &= flag[ichan] ;
2020          wflag_p += step ;
2021        }
2022        sp.putStorage( sp_p, bsp ) ;
2023        spectraCol.put( irow, sp ) ;
2024        flagtraCol.put( irow, flag ) ;
2025        flagRowCol.put( irow, ((rflag>0) ? 1 : 0) ) ;
2026        directionCol.put( irow, dir ) ;
2027        polnoCol.put( irow, pollist_[ipol] ) ;
2028        scannoCol.put( irow, scanno ) ;
2029        irow++ ;
2030      }
2031      scanno++ ;
2032    }
2033  }
2034  data_.freeStorage( data_p, bdata ) ;
2035  flag_.freeStorage( flag_p, bflag ) ;
2036
2037  fillMainColumns( tab ) ;
2038}
2039
2040void STGrid::fillMainColumns( Table &tab )
2041{
2042  // values for fill
2043  //Table t( infileList_[0], Table::Old ) ;
2044  Table t ;
2045  table( t, 0 ) ;
2046  Table tsel = t( t.col( "IFNO" ) == (uInt)ifno_, 1 ) ;
2047  ROTableRow row( tsel ) ;
2048  row.get( 0 ) ;
2049  const TableRecord &rec = row.record() ;
2050  uInt freqId = rec.asuInt( "FREQ_ID" ) ;
2051  uInt molId = rec.asuInt( "MOLECULE_ID" ) ;
2052  uInt tcalId = rec.asuInt( "TCAL_ID" ) ;
2053  uInt focusId = rec.asuInt( "FOCUS_ID" ) ;
2054  uInt weatherId = rec.asuInt( "WEATHER_ID" ) ;
2055  String srcname = rec.asString( "SRCNAME" ) ;
2056  String fieldname = rec.asString( "FIELDNAME" ) ;
2057  Vector<Float> defaultTsys( 1, 1.0 ) ;
2058  // @todo how to set flagtra for gridded spectra?
2059  Vector<uChar> flagtra = rec.asArrayuChar( "FLAGTRA" ) ;
2060  flagtra = (uChar)0 ;
2061  Float opacity = rec.asFloat( "OPACITY" ) ;
2062  Double srcvel = rec.asDouble( "SRCVELOCITY" ) ;
2063  Vector<Double> srcpm = rec.asArrayDouble( "SRCPROPERMOTION" ) ;
2064  Vector<Double> srcdir = rec.asArrayDouble( "SRCDIRECTION" ) ;
2065  Vector<Double> scanrate = rec.asArrayDouble( "SCANRATE" ) ;
2066  Double time = rec.asDouble( "TIME" ) ;
2067  Double interval = rec.asDouble( "INTERVAL" ) ;
2068
2069  // fill columns
2070  Int nrow = tab.nrow() ;
2071  ScalarColumn<uInt> ifnoCol( tab, "IFNO" ) ;
2072  ScalarColumn<uInt> beamnoCol(tab, "BEAMNO");
2073  ScalarColumn<uInt> freqIdCol( tab, "FREQ_ID" ) ;
2074  ScalarColumn<uInt> molIdCol( tab, "MOLECULE_ID" ) ;
2075  ScalarColumn<uInt> tcalidCol( tab, "TCAL_ID" ) ;
2076  ScalarColumn<Int> fitidCol( tab, "FIT_ID" ) ;
2077  ScalarColumn<uInt> focusidCol( tab, "FOCUS_ID" ) ;
2078  ScalarColumn<uInt> weatheridCol( tab, "WEATHER_ID" ) ;
2079  ArrayColumn<uChar> flagtraCol( tab, "FLAGTRA" ) ;
2080  ScalarColumn<uInt> rflagCol( tab, "FLAGROW" ) ;
2081  ArrayColumn<Float> tsysCol( tab, "TSYS" ) ;
2082  ScalarColumn<String> srcnameCol( tab, "SRCNAME" ) ;
2083  ScalarColumn<String> fieldnameCol( tab, "FIELDNAME" ) ;
2084  ScalarColumn<Int> srctypeCol( tab, "SRCTYPE" ) ;
2085  ScalarColumn<Float> opacityCol( tab, "OPACITY" ) ;
2086  ScalarColumn<Double> srcvelCol( tab, "SRCVELOCITY" ) ;
2087  ArrayColumn<Double> srcpmCol( tab, "SRCPROPERMOTION" ) ;
2088  ArrayColumn<Double> srcdirCol( tab, "SRCDIRECTION" ) ;
2089  ArrayColumn<Double> scanrateCol( tab, "SCANRATE" ) ;
2090  ScalarColumn<Double> timeCol( tab, "TIME" ) ;
2091  ScalarColumn<Double> intervalCol( tab, "INTERVAL" ) ;
2092  for ( Int i = 0 ; i < nrow ; i++ ) {
2093    ifnoCol.put( i, (uInt)ifno_ ) ;
2094    beamnoCol.put(i, 0);
2095    freqIdCol.put( i, freqId ) ;
2096    molIdCol.put( i, molId ) ;
2097    tcalidCol.put( i, tcalId ) ;
2098    fitidCol.put( i, -1 ) ;
2099    focusidCol.put( i, focusId ) ;
2100    weatheridCol.put( i, weatherId ) ;
2101    //flagtraCol.put( i, flagtra ) ;
2102    //rflagCol.put( i, 0 ) ;
2103    tsysCol.put( i, defaultTsys ) ;
2104    srcnameCol.put( i, srcname ) ;
2105    fieldnameCol.put( i, fieldname ) ;
2106    srctypeCol.put( i, (Int)SrcType::PSON ) ;
2107    opacityCol.put( i, opacity ) ;
2108    srcvelCol.put( i, srcvel ) ;
2109    srcpmCol.put( i, srcpm ) ;
2110    srcdirCol.put( i, srcdir ) ;
2111    scanrateCol.put( i, scanrate ) ;
2112    timeCol.put( i, time ) ;
2113    intervalCol.put( i, interval ) ;
2114    if ( (i + 1) % npol_ == 0 ) {
2115      time += interval / 86400.0;
2116    }
2117  }
2118}
2119
2120vector<int> STGrid::getResultantMapSize()
2121{
2122  vector<int> r(2);
2123  r[0] = nx_;
2124  r[1] = ny_;
2125  return r;
2126}
2127
2128vector<double> STGrid::getResultantCellSize()
2129{
2130  vector<double> r(2);
2131  r[0] = cellx_;
2132  r[1] = celly_;
2133  return r;
2134}
2135
2136// STGrid2
2137STGrid2::STGrid2()
2138  : STGrid()
2139{
2140}
2141
2142STGrid2::STGrid2( const ScantableWrapper &s )
2143  : STGrid()
2144{
2145  setScantable( s ) ;
2146}
2147
2148STGrid2::STGrid2( const vector<ScantableWrapper> &v )
2149  : STGrid()
2150{
2151  setScantableList( v ) ;
2152}
2153
2154void STGrid2::setScantable( const ScantableWrapper &s )
2155{
2156  nfile_ = 1 ;
2157  dataList_.resize( nfile_ ) ;
2158  dataList_[0] = s ;
2159  infileList_.resize( nfile_ ) ;
2160  infileList_[0] = s.getCP()->table().tableName() ;
2161}
2162
2163void STGrid2::setScantableList( const vector<ScantableWrapper> &v )
2164{
2165  nfile_ = v.size() ;
2166  dataList_.resize( nfile_ ) ;
2167  infileList_.resize( nfile_ ) ;
2168  for ( uInt i = 0 ; i < nfile_ ; i++ ) {
2169    dataList_[i] = v[i] ;
2170    infileList_[i] = v[i].getCP()->table().tableName() ;
2171  }
2172}
2173
2174ScantableWrapper STGrid2::getResultAsScantable( int tp )
2175{
2176  ScantableWrapper sw( tp ) ;
2177  CountedPtr<Scantable> s = sw.getCP() ;
2178  s->setHeader( dataList_[0].getCP()->getHeader() ) ;
2179  Table tout, tin ;
2180  String subt[] = { "FREQUENCIES", "FOCUS", "WEATHER",
2181                    "TCAL", "MOLECULES", "HISTORY", "FIT" } ;
2182  for ( uInt i = 0 ; i < 7 ; i++ ) {
2183    tout = s->table().rwKeywordSet().asTable(subt[i]) ;
2184    tin = dataList_[0].getCP()->table().rwKeywordSet().asTable(subt[i]) ;
2185    TableCopy::copyRows( tout, tin ) ;
2186    tout.rwKeywordSet() = tin.rwKeywordSet();
2187  }
2188  fillTable( s->table() ) ;
2189  return sw ;
2190}
2191
2192void STGrid2::table( Table &tab, uInt i )
2193{
2194  if ( i < nfile_ )
2195    tab = dataList_[i].getCP()->table() ;
2196}
2197
2198}
Note: See TracBrowser for help on using the repository browser.