source: trunk/src/STGrid.cpp@ 3129

Last change on this file since 3129 was 3106, checked in by Takeshi Nakazato, 8 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes/No

Interface Changes: Yes/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...


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