source: trunk/src/STGrid.cpp@ 2669

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

New Development: No

JIRA Issue: Yes CAS-4429

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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