source: trunk/src/STGrid.cpp @ 2393

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

New Development: No

JIRA Issue: Yes CAS-2816

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

Parallel version of STGrid. Implemented using concurrent module.
Two threads are running in parallel. Those are responsible for
retrieving data and processing data respectively.


File size: 39.1 KB
Line 
1#include <iostream>
2#include <fstream>
3
4#include <casa/BasicSL/String.h>
5#include <casa/Arrays/Vector.h>
6#include <casa/Arrays/Matrix.h>
7#include <casa/Arrays/Cube.h>
8#include <casa/Arrays/ArrayMath.h>
9#include <casa/Arrays/ArrayIter.h>
10#include <casa/Quanta/Quantum.h>
11#include <casa/Quanta/QuantumHolder.h>
12#include <casa/Utilities/CountedPtr.h>
13#include <casa/Logging/LogIO.h>
14
15#include <tables/Tables/Table.h>
16#include <tables/Tables/TableRecord.h>
17#include <tables/Tables/ExprNode.h>
18#include <tables/Tables/ScalarColumn.h>
19#include <tables/Tables/ArrayColumn.h>
20#include <tables/Tables/TableIter.h>
21
22#include <measures/Measures/MDirection.h>
23
24#include <MathUtils.h>
25
26#include "STGrid.h"
27
28using namespace std ;
29using namespace concurrent ;
30using namespace casa ;
31using namespace asap ;
32
33namespace asap {
34
35// for performance check
36double eToInt = 0.0 ;
37double eGetWeight = 0.0 ;
38
39// constructor
40STGrid::STGrid()
41  : vshape_( 1 ), wshape_( 2 ), dshape_( 2 )
42{
43  init() ;
44}
45
46STGrid::STGrid( const string infile )
47  : vshape_( 1 ), wshape_( 2 ), dshape_( 2 )
48{
49  init() ;
50
51  setFileIn( infile ) ;
52}
53
54STGrid::STGrid( const vector<string> infile )
55{
56  init() ;
57
58  setFileList( infile ) ;
59}
60
61void  STGrid::init()
62{
63  ifno_ = -1 ;
64  nx_ = -1 ;
65  ny_ = -1 ;
66  npol_ = 0 ;
67  nchan_ = 0 ;
68  nrow_ = 0 ;
69  cellx_ = 0.0 ;
70  celly_ = 0.0 ;
71  center_ = Vector<Double> ( 2, 0.0 ) ;
72  convType_ = "BOX" ;
73  wtype_ = "UNIFORM" ;
74  convSupport_ = -1 ;
75  userSupport_ = -1 ;
76  convSampling_ = 100 ;
77  nprocessed_ = 0 ;
78  nchunk_ = 0 ;
79
80  // initialize user input
81  nxUI_ = -1 ;
82  nyUI_ = -1 ;
83  cellxUI_ = "" ;
84  cellyUI_ = "" ;
85  centerUI_ = "" ;
86}
87
88void STGrid::setFileIn( const string infile )
89{
90  nfile_ = 1 ;
91  String name( infile ) ;
92  infileList_.resize( nfile_ ) ;
93  infileList_[0] = String(infile) ;
94}
95
96void STGrid::setFileList( const vector<string> infile )
97{
98  nfile_ = infile.size() ;
99  infileList_.resize( nfile_ ) ;
100  for ( uInt i = 0 ; i < nfile_ ; i++ ) {
101    infileList_[i] = infile[i] ;
102  }
103}
104
105void STGrid::setPolList( vector<unsigned int> pols )
106{
107  pollist_.assign( Vector<uInt>( pols ) ) ;
108}
109
110void STGrid::setScanList( vector<unsigned int> scans )
111{
112  scanlist_.assign( Vector<uInt>( scans ) ) ;
113}
114
115void STGrid::setWeight( const string wType )
116{
117  wtype_ = String( wType ) ;
118  wtype_.upcase() ;
119}
120
121void STGrid::defineImage( int nx,
122                          int ny,
123                          string scellx,
124                          string scelly,
125                          string scenter )
126{
127  nxUI_ = (Int)nx ;
128  nyUI_ = (Int)ny ;
129  cellxUI_ = String( scellx ) ;
130  cellyUI_ = String( scelly ) ;
131  centerUI_ = String( scenter ) ;
132}
133 
134void STGrid::setFunc( string convType,
135                      int convSupport )
136{
137  convType_ = String( convType ) ;
138  convType_.upcase() ;
139  userSupport_ = (Int)convSupport ;
140}
141
142#define NEED_UNDERSCORES
143#if defined(NEED_UNDERSCORES)
144#define ggridsd ggridsd_
145#endif
146extern "C" {
147   void ggridsd(Double*,
148                const Complex*,
149                Int*,
150                Int*,
151                Int*,
152                const Int*,
153                const Int*,
154                const Float*,
155                Int*,
156                Int*,
157                Complex*,
158                Float*,
159                Int*,
160                Int*,
161                Int *,
162                Int *,
163                Int*,
164                Int*,
165                Float*,
166                Int*,
167                Int*,
168                Double*);
169}
170void STGrid::call_ggridsd( Array<Double> &xypos,
171                           Array<Complex> &spectra,
172                           Int &nvispol,
173                           Int &nvischan,
174                           Array<Int> &flagtra,
175                           Array<Int> &flagrow,
176                           Array<Float> &weight,
177                           Int &nrow,
178                           Int &irow,
179                           Array<Complex> &gdata,
180                           Array<Float> &gwgt,
181                           Int &nx,
182                           Int &ny,
183                           Int &npol,
184                           Int &nchan,
185                           Int &support,
186                           Int &sampling,
187                           Vector<Float> &convFunc,
188                           Int *chanMap,
189                           Int *polMap )
190{
191  // parameters for gridding
192  Int idopsf = 0 ;
193  Int len = npol*nchan ;
194  Double *sumw_p = new Double[len] ;
195  {
196    Double *work_p = sumw_p ;
197    for ( Int i = 0 ; i < len ; i++ ) {
198      *work_p = 0.0 ;
199      work_p++ ;
200    }
201  }
202
203  // prepare pointer
204  Bool deletePos, deleteData, deleteWgt, deleteFlag, deleteFlagR, deleteConv, deleteDataG, deleteWgtG ;
205  Double *xy_p = xypos.getStorage( deletePos ) ;
206  const Complex *values_p = spectra.getStorage( deleteData ) ;
207  const Int *flag_p = flagtra.getStorage( deleteFlag ) ;
208  const Int *rflag_p = flagrow.getStorage( deleteFlagR ) ;
209  const Float *wgt_p = weight.getStorage( deleteWgt ) ;
210  Complex *grid_p = gdata.getStorage( deleteDataG ) ;
211  Float *wgrid_p = gwgt.getStorage( deleteWgtG ) ;
212  Float *conv_p = convFunc.getStorage( deleteConv ) ;
213
214  // pass copy of irow to ggridsd since it will be modified in theroutine
215  Int irowCopy = irow ;
216     
217  // call ggridsd
218  ggridsd( xy_p,
219           values_p,
220           &nvispol,
221           &nvischan,
222           &idopsf,
223           flag_p,
224           rflag_p,
225           wgt_p,
226           &nrow,
227           &irowCopy,
228           grid_p,
229           wgrid_p,
230           &nx,
231           &ny,
232           &npol,
233           &nchan,
234           &support,
235           &sampling,
236           conv_p,
237           chanMap,
238           polMap,
239           sumw_p ) ;
240
241  // finalization
242  xypos.putStorage( xy_p, deletePos ) ;
243  spectra.freeStorage( values_p, deleteData ) ;
244  flagtra.freeStorage( flag_p, deleteFlag ) ;
245  flagrow.freeStorage( rflag_p, deleteFlagR ) ;
246  weight.freeStorage( wgt_p, deleteWgt ) ;
247  gdata.putStorage( grid_p, deleteDataG ) ;
248  gwgt.putStorage( wgrid_p, deleteWgtG ) ;
249  convFunc.putStorage( conv_p, deleteConv ) ;
250  delete sumw_p ;
251}
252
253
254void STGrid::grid()
255{
256  LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
257  double t0,t1 ;
258
259  // data selection
260  t0 = mathutil::gettimeofday_sec() ;
261  selectData() ;
262  t1 = mathutil::gettimeofday_sec() ;
263  os << "selectData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
264
265  setupGrid() ;
266  setupArray() ;
267
268  if ( wtype_.compare("UNIFORM") != 0 &&
269       wtype_.compare("TINT") != 0 &&
270       wtype_.compare("TSYS") != 0 &&
271       wtype_.compare("TINTSYS") != 0 ) {
272    LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
273    os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
274    wtype_ = "UNIFORM" ;
275  }
276
277  // grid parameter
278  os << LogIO::DEBUGGING ;
279  os << "----------" << endl ;
280  os << "Data selection summary" << endl ;
281  os << "   ifno = " << ifno_ << endl ;
282  os << "   pollist = " << pollist_ << endl ;
283  os << "   scanlist = " << scanlist_ << endl ;
284  os << "----------" << endl ;
285  os << "Grid parameter summary" << endl ;
286  os << "   (nx,ny) = (" << nx_ << "," << ny_ << ")" << endl ;
287  os << "   (cellx,celly) = (" << cellx_ << "," << celly_ << ")" << endl ;
288  os << "   center = " << center_ << endl ;
289  os << "   weighting = " << wtype_ << endl ;
290  os << "   convfunc = " << convType_ << " with support " << convSupport_ << endl ;
291  os << "----------" << LogIO::POST ;
292  os << LogIO::NORMAL ;
293
294  Bool doAll = examine() ;
295
296  if ( doAll )
297    gridPerPol() ;
298  else
299    gridPerRow() ;
300}
301
302Bool STGrid::examine()
303{
304  // TODO: nchunk_ must be determined from nchan_, npol_, and (nx_,ny_)
305  //       by considering data size to be allocated for ggridsd input/output
306  nchunk_ = 400 ;
307  Bool b = nchunk_ >= nrow_ ;
308  nchunk_ = min( nchunk_, nrow_ ) ;
309  vshape_ = IPosition( 1, nchunk_ ) ;
310  wshape_ = IPosition( 2, nchan_, nchunk_ ) ;
311  dshape_ = IPosition( 2, 2, nchunk_ ) ;
312  return b ;
313}
314
315struct STGChunk {
316  Int nrow ;
317  Array<Complex> spectra;
318  Array<Int> flagtra;
319  Array<Int> rflag;
320  Array<Float> weight;
321  Array<Double> direction;
322  STGChunk(IPosition const &wshape, IPosition const &vshape,
323           IPosition const &dshape)
324    : spectra(wshape), flagtra(wshape), rflag(vshape), weight(wshape),
325      direction(dshape)
326  { }
327};
328
329struct STCommonData {
330  Int gnx;
331  Int gny;
332  Int *chanMap;
333  Vector<Float> convFunc ;
334  Array<Complex> gdataArrC;
335  Array<Float> gwgtArr;
336  STCommonData(IPosition const &gshape, Array<Float> const &data)
337    : gdataArrC(gshape, 0.0), gwgtArr(data) {}
338};
339
340#define DO_AHEAD 3
341
342struct STContext {
343  STCommonData &common;
344  FIFO<STGChunk *, DO_AHEAD> queue;
345  STGrid *const self;
346  const Int pol;
347  STContext(STGrid *obj, STCommonData &common, Int pol)
348    : self(obj), common(common), pol(pol) {}
349};
350
351
352bool STGrid::produceChunk(void *ctx) throw(PCException)
353{
354  STContext &context = *(STContext *)ctx;
355  if ( context.self->nprocessed_ >= context.self->nrow_ ) {
356    return false;
357  }
358  STGChunk *chunk = new STGChunk(context.self->wshape_,
359                                 context.self->vshape_,
360                                 context.self->dshape_);
361
362  double t0 = mathutil::gettimeofday_sec() ;
363  chunk->nrow = context.self->getDataChunk(
364        context.self->wshape_, context.self->vshape_, context.self->dshape_,
365        chunk->spectra, chunk->direction,
366        chunk->flagtra, chunk->rflag, chunk->weight);
367  double t1 = mathutil::gettimeofday_sec() ;
368  context.self->eGetData_ += t1-t0 ;
369
370  context.queue.lock();
371  context.queue.put(chunk);
372  context.queue.unlock();
373  return true;
374}
375
376void STGrid::consumeChunk(void *ctx) throw(PCException)
377{
378  STContext &context = *(STContext *)ctx;
379  STGChunk *chunk = NULL;
380  try {
381    context.queue.lock();
382    chunk = context.queue.get();
383    context.queue.unlock();
384  } catch (FullException &e) {
385    context.queue.unlock();
386    // TODO: log error
387    throw PCException();
388  }
389
390  double t0, t1 ;
391  // world -> pixel
392  Array<Double> xypos( context.self->dshape_ ) ;
393  t0 = mathutil::gettimeofday_sec() ;
394  context.self->toPixel( chunk->direction, xypos ) ;
395  t1 = mathutil::gettimeofday_sec() ;
396  context.self->eToPixel_ += t1-t0 ;
397   
398  // call ggridsd
399  Int nvispol = 1 ;
400  Int irow = -1 ;
401  t0 = mathutil::gettimeofday_sec() ;
402  context.self->call_ggridsd( xypos,
403                chunk->spectra,
404                nvispol,
405                context.self->nchan_,
406                chunk->flagtra,
407                chunk->rflag,
408                chunk->weight,
409                chunk->nrow,
410                irow,
411                context.common.gdataArrC,
412                context.common.gwgtArr,
413                context.common.gnx,
414                context.common.gny,
415                context.self->npol_,
416                context.self->nchan_,
417                context.self->convSupport_,
418                context.self->convSampling_,
419                context.common.convFunc,
420                context.common.chanMap,
421                (Int*)&context.pol ) ;
422  t1 = mathutil::gettimeofday_sec() ;
423  context.self->eGGridSD_ += t1-t0 ;
424 
425  delete chunk;
426}
427
428void STGrid::gridPerRow()
429{
430  LogIO os( LogOrigin("STGrid", "gridPerRow", WHERE) ) ;
431  double t0, t1 ;
432
433
434  // grid data
435  // Extend grid plane with convSupport_
436  //   Int gnx = nx_+convSupport_*2 ;
437  //   Int gny = ny_+convSupport_*2 ;
438  Int gnx = nx_;
439  Int gny = ny_;
440
441  IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
442  // 2011/12/20 TN
443  // data_ and gwgtArr share storage
444  data_.resize( gshape ) ;
445  data_ = 0.0 ;
446  STCommonData common = STCommonData(gshape, data_);
447  common.gnx = gnx ;
448  common.gny = gny ;
449
450  // parameters for gridding
451  Int *chanMap = new Int[nchan_] ;
452  for ( Int i = 0 ; i < nchan_ ; i++ ) {
453    chanMap[i] = i ;
454  }
455  common.chanMap = chanMap;
456
457  // convolution kernel
458  t0 = mathutil::gettimeofday_sec() ;
459  setConvFunc( common.convFunc ) ;
460  t1 = mathutil::gettimeofday_sec() ;
461  os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
462
463  // for performance check
464  eGetData_ = 0.0 ;
465  eToPixel_ = 0.0 ;
466  eGGridSD_ = 0.0 ;
467  double eInitPol = 0.0 ;
468
469  Broker broker = Broker(produceChunk, consumeChunk);
470  for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) {
471    initTable( ifile ) ;
472
473    os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ;   
474    for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
475      t0 = mathutil::gettimeofday_sec() ;
476      initPol( ipol ) ; // set ptab_ and attach()
477      t1 = mathutil::gettimeofday_sec() ;
478      eInitPol += t1-t0 ;
479     
480      STContext context(this, common, ipol);
481     
482      os << "start pol " << ipol << LogIO::POST ;
483     
484      nprocessed_ = 0 ;
485#if 1
486      broker.runProducerAsMasterThread(&context, DO_AHEAD);
487#else
488      for (;;) {
489        bool produced = produceChunk(&context);
490        if (! produced) {
491          break;
492        }
493        consumeChunk(&context);
494      }
495#endif
496
497      os << "end pol " << ipol << LogIO::POST ;
498
499    }
500    os << "end table " << ifile << LogIO::POST ;   
501  }
502  os << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ;
503  os << "getData: elapsed time is " << eGetData_-eToInt-eGetWeight << " sec." << LogIO::POST ;
504  os << "toPixel: elapsed time is " << eToPixel_ << " sec." << LogIO::POST ;
505  os << "ggridsd: elapsed time is " << eGGridSD_ << " sec." << LogIO::POST ;
506  os << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
507  os << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
508 
509  delete chanMap ;
510
511  // set data
512  setData( common.gdataArrC, common.gwgtArr ) ;
513
514}
515
516void STGrid::gridPerPol()
517{
518  LogIO os( LogOrigin("STGrid", "gridPerPol", WHERE) ) ;
519  double t0, t1 ;
520
521  // convolution kernel
522  Vector<Float> convFunc ;
523  t0 = mathutil::gettimeofday_sec() ;
524  setConvFunc( convFunc ) ;
525  t1 = mathutil::gettimeofday_sec() ;
526  os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
527
528  // prepare grid data storage
529  Int gnx = nx_ ;
530  Int gny = ny_ ;
531//   // Extend grid plane with convSupport_
532//   Int gnx = nx_+convSupport_*2 ;
533//   Int gny = ny_+convSupport_*2 ;
534  IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
535  Array<Complex> gdataArrC( gshape, 0.0 ) ;
536  //Array<Float> gwgtArr( gshape, 0.0 ) ;
537  // 2011/12/20 TN
538  // data_ and weight array shares storage
539  data_.resize( gshape ) ;
540  data_ = 0.0 ;
541  Array<Float> gwgtArr( data_ ) ;
542
543  // maps
544  Int *chanMap = new Int[nchan_] ;
545  {
546    Int *work_p = chanMap ;
547    for ( Int i = 0 ; i < nchan_ ; i++ ) {
548      *work_p = i ;
549      work_p++ ;
550    }
551  }
552  Int polMap[1] ;
553
554  // some parameters for ggridsd
555  Int nvispol = 1 ;
556  Int irow = -1 ;
557 
558  // for performance check
559  double eInitPol = 0.0 ;
560  double eGetData = 0.0 ;
561  double eToPixel = 0.0 ;
562  double eGGridSD = 0.0 ;
563  double eToInt = 0.0 ;
564
565  for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) {
566    initTable( ifile ) ;
567
568    os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ;   
569    // to read data from the table
570    IPosition mshape( 2, nchan_, nrow_ ) ;
571    IPosition dshape( 2, 2, nrow_ ) ;
572    Array<Complex> spectra( mshape ) ;
573    Array<Double> direction( dshape ) ;
574    Array<Int> flagtra( mshape ) ;
575    Array<Int> rflag( IPosition(1,nrow_) ) ;
576    Array<Float> weight( mshape ) ;
577    Array<Double> xypos( dshape ) ;
578
579    for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
580      t0 = mathutil::gettimeofday_sec() ;
581      initPol( ipol ) ;
582      t1 = mathutil::gettimeofday_sec() ;
583      eInitPol += t1-t0 ;
584     
585      os << "start pol " << ipol << LogIO::POST ;
586     
587      // retrieve data
588      t0 = mathutil::gettimeofday_sec() ;
589      getDataPerPol( spectra, direction, flagtra, rflag, weight ) ;
590      t1 = mathutil::gettimeofday_sec() ;
591      eGetData += t1-t0 ;
592     
593      // world -> pixel
594      t0 = mathutil::gettimeofday_sec() ;
595      toPixel( direction, xypos ) ; 
596      t1 = mathutil::gettimeofday_sec() ;
597      eToPixel += t1-t0 ;
598     
599      // call ggridsd
600      polMap[0] = ipol ;
601      t0 = mathutil::gettimeofday_sec() ;
602      call_ggridsd( xypos,
603                    spectra,
604                    nvispol,
605                    nchan_,
606                    flagtra,
607                    rflag,
608                    weight,
609                    nrow_,
610                    irow,
611                    gdataArrC,
612                    gwgtArr,
613                    gnx,
614                    gny,
615                    npol_,
616                    nchan_,
617                    convSupport_,
618                    convSampling_,
619                    convFunc,
620                    chanMap,
621                    polMap ) ;
622      t1 = mathutil::gettimeofday_sec() ;
623      eGGridSD += t1-t0 ;
624
625      os << "end pol " << ipol << LogIO::POST ;
626
627    }
628    os << "end table " << ifile << LogIO::POST ;   
629  }
630  os << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ;
631  os << "getData: elapsed time is " << eGetData-eToInt-eGetWeight << " sec." << LogIO::POST ;
632  os << "toPixel: elapsed time is " << eToPixel << " sec." << LogIO::POST ;
633  os << "ggridsd: elapsed time is " << eGGridSD << " sec." << LogIO::POST ;
634  os << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
635  os << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
636
637  // delete maps
638  delete chanMap ;
639
640  setData( gdataArrC, gwgtArr ) ;
641//   os << "gdataArr = " << gdataArr << LogIO::POST ;
642//   os << "gwgtArr = " << gwgtArr << LogIO::POST ;
643//   os << "data_ " << data_ << LogIO::POST ;
644}
645
646void STGrid::initPol( Int ipol )
647{
648  LogIO os( LogOrigin("STGrid","initPol",WHERE) ) ;
649  if ( npolOrg_ == 1 ) {
650    os << "single polarization data." << LogIO::POST ;
651    ptab_ = tab_ ;
652  }
653  else
654    ptab_ = tab_( tab_.col("POLNO") == (uInt)ipol ) ;
655
656  attach( ptab_ ) ;
657}
658
659void STGrid::initTable( uInt idx )
660{
661  tab_ = tableList_[idx] ;
662  nrow_ = rows_[idx] ;
663}
664
665void STGrid::setData( Array<Complex> &gdata,
666                      Array<Float> &gwgt )
667{
668  // 2011/12/20 TN
669  // gwgt and data_ share storage
670  LogIO os( LogOrigin("STGrid","setData",WHERE) ) ;
671  double t0, t1 ;
672  t0 = mathutil::gettimeofday_sec() ;
673  uInt len = data_.nelements() ;
674  const Complex *w1_p ;
675  Float *w2_p ;
676  Bool b1, b2 ;
677  const Complex *gdata_p = gdata.getStorage( b1 ) ;
678  Float *gwgt_p = data_.getStorage( b2 ) ;
679  w1_p = gdata_p ;
680  w2_p = gwgt_p ;
681  for ( uInt i = 0 ; i < len ; i++ ) {
682    if ( *w2_p > 0.0 ) *w2_p = (*w1_p).real() / *w2_p ;
683    w1_p++ ;
684    w2_p++ ;
685  }
686  gdata.freeStorage( gdata_p, b1 ) ;
687  data_.putStorage( gwgt_p, b2 ) ;
688  t1 = mathutil::gettimeofday_sec() ;
689  os << "setData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
690}
691
692void STGrid::setupGrid()
693{
694  Double xmin,xmax,ymin,ymax ;
695  mapExtent( xmin, xmax, ymin, ymax ) ;
696 
697  setupGrid( nxUI_, nyUI_, cellxUI_, cellyUI_,
698             xmin, xmax, ymin, ymax, centerUI_ ) ;
699}
700
701void STGrid::setupGrid( Int &nx,
702                        Int &ny,
703                        String &cellx,
704                        String &celly,
705                        Double &xmin,
706                        Double &xmax,
707                        Double &ymin,
708                        Double &ymax,
709                        String &center )
710{
711  LogIO os( LogOrigin("STGrid","setupGrid",WHERE) ) ;
712  //cout << "nx=" << nx << ", ny=" << ny << endl ;
713
714  // center position
715  if ( center.size() == 0 ) {
716    center_(0) = 0.5 * ( xmin + xmax ) ;
717    center_(1) = 0.5 * ( ymin + ymax ) ;
718  }
719  else {
720    String::size_type pos0 = center.find( " " ) ;
721    if ( pos0 == String::npos ) {
722      throw AipsError( "bad string format in parameter center" ) ;
723    }
724    String::size_type pos1 = center.find( " ", pos0+1 ) ;
725    String typestr, xstr, ystr ;
726    if ( pos1 != String::npos ) {
727      typestr = center.substr( 0, pos0 ) ;
728      xstr = center.substr( pos0+1, pos1-pos0 ) ;
729      ystr = center.substr( pos1+1 ) ;
730      // todo: convert to J2000 (or direction ref for DIRECTION column)
731    }
732    else {
733      typestr = "J2000" ;
734      xstr = center.substr( 0, pos0 ) ;
735      ystr = center.substr( pos0+1 ) ;
736    }
737    QuantumHolder qh ;
738    String err ;
739    qh.fromString( err, xstr ) ;
740    Quantum<Double> xcen = qh.asQuantumDouble() ;
741    qh.fromString( err, ystr ) ;
742    Quantum<Double> ycen = qh.asQuantumDouble() ;
743    center_(0) = xcen.getValue( "rad" ) ;
744    center_(1) = ycen.getValue( "rad" ) ;
745  }
746
747
748  nx_ = nx ;
749  ny_ = ny ;
750  if ( nx < 0 && ny > 0 ) {
751    nx_ = ny ;
752    ny_ = ny ;
753  }
754  if ( ny < 0 && nx > 0 ) {
755    nx_ = nx ;
756    ny_ = nx ;
757  }
758
759  //Double wx = xmax - xmin ;
760  //Double wy = ymax - ymin ;
761  Double wx = max( abs(xmax-center_(0)), abs(xmin-center_(0)) ) * 2 ;
762  Double wy = max( abs(ymax-center_(1)), abs(ymin-center_(1)) ) * 2 ;
763  // take 10% margin
764  wx *= 1.10 ;
765  wy *= 1.10 ;
766
767  Quantum<Double> qcellx ;
768  Quantum<Double> qcelly ;
769  //cout << "nx_ = " << nx_ << ",  ny_ = " << ny_ << endl ;
770  if ( cellx.size() != 0 && celly.size() != 0 ) {
771    readQuantity( qcellx, cellx ) ;
772    readQuantity( qcelly, celly ) ;
773  }
774  else if ( celly.size() != 0 ) {
775    os << "Using celly to x-axis..." << LogIO::POST ;
776    readQuantity( qcelly, celly ) ;
777    qcellx = qcelly ;
778  }
779  else if ( cellx.size() != 0 ) {
780    os << "Using cellx to y-axis..." << LogIO::POST ;
781    readQuantity( qcellx, cellx ) ;
782    qcelly = qcellx ;
783  }
784  else {
785    if ( nx_ < 0 ) {
786      os << "No user preference in grid setting. Using default..." << LogIO::POST ;
787      readQuantity( qcellx, "1.0arcmin" ) ;
788      qcelly = qcellx ;
789    }
790    else {
791      if ( wx == 0.0 ) {
792        os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
793        wx = 0.00290888 ;
794      }
795      if ( wy == 0.0 ) {
796        os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
797        wy = 0.00290888 ;
798      }
799      qcellx = Quantum<Double>( wx/nx_, "rad" ) ;
800      qcelly = Quantum<Double>( wy/ny_, "rad" ) ;
801    }
802  }
803  cellx_ = qcellx.getValue( "rad" ) ;
804  celly_ = qcelly.getValue( "rad" ) ;
805  if ( nx_ < 0 ) {
806    if ( wx == 0.0 ) {
807      os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
808      wx = 0.00290888 ;
809    }
810    if ( wy == 0.0 ) {
811      os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
812      wy = 0.00290888 ;
813    }
814    nx_ = Int( ceil( wx/cellx_ ) ) ;
815    ny_ = Int( ceil( wy/celly_ ) ) ;
816  }
817}
818
819void STGrid::mapExtent( Double &xmin, Double &xmax,
820                        Double &ymin, Double &ymax )
821{
822  //LogIO os( LogOrigin("STGrid","mapExtent",WHERE) ) ;
823  directionCol_.attach( tableList_[0], "DIRECTION" ) ;
824  Matrix<Double> direction = directionCol_.getColumn() ;
825  //os << "dirCol.nrow() = " << dirCol.nrow() << LogIO::POST ;
826  minMax( xmin, xmax, direction.row( 0 ) ) ;
827  minMax( ymin, ymax, direction.row( 1 ) ) ;
828  Double amin, amax, bmin, bmax ;
829  for ( uInt i = 1 ; i < nfile_ ; i++ ) {
830    directionCol_.attach( tableList_[i], "DIRECTION" ) ;
831    direction.assign( directionCol_.getColumn() ) ;
832    //os << "dirCol.nrow() = " << dirCol.nrow() << LogIO::POST ;
833    minMax( amin, amax, direction.row( 0 ) ) ;
834    minMax( bmin, bmax, direction.row( 1 ) ) ;
835    xmin = min( xmin, amin ) ;
836    xmax = max( xmax, amax ) ;
837    ymin = min( ymin, bmin ) ;
838    ymax = max( ymax, bmax ) ;
839  }
840  //os << "(xmin,xmax)=(" << xmin << "," << xmax << ")" << LogIO::POST ;
841  //os << "(ymin,ymax)=(" << ymin << "," << ymax << ")" << LogIO::POST ;
842}
843
844void STGrid::selectData()
845{
846  LogIO os( LogOrigin("STGrid","selectData",WHERE) ) ;   
847  Int ifno = ifno_ ;
848  tableList_.resize( nfile_ ) ;
849  if ( ifno == -1 ) {
850    Table taborg( infileList_[0] ) ;
851    ROScalarColumn<uInt> ifnoCol( taborg, "IFNO" ) ;
852    ifno = ifnoCol( 0 ) ;
853    os << LogIO::WARN
854       << "IFNO is not given. Using default IFNO: " << ifno << LogIO::POST ;
855  }
856  for ( uInt i = 0 ; i < nfile_ ; i++ ) {
857    Table taborg( infileList_[i] ) ;
858    TableExprNode node ;
859    if ( isMultiIF( taborg ) ) {
860      os << "apply selection on IFNO" << LogIO::POST ;
861      node = taborg.col("IFNO") == ifno ;
862    }
863    if ( scanlist_.size() > 0 ) {
864      os << "apply selection on SCANNO" << LogIO::POST ;
865      node = node && taborg.col("SCANNO").in( scanlist_ ) ;
866    }
867    if ( node.isNull() ) {
868      tableList_[i] = taborg ;
869    }
870    else {
871      tableList_[i] = taborg( node ) ;
872    }
873    os << "tableList_[" << i << "].nrow()=" << tableList_[i].nrow() << LogIO::POST ;
874    if ( tableList_[i].nrow() == 0 ) {
875      os << LogIO::SEVERE
876         << "No corresponding rows for given selection: IFNO " << ifno ;
877      if ( scanlist_.size() > 0 )
878        os << " SCANNO " << scanlist_ ;
879      os << LogIO::EXCEPTION ;
880    }
881  }
882}
883
884Bool STGrid::isMultiIF( Table &tab )
885{
886  ROScalarColumn<uInt> ifnoCol( tab, "IFNO" ) ;
887  Vector<uInt> ifnos = ifnoCol.getColumn() ;
888  return anyNE( ifnos, ifnos[0] ) ;
889}
890
891void STGrid::attach( Table &tab )
892{
893  // attach to table
894  spectraCol_.attach( tab, "SPECTRA" ) ;
895  flagtraCol_.attach( tab, "FLAGTRA" ) ;
896  directionCol_.attach( tab, "DIRECTION" ) ;
897  flagRowCol_.attach( tab, "FLAGROW" ) ;
898  tsysCol_.attach( tab, "TSYS" ) ;
899  intervalCol_.attach( tab, "INTERVAL" ) ;
900}
901
902void STGrid::getDataPerPol( Array<Float> &spectra,
903                            Array<Double> &direction,
904                            Array<uChar> &flagtra,
905                            Array<uInt> &rflag,
906                            Array<Float> &weight )
907{
908  LogIO os( LogOrigin("STGrid","getDataPerPol",WHERE) ) ;
909
910  // 2011/12/22 TN
911  // weight and tsys shares its storage
912  Array<Float> tsys( weight ) ;
913  Array<Double> tint( rflag.shape() ) ;
914
915  Vector<uInt> rflagVec( rflag ) ;
916  Vector<Double> tintVec( tint ) ;
917
918  spectraCol_.getColumn( spectra ) ;
919  flagtraCol_.getColumn( flagtra ) ;
920  flagRowCol_.getColumn( rflagVec ) ;
921  directionCol_.getColumn( direction ) ;
922  intervalCol_.getColumn( tintVec ) ;
923  Vector<Float> tsysTemp = tsysCol_( 0 ) ;
924  if ( tsysTemp.nelements() == (uInt)nchan_ ) {
925    tsysCol_.getColumn( tsys ) ;
926  }
927  else {
928    tsys = tsysTemp[0] ;
929  }
930 
931  double t0,t1 ;
932  t0 = mathutil::gettimeofday_sec() ;
933  getWeight( weight, tsys, tint ) ;
934  t1 = mathutil::gettimeofday_sec() ;
935  eGetWeight += t1-t0 ;
936}
937
938void STGrid::getDataPerPol( Array<Complex> &spectra,
939                            Array<Double> &direction,
940                            Array<Int> &flagtra,
941                            Array<Int> &rflag,
942                            Array<Float> &weight )
943{
944  LogIO os( LogOrigin("STGrid","getDataPerPol",WHERE) ) ;
945  double t0, t1 ;
946
947  Array<uChar> flagUC( flagtra.shape() ) ;
948  Array<uInt> rflagUI( rflag.shape() ) ;
949  Array<Float> spectraF( spectra.shape() ) ;
950  getDataPerPol( spectraF, direction, flagUC, rflagUI, weight ) ;
951
952  convertArray( spectra, spectraF ) ;
953  t0 = mathutil::gettimeofday_sec() ;
954  toInt( flagUC, flagtra ) ;
955  toInt( rflagUI, rflag ) ;
956  t1 = mathutil::gettimeofday_sec() ;
957  eToInt += t1-t0 ;
958  //os << "toInt: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
959}
960
961Int STGrid::getDataChunk(
962                         IPosition const &wshape,
963                         IPosition const &vshape,
964                         IPosition const &dshape,
965                         Array<Complex> &spectra,
966                         Array<Double> &direction,
967                         Array<Int> &flagtra,
968                         Array<Int> &rflag,
969                         Array<Float> &weight )
970{
971  LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
972
973  Array<Float> spectraF_(wshape);
974  Array<uChar> flagtraUC_(wshape);
975  Array<uInt> rflagUI_(vshape);
976  Int nrow = getDataChunk( spectraF_, direction, flagtraUC_, rflagUI_, weight ) ;
977  if ( nrow < nchunk_ ) {
978    spectra.resize( spectraF_.shape() ) ;
979    flagtra.resize( flagtraUC_.shape() ) ;
980    rflag.resize( rflagUI_.shape() ) ;
981  }
982  double t0, t1 ;
983  t0 = mathutil::gettimeofday_sec() ;
984  convertArray( spectra, spectraF_ ) ;
985  toInt( flagtraUC_, flagtra ) ;
986  toInt( rflagUI_, rflag ) ;
987  t1 = mathutil::gettimeofday_sec() ;
988  eToInt = t1 - t0 ;
989 
990  return nrow ;
991}
992
993#if 0
994Int STGrid::getDataChunk( Array<Complex> &spectra,
995                          Array<Double> &direction,
996                          Array<Int> &flagtra,
997                          Array<Int> &rflag,
998                          Array<Float> &weight )
999{
1000  LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
1001  Int nrow = getDataChunk( spectraF_, direction, flagtraUC_, rflagUI_, weight ) ;
1002  if ( nrow < nchunk_ ) {
1003    spectra.resize( spectraF_.shape() ) ;
1004    flagtra.resize( flagtraUC_.shape() ) ;
1005    rflag.resize( rflagUI_.shape() ) ;
1006  }
1007  double t0, t1 ;
1008  t0 = mathutil::gettimeofday_sec() ;
1009  convertArray( spectra, spectraF_ ) ;
1010  toInt( flagtraUC_, flagtra ) ;
1011  toInt( rflagUI_, rflag ) ;
1012  t1 = mathutil::gettimeofday_sec() ;
1013  eToInt = t1 - t0 ;
1014 
1015  return nrow ;
1016}
1017#endif
1018
1019Int STGrid::getDataChunk( Array<Float> &spectra,
1020                          Array<Double> &direction,
1021                          Array<uChar> &flagtra,
1022                          Array<uInt> &rflag,
1023                          Array<Float> &weight )
1024{
1025  LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
1026  Int nrow = spectra.shape()[1] ;
1027  Int remainingRow = nrow_ - nprocessed_ ;
1028  if ( remainingRow < nrow ) {
1029    nrow = remainingRow ;
1030    IPosition mshape( 2, nchan_, nrow ) ;
1031    IPosition vshape( 1, nrow ) ;
1032    spectra.resize( mshape ) ;
1033    flagtra.resize( mshape ) ;
1034    direction.resize( IPosition(2,2,nrow) ) ;
1035    rflag.resize( vshape ) ;
1036    weight.resize( mshape ) ;
1037  }
1038  // 2011/12/22 TN
1039  // tsys shares its storage with weight
1040  Array<Float> tsys( weight ) ;
1041  Array<Double> tint( rflag.shape() ) ;
1042
1043  Vector<uInt> rflagVec( rflag ) ;
1044  Vector<Double> tintVec( tint ) ;
1045
1046  RefRows rows( nprocessed_, nprocessed_+nrow-1, 1 ) ;
1047  //os<<LogIO::DEBUGGING<<"nprocessed_="<<nprocessed_<<": rows.nrows()="<<rows.nrows()<<LogIO::POST ;
1048  spectraCol_.getColumnCells( rows, spectra ) ;
1049  flagtraCol_.getColumnCells( rows, flagtra ) ;
1050  directionCol_.getColumnCells( rows, direction ) ;
1051  flagRowCol_.getColumnCells( rows, rflagVec ) ;
1052  intervalCol_.getColumnCells( rows, tintVec ) ;
1053  Vector<Float> tsysTemp = tsysCol_( nprocessed_ ) ;
1054  if ( tsysTemp.nelements() == (uInt)nchan_ )
1055    tsysCol_.getColumnCells( rows, tsys ) ;
1056  else
1057    tsys = tsysTemp[0] ;
1058
1059  double t0,t1 ;
1060  t0 = mathutil::gettimeofday_sec() ;
1061  getWeight( weight, tsys, tint ) ;
1062  t1 = mathutil::gettimeofday_sec() ;
1063  eGetWeight += t1-t0 ;
1064
1065  nprocessed_ += nrow ;
1066 
1067  return nrow ;
1068}
1069
1070void STGrid::setupArray()
1071{
1072  LogIO os( LogOrigin("STGrid","setupArray",WHERE) ) ;
1073  ROScalarColumn<uInt> polnoCol( tableList_[0], "POLNO" ) ;
1074  Vector<uInt> pols = polnoCol.getColumn() ;
1075  //os << pols << LogIO::POST ;
1076  Vector<uInt> pollistOrg ;
1077  npolOrg_ = 0 ;
1078  uInt polno ;
1079  for ( uInt i = 0 ; i < polnoCol.nrow() ; i++ ) {
1080    //polno = polnoCol( i ) ;
1081    polno = pols( i ) ;
1082    if ( allNE( pollistOrg, polno ) ) {
1083      pollistOrg.resize( npolOrg_+1, True ) ;
1084      pollistOrg[npolOrg_] = polno ;
1085      npolOrg_++ ;
1086    }
1087  }
1088  if ( pollist_.size() == 0 )
1089    pollist_ = pollistOrg ;
1090  else {
1091    Vector<uInt> newlist ;
1092    uInt newsize = 0 ;
1093    for ( uInt i = 0 ; i < pollist_.size() ; i++ ) {
1094      if ( anyEQ( pollistOrg, pollist_[i] ) ) {
1095        newlist.resize( newsize+1, True ) ;
1096        newlist[newsize] = pollist_[i] ;
1097        newsize++ ;
1098      }
1099    }
1100    pollist_.assign( newlist ) ;
1101  }
1102  npol_ = pollist_.size() ;
1103  if ( npol_ == 0 ) {
1104    os << LogIO::SEVERE << "Empty pollist" << LogIO::EXCEPTION ;
1105  }
1106  rows_.resize( nfile_ ) ;
1107  for ( uInt i = 0 ; i < nfile_ ; i++ ) {
1108    rows_[i] = tableList_[i].nrow() / npolOrg_ ;
1109    if ( nrow_ < rows_[i] )
1110      nrow_ = rows_[i] ;
1111  }
1112  flagtraCol_.attach( tableList_[0], "FLAGTRA" ) ;
1113  nchan_ = flagtraCol_( 0 ).nelements() ;
1114//   os << "npol_ = " << npol_ << "(" << pollist_ << ")" << endl
1115//      << "nchan_ = " << nchan_ << endl
1116//      << "nrow_ = " << nrow_ << LogIO::POST ;
1117}
1118
1119void STGrid::getWeight( Array<Float> &w,
1120                              Array<Float> &tsys,
1121                              Array<Double> &tint )
1122{
1123  LogIO os( LogOrigin("STGrid","getWeight",WHERE) ) ;
1124
1125  // 2011/12/22 TN
1126  // w (weight) and tsys share storage
1127  IPosition refShape = tsys.shape() ;
1128  Int nchan = refShape[0] ;
1129  Int nrow = refShape[1] ;
1130//   os << "nchan=" << nchan << ", nrow=" << nrow << LogIO::POST ;
1131//   os << "w.shape()=" << w.shape() << endl
1132//      << "tsys.shape()=" << tsys.shape() << endl
1133//      << "tint.shape()=" << tint.shape() << LogIO::POST ;
1134
1135  // set weight
1136  if ( wtype_.compare( "UNIFORM" ) == 0 ) {
1137    w = 1.0 ;
1138  }
1139  else if ( wtype_.compare( "TINT" ) == 0 ) {
1140    Bool b0, b1 ;
1141    Float *w_p = w.getStorage( b0 ) ;
1142    Float *w0_p = w_p ;
1143    const Double *ti_p = tint.getStorage( b1 ) ;
1144    const Double *w1_p = ti_p ;
1145    for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1146      for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1147        *w0_p = *w1_p ;
1148        w0_p++ ;
1149      }
1150      w1_p++ ;
1151    }
1152    w.putStorage( w_p, b0 ) ;
1153    tint.freeStorage( ti_p, b1 ) ;
1154  }
1155  else if ( wtype_.compare( "TSYS" ) == 0 ) {
1156    Bool b0 ;
1157    Float *w_p = w.getStorage( b0 ) ;
1158    Float *w0_p = w_p ;
1159    for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1160      for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1161        Float temp = *w0_p ;
1162        *w0_p = 1.0 / ( temp * temp ) ;
1163        w0_p++ ;
1164      }
1165    }
1166    w.putStorage( w_p, b0 ) ;
1167  }
1168  else if ( wtype_.compare( "TINTSYS" ) == 0 ) {
1169    Bool b0, b1 ;
1170    Float *w_p = w.getStorage( b0 ) ;
1171    Float *w0_p = w_p ;
1172    const Double *ti_p = tint.getStorage( b1 ) ;
1173    const Double *w1_p = ti_p ;
1174    for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1175      Float interval = *w1_p ;
1176      for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1177        Float temp = *w0_p ;
1178        *w0_p = interval / ( temp * temp ) ;
1179        w0_p++ ;
1180      }
1181      w1_p++ ;
1182    }
1183    w.putStorage( w_p, b0 ) ;
1184    tint.freeStorage( ti_p, b1 ) ;
1185  }
1186  else {
1187    //LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
1188    //os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
1189    w = 1.0 ;
1190  }
1191}
1192
1193void STGrid::toInt( Array<uChar> &u, Array<Int> &v )
1194{
1195  uInt len = u.nelements() ;
1196  Int *int_p = new Int[len] ;
1197  Bool deleteIt ;
1198  const uChar *data_p = u.getStorage( deleteIt ) ;
1199  Int *i_p = int_p ;
1200  const uChar *u_p = data_p ;
1201  for ( uInt i = 0 ; i < len ; i++ ) {
1202    *i_p = ( *u_p == 0 ) ? 0 : 1 ;
1203    i_p++ ;
1204    u_p++ ;
1205  }
1206  u.freeStorage( data_p, deleteIt ) ;
1207  v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
1208}
1209
1210void STGrid::toInt( Array<uInt> &u, Array<Int> &v )
1211{
1212  uInt len = u.nelements() ;
1213  Int *int_p = new Int[len] ;
1214  Bool deleteIt ;
1215  const uInt *data_p = u.getStorage( deleteIt ) ;
1216  Int *i_p = int_p ;
1217  const uInt *u_p = data_p ;
1218  for ( uInt i = 0 ; i < len ; i++ ) {
1219    *i_p = ( *u_p == 0 ) ? 0 : 1 ;
1220    i_p++ ;
1221    u_p++ ;
1222  }
1223  u.freeStorage( data_p, deleteIt ) ;
1224  v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
1225}
1226
1227void STGrid::toPixel( Array<Double> &world, Array<Double> &pixel )
1228{
1229  // gridding will be done on (nx_+2*convSupport_) x (ny_+2*convSupport_)
1230  // grid plane to avoid unexpected behavior on grid edge
1231  Block<Double> pixc( 2 ) ;
1232  pixc[0] = Double( nx_-1 ) * 0.5 ;
1233  pixc[1] = Double( ny_-1 ) * 0.5 ;
1234//   pixc[0] = Double( nx_+2*convSupport_-1 ) * 0.5 ;
1235//   pixc[1] = Double( ny_+2*convSupport_-1 ) * 0.5 ;
1236  uInt nrow = world.shape()[1] ;
1237  Bool bw, bp ;
1238  const Double *w_p = world.getStorage( bw ) ;
1239  Double *p_p = pixel.getStorage( bp ) ;
1240  const Double *ww_p = w_p ;
1241  Double *wp_p = p_p ;
1242  for ( uInt i = 0 ; i < nrow ; i++ ) {
1243    *wp_p = pixc[0] + ( *ww_p - center_[0] ) / cellx_ ;
1244    wp_p++ ;
1245    ww_p++ ;
1246    *wp_p = pixc[1] + ( *ww_p - center_[1] ) / celly_ ;
1247    wp_p++ ;
1248    ww_p++ ;
1249  }
1250  world.freeStorage( w_p, bw ) ;
1251  pixel.putStorage( p_p, bp ) ; 
1252//   String gridfile = "grid."+convType_+"."+String::toString(convSupport_)+".dat" ;
1253//   ofstream ofs( gridfile.c_str(), ios::out ) ;
1254//   ofs << "center " << center_(0) << " " << pixc(0)
1255//       << " " << center_(1) << " " << pixc(1) << endl ;
1256//   for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
1257//     ofs << irow ;
1258//     for ( uInt i = 0 ; i < 2 ; i++ ) {
1259//       ofs << " " << world(i, irow) << " " << pixel(i, irow) ;
1260//     }
1261//     ofs << endl ;
1262//   }
1263//   ofs.close() ;
1264}
1265
1266void STGrid::boxFunc( Vector<Float> &convFunc, Int &convSize )
1267{
1268  convFunc = 0.0 ;
1269  for ( Int i = 0 ; i < convSize/2 ; i++ )
1270    convFunc(i) = 1.0 ;
1271}
1272
1273#define NEED_UNDERSCORES
1274#if defined(NEED_UNDERSCORES)
1275#define grdsf grdsf_
1276#endif
1277extern "C" {
1278   void grdsf(Double*, Double*);
1279}
1280void STGrid::spheroidalFunc( Vector<Float> &convFunc )
1281{
1282  convFunc = 0.0 ;
1283  for ( Int i = 0 ; i < convSampling_*convSupport_ ; i++ ) {
1284    Double nu = Double(i) / Double(convSupport_*convSampling_) ;
1285    Double val ;
1286    grdsf( &nu, &val ) ;
1287    convFunc(i) = ( 1.0 - nu * nu ) * val ;
1288  }
1289}
1290
1291void STGrid::gaussFunc( Vector<Float> &convFunc )
1292{
1293  convFunc = 0.0 ;
1294  // HWHM of the Gaussian is convSupport_ / 4
1295  // To take into account Gaussian tail, kernel cutoff is set to 4 * HWHM
1296  Int len = convSampling_ * convSupport_ ;
1297  Double hwhm = len * 0.25 ;
1298  for ( Int i = 0 ; i < len ; i++ ) {
1299    Double val = Double(i) / hwhm ;
1300    convFunc(i) = exp( -log(2)*val*val ) ;
1301  }
1302}
1303
1304void STGrid::pbFunc( Vector<Float> &convFunc )
1305{
1306  convFunc = 0.0 ;
1307}
1308
1309void STGrid::setConvFunc( Vector<Float> &convFunc )
1310{
1311  convSupport_ = userSupport_ ;
1312  if ( convType_ == "BOX" ) {
1313    if ( convSupport_ < 0 )
1314      convSupport_ = 0 ;
1315    Int convSize = convSampling_ * ( 2 * convSupport_ + 2 )  ;
1316    convFunc.resize( convSize ) ;
1317    boxFunc( convFunc, convSize ) ;
1318  }
1319  else if ( convType_ == "SF" ) {
1320    if ( convSupport_ < 0 )
1321      convSupport_ = 3 ;
1322    Int convSize = convSampling_ * ( 2 * convSupport_ + 2 )  ;
1323    convFunc.resize( convSize ) ;
1324    spheroidalFunc( convFunc ) ;
1325  }
1326  else if ( convType_ == "GAUSS" ) {
1327    // to take into account Gaussian tail
1328    if ( convSupport_ < 0 )
1329      convSupport_ = 12 ; // 3 * 4
1330    else {
1331      convSupport_ = userSupport_ * 4 ;
1332    }
1333    Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
1334    convFunc.resize( convSize ) ;
1335    gaussFunc( convFunc ) ;
1336  }
1337  else if ( convType_ == "PB" ) {
1338    if ( convSupport_ < 0 )
1339      convSupport_ = 0 ;
1340    pbFunc( convFunc ) ;
1341  }
1342  else {
1343    throw AipsError( "Unsupported convolution function" ) ;
1344  }
1345}
1346
1347string STGrid::saveData( string outfile )
1348{
1349  LogIO os( LogOrigin("STGrid", "saveData", WHERE) ) ;
1350  double t0, t1 ;
1351  t0 = mathutil::gettimeofday_sec() ;
1352
1353  //Int polno = 0 ;
1354  String outfile_ ;
1355  if ( outfile.size() == 0 ) {
1356    if ( infileList_[0].lastchar() == '/' ) {
1357      outfile_ = infileList_[0].substr( 0, infileList_[0].size()-1 ) ;
1358    }
1359    else {
1360      outfile_ = infileList_[0] ;
1361    }
1362    outfile_ += ".grid" ;
1363  }
1364  else {
1365    outfile_ = outfile ;
1366  }
1367  Table tab ;
1368  prepareTable( tab, outfile_ ) ;
1369  IPosition dshape = data_.shape() ;
1370  Int nrow = nx_ * ny_ * npol_ ;
1371  tab.rwKeywordSet().define( "nPol", npol_ ) ;
1372  tab.addRow( nrow ) ;
1373  Vector<Double> cpix( 2 ) ;
1374  cpix(0) = Double( nx_ - 1 ) * 0.5 ;
1375  cpix(1) = Double( ny_ - 1 ) * 0.5 ;
1376  Vector<Double> dir( 2 ) ;
1377  ArrayColumn<Double> directionCol( tab, "DIRECTION" ) ;
1378  ArrayColumn<Float> spectraCol( tab, "SPECTRA" ) ;
1379  ScalarColumn<uInt> polnoCol( tab, "POLNO" ) ;
1380  Int irow = 0 ;
1381  Vector<Float> sp( nchan_ ) ;
1382  Bool bsp, bdata ;
1383  const Float *data_p = data_.getStorage( bdata ) ;
1384  Float *wsp_p, *sp_p ;
1385  const Float *wdata_p = data_p ;
1386  long step = nx_ * ny_ * npol_ ;
1387  long offset ;
1388  for ( Int iy = 0 ; iy < ny_ ; iy++ ) {
1389    dir(1) = center_(1) - ( cpix(1) - (Double)iy ) * celly_ ;
1390    for ( Int ix = 0 ; ix < nx_ ; ix++ ) {
1391      dir(0) = center_(0) - ( cpix(0) - (Double)ix ) * cellx_ ;
1392      for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
1393        offset = ix + iy * nx_ + ipol * nx_ * ny_ ;
1394        //os << "offset = " << offset << LogIO::POST ;
1395        sp_p = sp.getStorage( bsp ) ;
1396        wsp_p = sp_p ;
1397        wdata_p = data_p + offset ;
1398        for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
1399          *wsp_p = *wdata_p ;
1400          wsp_p++ ;
1401          wdata_p += step ;
1402        }
1403        sp.putStorage( sp_p, bsp ) ;
1404        spectraCol.put( irow, sp ) ;
1405        directionCol.put( irow, dir ) ;
1406        polnoCol.put( irow, pollist_[ipol] ) ;
1407        irow++ ;
1408      }
1409    }
1410  }
1411  data_.freeStorage( data_p, bdata ) ;
1412
1413  t1 = mathutil::gettimeofday_sec() ;
1414  os << "saveData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
1415
1416  return outfile_ ;
1417}
1418
1419void STGrid::prepareTable( Table &tab, String &name )
1420{
1421  Table t( infileList_[0], Table::Old ) ;
1422  t.deepCopy( name, Table::New, False, t.endianFormat(), True ) ;
1423  tab = Table( name, Table::Update ) ;
1424}
1425
1426}
Note: See TracBrowser for help on using the repository browser.