source: trunk/src/STGrid.cpp@ 2671

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

New Development: No

JIRA Issue: Yes CAS-4429

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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

Also, defined helper function that plots current grid function
against radius in pixel.


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