source: trunk/src/STGrid.cpp@ 2631

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

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: test_sdgrid

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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