source: trunk/src/STGrid.cpp@ 3026

Last change on this file since 3026 was 3015, checked in by Takeshi Nakazato, 10 years ago

New Development: No

JIRA Issue: Yes CAS-4362

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

Fix for the issue on sdgrid that resulting map extent expands or shrinks when cell size is calculated by the task.
A cause of the issue is that the task always takes 10% margin when it calculates cell size automatically. In this fix, margin is set to 1 pixel.
File size: 63.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 = " << userSupport_ << 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 : common(common), self(obj), 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 : common(common), self(obj), 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 flag_.resize( gshape ) ;
664 flag_ = (uChar)0;
665 STCommonData common = STCommonData(gshape, data_);
666 common.gnx = gnx ;
667 common.gny = gny ;
668
669 // parameters for gridding
670 Int *chanMap = new Int[nchan_] ;
671 for ( Int i = 0 ; i < nchan_ ; i++ ) {
672 chanMap[i] = i ;
673 }
674 common.chanMap = chanMap;
675
676 // convolution kernel
677 t0 = mathutil::gettimeofday_sec() ;
678 setConvFunc( common.convFunc ) ;
679 t1 = mathutil::gettimeofday_sec() ;
680 os << LogIO::DEBUGGING << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
681
682 // for performance check
683 eGetData_ = 0.0 ;
684 eToPixel_ = 0.0 ;
685 eGGridSD_ = 0.0 ;
686 double eInitPol = 0.0 ;
687
688 for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) {
689 initTable( ifile ) ;
690
691 os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ;
692 Broker broker = Broker(produceChunk, consumeChunk);
693 for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
694 t0 = mathutil::gettimeofday_sec() ;
695 initPol( ipol ) ; // set ptab_ and attach()
696 t1 = mathutil::gettimeofday_sec() ;
697 eInitPol += t1-t0 ;
698
699 STContext context(this, common, ipol);
700
701 os << "start pol " << ipol << LogIO::POST ;
702
703 nprocessed_ = 0 ;
704#if 1
705 broker.runProducerAsMasterThread(&context, DO_AHEAD);
706#else
707 for (;;) {
708 bool produced = produceChunk(&context);
709 if (! produced) {
710 break;
711 }
712 consumeChunk(&context);
713 }
714#endif
715
716 os << "end pol " << ipol << LogIO::POST ;
717
718 }
719 os << "end table " << ifile << LogIO::POST ;
720 }
721 os << LogIO::DEBUGGING << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ;
722 os << LogIO::DEBUGGING << "getData: elapsed time is " << eGetData_-eToInt-eGetWeight << " sec." << LogIO::POST ;
723 os << LogIO::DEBUGGING << "toPixel: elapsed time is " << eToPixel_ << " sec." << LogIO::POST ;
724 os << LogIO::DEBUGGING << "ggridsd: elapsed time is " << eGGridSD_ << " sec." << LogIO::POST ;
725 os << LogIO::DEBUGGING << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
726 os << LogIO::DEBUGGING << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
727
728 delete chanMap ;
729
730 // set data
731 setData( common.gdataArrC, common.gwgtArr ) ;
732
733}
734
735void STGrid::consumeChunkWithClipping(void *ctx) throw(PCException)
736{
737 STContextWithClipping &context = *(STContextWithClipping *)ctx;
738 STGChunk *chunk = NULL;
739 try {
740 context.queue.lock();
741 chunk = context.queue.get();
742 context.queue.unlock();
743 } catch (FullException &e) {
744 context.queue.unlock();
745 // TODO: log error
746 throw PCException();
747 }
748
749 double t0, t1 ;
750 // world -> pixel
751 Array<Double> xypos( context.self->dshape_ ) ;
752 t0 = mathutil::gettimeofday_sec() ;
753 context.self->toPixel( chunk->direction, xypos ) ;
754 t1 = mathutil::gettimeofday_sec() ;
755 context.self->eToPixel_ += t1-t0 ;
756
757 // call ggridsd
758 Int nvispol = 1 ;
759 Int irow = -1 ;
760 t0 = mathutil::gettimeofday_sec() ;
761 context.self->call_ggridsd2( xypos,
762 chunk->spectra,
763 nvispol,
764 context.self->nchan_,
765 chunk->flagtra,
766 chunk->rflag,
767 chunk->weight,
768 chunk->nrow,
769 irow,
770 context.common.gdataArrC,
771 context.common.gwgtArr,
772 context.common.npoints,
773 context.common.clipMin,
774 context.common.clipWMin,
775 context.common.clipCMin,
776 context.common.clipMax,
777 context.common.clipWMax,
778 context.common.clipCMax,
779 context.common.gnx,
780 context.common.gny,
781 context.self->npol_,
782 context.self->nchan_,
783 context.self->convSupport_,
784 context.self->convSampling_,
785 context.common.convFunc,
786 context.common.chanMap,
787 (Int*)&context.pol ) ;
788 t1 = mathutil::gettimeofday_sec() ;
789 context.self->eGGridSD_ += t1-t0 ;
790
791 delete chunk;
792}
793
794void STGrid::gridPerRowWithClipping()
795{
796 LogIO os( LogOrigin("STGrid", "gridPerRowWithClipping", WHERE) ) ;
797 double t0, t1 ;
798
799
800 // grid data
801 // Extend grid plane with convSupport_
802 // Int gnx = nx_+convSupport_*2 ;
803 // Int gny = ny_+convSupport_*2 ;
804 Int gnx = nx_;
805 Int gny = ny_;
806
807 IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
808 IPosition pshape( 3, gnx, gny, npol_ ) ;
809 // 2011/12/20 TN
810 // data_ and gwgtArr share storage
811 data_.resize( gshape ) ;
812 data_ = 0.0 ;
813 flag_.resize( gshape ) ;
814 flag_ = (uChar)0;
815 STCommonDataWithClipping common = STCommonDataWithClipping( gshape,
816 pshape,
817 data_ ) ;
818 common.gnx = gnx ;
819 common.gny = gny ;
820
821 // parameters for gridding
822 Int *chanMap = new Int[nchan_] ;
823 for ( Int i = 0 ; i < nchan_ ; i++ ) {
824 chanMap[i] = i ;
825 }
826 common.chanMap = chanMap;
827
828 // convolution kernel
829 t0 = mathutil::gettimeofday_sec() ;
830 setConvFunc( common.convFunc ) ;
831 t1 = mathutil::gettimeofday_sec() ;
832 os << LogIO::DEBUGGING << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
833
834 // for performance check
835 eGetData_ = 0.0 ;
836 eToPixel_ = 0.0 ;
837 eGGridSD_ = 0.0 ;
838 double eInitPol = 0.0 ;
839
840 for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) {
841 initTable( ifile ) ;
842
843 os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ;
844 Broker broker = Broker(produceChunk, consumeChunkWithClipping);
845 for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
846 t0 = mathutil::gettimeofday_sec() ;
847 initPol( ipol ) ; // set ptab_ and attach()
848 t1 = mathutil::gettimeofday_sec() ;
849 eInitPol += t1-t0 ;
850
851 STContextWithClipping context(this, common, ipol);
852
853 os << "start pol " << ipol << LogIO::POST ;
854
855 nprocessed_ = 0 ;
856#if 1
857 broker.runProducerAsMasterThread(&context, DO_AHEAD);
858#else
859 for (;;) {
860 bool produced = produceChunk(&context);
861 if (! produced) {
862 break;
863 }
864 consumeChunkWithClipping(&context);
865 }
866#endif
867
868 os << "end pol " << ipol << LogIO::POST ;
869
870 }
871 os << "end table " << ifile << LogIO::POST ;
872 }
873 os << LogIO::DEBUGGING << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ;
874 os << LogIO::DEBUGGING << "getData: elapsed time is " << eGetData_-eToInt-eGetWeight << " sec." << LogIO::POST ;
875 os << LogIO::DEBUGGING << "toPixel: elapsed time is " << eToPixel_ << " sec." << LogIO::POST ;
876 os << LogIO::DEBUGGING << "ggridsd2: elapsed time is " << eGGridSD_ << " sec." << LogIO::POST ;
877 os << LogIO::DEBUGGING << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
878 os << LogIO::DEBUGGING << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
879
880 delete chanMap ;
881
882 // clip min and max in each grid
883// os << "BEFORE CLIPPING" << LogIO::POST ;
884// os << "gdataArrC=" << common.gdataArrC << LogIO::POST ;
885// os << "gwgtArr=" << common.gwgtArr << LogIO::POST ;
886 t0 = mathutil::gettimeofday_sec() ;
887 clipMinMax( common.gdataArrC,
888 common.gwgtArr,
889 common.npoints,
890 common.clipMin,
891 common.clipWMin,
892 common.clipCMin,
893 common.clipMax,
894 common.clipWMax,
895 common.clipCMax ) ;
896 t1 = mathutil::gettimeofday_sec() ;
897 os << LogIO::DEBUGGING << "clipMinMax: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
898// os << "AFTER CLIPPING" << LogIO::POST ;
899// os << "gdataArrC=" << common.gdataArrC << LogIO::POST ;
900// os << "gwgtArr=" << common.gwgtArr << LogIO::POST ;
901
902 // set data
903 setData( common.gdataArrC, common.gwgtArr ) ;
904
905}
906
907void STGrid::clipMinMax( Array<Complex> &grid,
908 Array<Float> &weight,
909 Array<Int> &npoints,
910 Array<Complex> &clipmin,
911 Array<Float> &clipwmin,
912 Array<Float> &clipcmin,
913 Array<Complex> &clipmax,
914 Array<Float> &clipwmax,
915 Array<Float> &clipcmax )
916{
917 //LogIO os( LogOrigin("STGrid","clipMinMax",WHERE) ) ;
918
919 // prepare pointers
920 Bool delG, delW, delNP, delCMin, delCWMin, delCCMin, delCMax, delCWMax, delCCMax ;
921 Complex *grid_p = grid.getStorage( delG ) ;
922 Float *wgt_p = weight.getStorage( delW ) ;
923 const Int *npts_p = npoints.getStorage( delNP ) ;
924 const Complex *cmin_p = clipmin.getStorage( delCMin ) ;
925 const Float *cwmin_p = clipwmin.getStorage( delCWMin ) ;
926 const Float *ccmin_p = clipcmin.getStorage( delCCMin ) ;
927 const Complex *cmax_p = clipmax.getStorage( delCMax ) ;
928 const Float *cwmax_p = clipwmax.getStorage( delCWMax ) ;
929 const Float *ccmax_p = clipcmax.getStorage( delCCMax ) ;
930
931 const IPosition &gshape = grid.shape() ;
932 long offset = gshape[0] * gshape[1] * gshape[2] ; // nx * ny * npol
933 Int nchan = gshape[3] ;
934 long origin = nchan * offset ;
935 for ( long i = 0 ; i < offset ; i++ ) {
936 if ( *npts_p > 2 ) {
937 for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
938 // clip minimum and maximum
939 *grid_p -= (*cmin_p)*(*cwmin_p)*(*ccmin_p)
940 + (*cmax_p)*(*cwmax_p)*(*ccmax_p) ;
941 *wgt_p -= (*cwmin_p)*(*ccmin_p)
942 + (*cwmax_p)*(*ccmax_p) ;
943
944 grid_p += offset ;
945 wgt_p += offset ;
946 cmin_p += offset ;
947 cwmin_p += offset ;
948 ccmin_p += offset ;
949 cmax_p += offset ;
950 cwmax_p += offset ;
951 ccmax_p += offset ;
952 }
953 grid_p -= origin ;
954 wgt_p -= origin ;
955 cmin_p -= origin ;
956 cwmin_p -= origin ;
957 ccmin_p -= origin ;
958 cmax_p -= origin ;
959 cwmax_p -= origin ;
960 ccmax_p -= origin ;
961 }
962 grid_p++ ;
963 wgt_p++ ;
964 npts_p++ ;
965 cmin_p++ ;
966 cwmin_p++ ;
967 ccmin_p++ ;
968 cmax_p++ ;
969 cwmax_p++ ;
970 ccmax_p++ ;
971 }
972 grid_p -= offset ;
973 wgt_p -= offset ;
974 npts_p -= offset ;
975 cmin_p -= offset ;
976 cwmin_p -= offset ;
977 ccmin_p -= offset ;
978 cmax_p -= offset ;
979 cwmax_p -= offset ;
980 ccmax_p -= offset ;
981
982 // finalization
983 grid.putStorage( grid_p, delG ) ;
984 weight.putStorage( wgt_p, delW ) ;
985 npoints.freeStorage( npts_p, delNP ) ;
986 clipmin.freeStorage( cmin_p, delCMin ) ;
987 clipwmin.freeStorage( cwmin_p, delCWMin ) ;
988 clipcmin.freeStorage( ccmin_p, delCCMin ) ;
989 clipmax.freeStorage( cmax_p, delCMax ) ;
990 clipwmax.freeStorage( cwmax_p, delCWMax ) ;
991 clipcmax.freeStorage( ccmax_p, delCCMax ) ;
992}
993
994void STGrid::initPol( Int ipol )
995{
996 LogIO os( LogOrigin("STGrid","initPol",WHERE) ) ;
997 if ( npolOrg_ == 1 ) {
998 os << "single polarization data." << LogIO::POST ;
999 ptab_ = tab_ ;
1000 }
1001 else
1002 ptab_ = tab_( tab_.col("POLNO") == pollist_[ipol] ) ;
1003
1004 attach( ptab_ ) ;
1005}
1006
1007void STGrid::initTable( uInt idx )
1008{
1009 tab_ = tableList_[idx] ;
1010 nrow_ = rows_[idx] ;
1011 updateChunkShape() ;
1012}
1013
1014void STGrid::setData( Array<Complex> &gdata,
1015 Array<Float> &gwgt )
1016{
1017 // 2011/12/20 TN
1018 // gwgt and data_ share storage
1019 LogIO os( LogOrigin("STGrid","setData",WHERE) ) ;
1020 double t0, t1 ;
1021 t0 = mathutil::gettimeofday_sec() ;
1022 uInt len = data_.nelements() ;
1023 Bool b1, b2, b3 ;
1024 const Complex *gdata_p = gdata.getStorage( b1 ) ;
1025 Float *gwgt_p = gwgt.getStorage( b2 ) ; // storage shared with data_
1026 uChar *gflg_p = flag_.getStorage( b3 ) ;
1027 for ( uInt i = 0 ; i < len ; i++ ) {
1028 if (gwgt_p[i] > 0.0) {
1029 gwgt_p[i] = (gdata_p[i]).real() / gwgt_p[i];
1030 gflg_p[i] = (uChar)0;
1031 }
1032 else {
1033 gflg_p[i] = (uChar)1;
1034 }
1035 }
1036 gdata.freeStorage( gdata_p, b1 ) ;
1037 data_.putStorage( gwgt_p, b2 ) ;
1038 flag_.putStorage( gflg_p, b3 ) ;
1039 t1 = mathutil::gettimeofday_sec() ;
1040 os << LogIO::DEBUGGING << "setData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
1041}
1042
1043void STGrid::setupGrid()
1044{
1045 Double xmin,xmax,ymin,ymax ;
1046 mapExtent( xmin, xmax, ymin, ymax ) ;
1047
1048 setupGrid( nxUI_, nyUI_, cellxUI_, cellyUI_,
1049 xmin, xmax, ymin, ymax, centerUI_ ) ;
1050}
1051
1052void STGrid::setupGrid( Int &nx,
1053 Int &ny,
1054 String &cellx,
1055 String &celly,
1056 Double &xmin,
1057 Double &xmax,
1058 Double &ymin,
1059 Double &ymax,
1060 String &center )
1061{
1062 LogIO os( LogOrigin("STGrid","setupGrid",WHERE) ) ;
1063 //cout << "nx=" << nx << ", ny=" << ny << endl ;
1064
1065 // center position
1066 if ( center.size() == 0 ) {
1067 center_(0) = 0.5 * ( xmin + xmax ) ;
1068 center_(1) = 0.5 * ( ymin + ymax ) ;
1069 }
1070 else {
1071 String::size_type pos0 = center.find( " " ) ;
1072 if ( pos0 == String::npos ) {
1073 throw AipsError( "bad string format in parameter center" ) ;
1074 }
1075 String::size_type pos1 = center.find( " ", pos0+1 ) ;
1076 String typestr, xstr, ystr ;
1077 if ( pos1 != String::npos ) {
1078 typestr = center.substr( 0, pos0 ) ;
1079 xstr = center.substr( pos0+1, pos1-pos0 ) ;
1080 ystr = center.substr( pos1+1 ) ;
1081 // todo: convert to J2000 (or direction ref for DIRECTION column)
1082 }
1083 else {
1084 typestr = "J2000" ;
1085 xstr = center.substr( 0, pos0 ) ;
1086 ystr = center.substr( pos0+1 ) ;
1087 }
1088 QuantumHolder qh ;
1089 String err ;
1090 qh.fromString( err, xstr ) ;
1091 Quantum<Double> xcen = qh.asQuantumDouble() ;
1092 qh.fromString( err, ystr ) ;
1093 Quantum<Double> ycen = qh.asQuantumDouble() ;
1094 center_(0) = xcen.getValue( "rad" ) ;
1095 center_(1) = ycen.getValue( "rad" ) ;
1096 double base = 0.5 * (xmin + xmax) ;
1097 int maxrotate = 1 ;
1098 int nelem = 2 * maxrotate + 1 ;
1099 double *sep = new double[nelem] ;
1100 for ( int i = 0 ; i < nelem ; i++ )
1101 sep[i] = abs(base - center_[0] - (i-maxrotate) * C::_2pi) ;
1102// os << "sep[0]=" << sep[0] << endl
1103// << "sep[1]=" << sep[1] << endl
1104// << "sep[2]=" << sep[2] << LogIO::POST ;
1105 int idx = 0 ;
1106 base = sep[0] ;
1107 int nrotate = 0 ;
1108 while ( idx < nelem ) {
1109 if ( base > sep[idx] ) {
1110 base = sep[idx] ;
1111 nrotate = idx ;
1112 }
1113 idx++ ;
1114 }
1115 delete sep ;
1116 nrotate -= maxrotate ;
1117// os << "nrotate = " << nrotate << LogIO::POST ;
1118 center_[0] += nrotate * C::_2pi ;
1119 }
1120// os << "xmin=" << xmin << LogIO::POST ;
1121// os << "center_=" << center_ << LogIO::POST ;
1122
1123 nx_ = nx ;
1124 ny_ = ny ;
1125 if ( nx < 0 && ny > 0 ) {
1126 nx_ = ny ;
1127 ny_ = ny ;
1128 }
1129 if ( ny < 0 && nx > 0 ) {
1130 nx_ = nx ;
1131 ny_ = nx ;
1132 }
1133
1134 //Double wx = xmax - xmin ;
1135 //Double wy = ymax - ymin ;
1136 Double wx = max( abs(xmax-center_(0)), abs(xmin-center_(0)) ) * 2 ;
1137 Double wy = max( abs(ymax-center_(1)), abs(ymin-center_(1)) ) * 2 ;
1138
1139 Quantum<Double> qcellx ;
1140 Quantum<Double> qcelly ;
1141 //cout << "nx_ = " << nx_ << ", ny_ = " << ny_ << endl ;
1142 if ( cellx.size() != 0 && celly.size() != 0 ) {
1143 readQuantity( qcellx, cellx ) ;
1144 readQuantity( qcelly, celly ) ;
1145 }
1146 else if ( celly.size() != 0 ) {
1147 os << "Using celly to x-axis..." << LogIO::POST ;
1148 readQuantity( qcelly, celly ) ;
1149 qcellx = qcelly ;
1150 }
1151 else if ( cellx.size() != 0 ) {
1152 os << "Using cellx to y-axis..." << LogIO::POST ;
1153 readQuantity( qcellx, cellx ) ;
1154 qcelly = qcellx ;
1155 }
1156 else {
1157 if ( nx_ < 0 ) {
1158 os << "No user preference in grid setting. Using default..." << LogIO::POST ;
1159 readQuantity( qcellx, "1.0arcmin" ) ;
1160 qcelly = qcellx ;
1161 }
1162 else {
1163 if ( wx == 0.0 ) {
1164 os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
1165 wx = 0.00290888 ;
1166 }
1167 if ( wy == 0.0 ) {
1168 os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
1169 wy = 0.00290888 ;
1170 }
1171 if (nx_ > 1) {
1172 qcellx = Quantum<Double>(wx / (nx_ - 1) * cos(center_[1]), "rad");
1173 }
1174 else {
1175 qcellx = Quantum<Double>( 1.1f * wx / nx_ * cos(center_[1]), "rad" );
1176 }
1177 if (ny_ > 1) {
1178 qcelly = Quantum<Double>(wy / (ny_ - 1), "rad");
1179 }
1180 else {
1181 qcelly = Quantum<Double>( 1.1f * wy / ny_, "rad" );
1182 }
1183 }
1184 }
1185 cellx_ = qcellx.getValue( "rad" ) ;
1186 celly_ = qcelly.getValue( "rad" ) ;
1187 //os << "cellx_=" << cellx_ << ", celly_=" << celly_ << ", cos("<<center_(1)<<")=" << cos(center_(1)) << LogIO::POST ;
1188 if ( nx_ < 0 ) {
1189 if ( wx == 0.0 ) {
1190 os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
1191 wx = 0.00290888 ;
1192 }
1193 if ( wy == 0.0 ) {
1194 os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
1195 wy = 0.00290888 ;
1196 }
1197 nx_ = Int( ceil( wx/(cellx_/cos(center_[1])) ) ) ;
1198 ny_ = Int( ceil( wy/celly_ ) ) ;
1199 }
1200
1201 // create DirectionCoordinate
1202 Matrix<Double> xform(2,2) ;
1203 xform = 0.0 ;
1204 xform.diagonal() = 1.0 ;
1205 dircoord_ = new DirectionCoordinate(MDirection::J2000,
1206 Projection( Projection::SIN ),
1207 center_[0], center_[1],
1208 -cellx_, celly_,
1209 xform,
1210 0.5*Double(nx_-1),
1211 0.5*Double(ny_-1)) ;
1212}
1213
1214void STGrid::mapExtent( Double &xmin, Double &xmax,
1215 Double &ymin, Double &ymax )
1216{
1217 //LogIO os( LogOrigin("STGrid","mapExtent",WHERE) ) ;
1218 directionCol_.attach( tableList_[0], "DIRECTION" ) ;
1219 Matrix<Double> direction = directionCol_.getColumn() ;
1220 //os << "dirCol.nrow() = " << dirCol.nrow() << LogIO::POST ;
1221 Vector<Double> ra( direction.row(0) ) ;
1222 mathutil::rotateRA( ra ) ;
1223 minMax( xmin, xmax, direction.row( 0 ) ) ;
1224 minMax( ymin, ymax, direction.row( 1 ) ) ;
1225 Double amin, amax, bmin, bmax ;
1226 for ( uInt i = 1 ; i < nfile_ ; i++ ) {
1227 directionCol_.attach( tableList_[i], "DIRECTION" ) ;
1228 direction.assign( directionCol_.getColumn() ) ;
1229 //os << "dirCol.nrow() = " << dirCol.nrow() << LogIO::POST ;
1230 // to make contiguous RA distribution (no 2pi jump)
1231 Vector<Double> ra( direction.row(0) ) ;
1232 mathutil::rotateRA( ra ) ;
1233 minMax( amin, amax, direction.row( 0 ) ) ;
1234 minMax( bmin, bmax, direction.row( 1 ) ) ;
1235 xmin = min( xmin, amin ) ;
1236 xmax = max( xmax, amax ) ;
1237 ymin = min( ymin, bmin ) ;
1238 ymax = max( ymax, bmax ) ;
1239 }
1240 //os << "(xmin,xmax)=(" << xmin << "," << xmax << ")" << LogIO::POST ;
1241 //os << "(ymin,ymax)=(" << ymin << "," << ymax << ")" << LogIO::POST ;
1242}
1243
1244void STGrid::table( Table &tab, uInt i )
1245{
1246 if ( i < nfile_ )
1247 tab = Table( infileList_[i] ) ;
1248}
1249
1250void STGrid::selectData()
1251{
1252 LogIO os( LogOrigin("STGrid","selectData",WHERE) ) ;
1253 Int ifno = ifno_ ;
1254 tableList_.resize( nfile_ ) ;
1255 if ( ifno_ == -1 ) {
1256 //Table taborg( infileList_[0] ) ;
1257 Table taborg ;
1258 table( taborg, 0 ) ;
1259 ROScalarColumn<uInt> ifnoCol( taborg, "IFNO" ) ;
1260 ifno_ = ifnoCol( 0 ) ;
1261 os << LogIO::WARN
1262 << "IFNO is not given. Using default IFNO: " << ifno_ << LogIO::POST ;
1263 }
1264 for ( uInt i = 0 ; i < nfile_ ; i++ ) {
1265 //Table taborg( infileList_[i] ) ;
1266 Table taborg ;
1267 table( taborg, i ) ;
1268 TableExprNode node ;
1269 if ( ifno != -1 || isMultiIF( taborg ) ) {
1270 os << "apply selection on IFNO" << LogIO::POST ;
1271 node = taborg.col("IFNO") == ifno_ ;
1272 }
1273 if ( scanlist_.size() > 0 ) {
1274 os << "apply selection on SCANNO" << LogIO::POST ;
1275 node = node && taborg.col("SCANNO").in( scanlist_ ) ;
1276 }
1277 if ( node.isNull() ) {
1278 tableList_[i] = taborg ;
1279 }
1280 else {
1281 tableList_[i] = taborg( node ) ;
1282 }
1283 os << LogIO::DEBUGGING << "tableList_[" << i << "].nrow()=" << tableList_[i].nrow() << LogIO::POST ;
1284 if ( tableList_[i].nrow() == 0 ) {
1285 os << LogIO::SEVERE
1286 << "No corresponding rows for given selection: IFNO " << ifno_ ;
1287 if ( scanlist_.size() > 0 )
1288 os << " SCANNO " << scanlist_ ;
1289 os << LogIO::EXCEPTION ;
1290 }
1291 }
1292}
1293
1294Bool STGrid::isMultiIF( Table &tab )
1295{
1296 ROScalarColumn<uInt> ifnoCol( tab, "IFNO" ) ;
1297 Vector<uInt> ifnos = ifnoCol.getColumn() ;
1298 return anyNE( ifnos, ifnos[0] ) ;
1299}
1300
1301void STGrid::attach( Table &tab )
1302{
1303 // attach to table
1304 spectraCol_.attach( tab, "SPECTRA" ) ;
1305 flagtraCol_.attach( tab, "FLAGTRA" ) ;
1306 directionCol_.attach( tab, "DIRECTION" ) ;
1307 flagRowCol_.attach( tab, "FLAGROW" ) ;
1308 tsysCol_.attach( tab, "TSYS" ) ;
1309 intervalCol_.attach( tab, "INTERVAL" ) ;
1310}
1311
1312Int STGrid::getDataChunk(
1313 IPosition const &wshape,
1314 IPosition const &vshape,
1315 IPosition const &/* dshape */,
1316 Array<Complex> &spectra,
1317 Array<Double> &direction,
1318 Array<Int> &flagtra,
1319 Array<Int> &rflag,
1320 Array<Float> &weight )
1321{
1322 LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
1323
1324 Array<Float> spectraF_(wshape);
1325 Array<uChar> flagtraUC_(wshape);
1326 Array<uInt> rflagUI_(vshape);
1327 Int nrow = getDataChunk( spectraF_, direction, flagtraUC_, rflagUI_, weight ) ;
1328 if ( nrow < nchunk_ ) {
1329 spectra.resize( spectraF_.shape() ) ;
1330 flagtra.resize( flagtraUC_.shape() ) ;
1331 rflag.resize( rflagUI_.shape() ) ;
1332 }
1333 double t0, t1 ;
1334 t0 = mathutil::gettimeofday_sec() ;
1335 convertArray( spectra, spectraF_ ) ;
1336 toInt( flagtraUC_, flagtra ) ;
1337 toInt( rflagUI_, rflag ) ;
1338 t1 = mathutil::gettimeofday_sec() ;
1339 eToInt = t1 - t0 ;
1340
1341 return nrow ;
1342}
1343
1344#if 0
1345Int STGrid::getDataChunk( Array<Complex> &spectra,
1346 Array<Double> &direction,
1347 Array<Int> &flagtra,
1348 Array<Int> &rflag,
1349 Array<Float> &weight )
1350{
1351 LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
1352 Int nrow = getDataChunk( spectraF_, direction, flagtraUC_, rflagUI_, weight ) ;
1353 if ( nrow < nchunk_ ) {
1354 spectra.resize( spectraF_.shape() ) ;
1355 flagtra.resize( flagtraUC_.shape() ) ;
1356 rflag.resize( rflagUI_.shape() ) ;
1357 }
1358 double t0, t1 ;
1359 t0 = mathutil::gettimeofday_sec() ;
1360 convertArray( spectra, spectraF_ ) ;
1361 toInt( flagtraUC_, flagtra ) ;
1362 toInt( rflagUI_, rflag ) ;
1363 t1 = mathutil::gettimeofday_sec() ;
1364 eToInt = t1 - t0 ;
1365
1366 return nrow ;
1367}
1368#endif
1369
1370Int STGrid::getDataChunk( Array<Float> &spectra,
1371 Array<Double> &direction,
1372 Array<uChar> &flagtra,
1373 Array<uInt> &rflag,
1374 Array<Float> &weight )
1375{
1376 LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
1377 Int nrow = spectra.shape()[1] ;
1378 Int remainingRow = nrow_ - nprocessed_ ;
1379 if ( remainingRow < nrow ) {
1380 nrow = remainingRow ;
1381 IPosition mshape( 2, nchan_, nrow ) ;
1382 IPosition vshape( 1, nrow ) ;
1383 spectra.resize( mshape ) ;
1384 flagtra.resize( mshape ) ;
1385 direction.resize( IPosition(2,2,nrow) ) ;
1386 rflag.resize( vshape ) ;
1387 weight.resize( mshape ) ;
1388 }
1389 // 2011/12/22 TN
1390 // tsys shares its storage with weight
1391 Array<Float> tsys( weight ) ;
1392 Array<Double> tint( rflag.shape() ) ;
1393
1394 Vector<uInt> rflagVec( rflag ) ;
1395 Vector<Double> tintVec( tint ) ;
1396
1397 RefRows rows( nprocessed_, nprocessed_+nrow-1, 1 ) ;
1398 //os<<LogIO::DEBUGGING<<"nprocessed_="<<nprocessed_<<": rows.nrows()="<<rows.nrows()<<LogIO::POST ;
1399 spectraCol_.getColumnCells( rows, spectra ) ;
1400 flagtraCol_.getColumnCells( rows, flagtra ) ;
1401 directionCol_.getColumnCells( rows, direction ) ;
1402 // to make contiguous RA distribution (no 2pi jump)
1403 Vector<Double> v( Matrix<Double>(direction).row(0) ) ;
1404 mathutil::rotateRA( v ) ;
1405 flagRowCol_.getColumnCells( rows, rflagVec ) ;
1406 intervalCol_.getColumnCells( rows, tintVec ) ;
1407 Vector<Float> tsysTemp = tsysCol_( nprocessed_ ) ;
1408 if ( tsysTemp.nelements() == (uInt)nchan_ )
1409 tsysCol_.getColumnCells( rows, tsys ) ;
1410 else
1411 tsys = tsysTemp[0] ;
1412
1413 double t0,t1 ;
1414 t0 = mathutil::gettimeofday_sec() ;
1415 getWeight( weight, tsys, tint ) ;
1416 t1 = mathutil::gettimeofday_sec() ;
1417 eGetWeight += t1-t0 ;
1418
1419 nprocessed_ += nrow ;
1420
1421 return nrow ;
1422}
1423
1424void STGrid::setupArray()
1425{
1426 LogIO os( LogOrigin("STGrid","setupArray",WHERE) ) ;
1427 ROScalarColumn<uInt> polnoCol( tableList_[0], "POLNO" ) ;
1428 Vector<uInt> pols = polnoCol.getColumn() ;
1429 //os << pols << LogIO::POST ;
1430 Vector<uInt> pollistOrg ;
1431 npolOrg_ = 0 ;
1432 uInt polno ;
1433 for ( uInt i = 0 ; i < polnoCol.nrow() ; i++ ) {
1434 //polno = polnoCol( i ) ;
1435 polno = pols( i ) ;
1436 if ( allNE( pollistOrg, polno ) ) {
1437 pollistOrg.resize( npolOrg_+1, True ) ;
1438 pollistOrg[npolOrg_] = polno ;
1439 npolOrg_++ ;
1440 }
1441 }
1442 if ( pollist_.size() == 0 )
1443 pollist_ = pollistOrg ;
1444 else {
1445 Vector<uInt> newlist ;
1446 uInt newsize = 0 ;
1447 for ( uInt i = 0 ; i < pollist_.size() ; i++ ) {
1448 if ( anyEQ( pollistOrg, pollist_[i] ) ) {
1449 newlist.resize( newsize+1, True ) ;
1450 newlist[newsize] = pollist_[i] ;
1451 newsize++ ;
1452 }
1453 }
1454 pollist_.assign( newlist ) ;
1455 }
1456 npol_ = pollist_.size() ;
1457 if ( npol_ == 0 ) {
1458 os << LogIO::SEVERE << "Empty pollist" << LogIO::EXCEPTION ;
1459 }
1460 rows_.resize( nfile_ ) ;
1461 for ( uInt i = 0 ; i < nfile_ ; i++ ) {
1462 rows_[i] = tableList_[i].nrow() / npolOrg_ ;
1463 //if ( nrow_ < rows_[i] )
1464 // nrow_ = rows_[i] ;
1465 }
1466 flagtraCol_.attach( tableList_[0], "FLAGTRA" ) ;
1467 nchan_ = flagtraCol_( 0 ).nelements() ;
1468// os << "npol_ = " << npol_ << "(" << pollist_ << ")" << endl
1469// << "nchan_ = " << nchan_ << endl
1470// << "nrow_ = " << nrow_ << LogIO::POST ;
1471}
1472
1473void STGrid::getWeight( Array<Float> &w,
1474 Array<Float> &tsys,
1475 Array<Double> &tint )
1476{
1477 LogIO os( LogOrigin("STGrid","getWeight",WHERE) ) ;
1478
1479 // 2011/12/22 TN
1480 // w (weight) and tsys share storage
1481 IPosition refShape = tsys.shape() ;
1482 Int nchan = refShape[0] ;
1483 Int nrow = refShape[1] ;
1484// os << "nchan=" << nchan << ", nrow=" << nrow << LogIO::POST ;
1485// os << "w.shape()=" << w.shape() << endl
1486// << "tsys.shape()=" << tsys.shape() << endl
1487// << "tint.shape()=" << tint.shape() << LogIO::POST ;
1488
1489 // set weight
1490 if ( wtype_.compare( "UNIFORM" ) == 0 ) {
1491 w = 1.0 ;
1492 }
1493 else if ( wtype_.compare( "TINT" ) == 0 ) {
1494 Bool b0, b1 ;
1495 Float *w_p = w.getStorage( b0 ) ;
1496 Float *w0_p = w_p ;
1497 const Double *ti_p = tint.getStorage( b1 ) ;
1498 const Double *w1_p = ti_p ;
1499 for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1500 for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1501 *w0_p = *w1_p ;
1502 w0_p++ ;
1503 }
1504 w1_p++ ;
1505 }
1506 w.putStorage( w_p, b0 ) ;
1507 tint.freeStorage( ti_p, b1 ) ;
1508 }
1509 else if ( wtype_.compare( "TSYS" ) == 0 ) {
1510 Bool b0 ;
1511 Float *w_p = w.getStorage( b0 ) ;
1512 Float *w0_p = w_p ;
1513 for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1514 for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1515 Float temp = *w0_p ;
1516 *w0_p = 1.0 / ( temp * temp ) ;
1517 w0_p++ ;
1518 }
1519 }
1520 w.putStorage( w_p, b0 ) ;
1521 }
1522 else if ( wtype_.compare( "TINTSYS" ) == 0 ) {
1523 Bool b0, b1 ;
1524 Float *w_p = w.getStorage( b0 ) ;
1525 Float *w0_p = w_p ;
1526 const Double *ti_p = tint.getStorage( b1 ) ;
1527 const Double *w1_p = ti_p ;
1528 for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1529 Float interval = *w1_p ;
1530 for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1531 Float temp = *w0_p ;
1532 *w0_p = interval / ( temp * temp ) ;
1533 w0_p++ ;
1534 }
1535 w1_p++ ;
1536 }
1537 w.putStorage( w_p, b0 ) ;
1538 tint.freeStorage( ti_p, b1 ) ;
1539 }
1540 else {
1541 //LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
1542 //os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
1543 w = 1.0 ;
1544 }
1545}
1546
1547void STGrid::toInt( Array<uChar> &u, Array<Int> &v )
1548{
1549 uInt len = u.nelements() ;
1550 Int *int_p = new Int[len] ;
1551 Bool deleteIt ;
1552 const uChar *data_p = u.getStorage( deleteIt ) ;
1553 Int *i_p = int_p ;
1554 const uChar *u_p = data_p ;
1555 for ( uInt i = 0 ; i < len ; i++ ) {
1556 *i_p = ( *u_p == 0 ) ? 0 : 1 ;
1557 i_p++ ;
1558 u_p++ ;
1559 }
1560 u.freeStorage( data_p, deleteIt ) ;
1561 v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
1562}
1563
1564void STGrid::toInt( Array<uInt> &u, Array<Int> &v )
1565{
1566 uInt len = u.nelements() ;
1567 Int *int_p = new Int[len] ;
1568 Bool deleteIt ;
1569 const uInt *data_p = u.getStorage( deleteIt ) ;
1570 Int *i_p = int_p ;
1571 const uInt *u_p = data_p ;
1572 for ( uInt i = 0 ; i < len ; i++ ) {
1573 *i_p = ( *u_p == 0 ) ? 0 : 1 ;
1574 i_p++ ;
1575 u_p++ ;
1576 }
1577 u.freeStorage( data_p, deleteIt ) ;
1578 v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
1579}
1580
1581void STGrid::toPixel( Array<Double> &world, Array<Double> &pixel )
1582{
1583 uInt nrow = world.shape()[1] ;
1584 Bool bw, bp ;
1585 Double *w_p = world.getStorage( bw ) ;
1586 Double *p_p = pixel.getStorage( bp ) ;
1587 Double *ww_p = w_p ;
1588 Double *wp_p = p_p ;
1589 IPosition vshape( 1, 2 ) ;
1590 Vector<Double> _world, _pixel ;
1591 for ( uInt i = 0 ; i < nrow ; i++ ) {
1592 _world.takeStorage( vshape, ww_p, SHARE ) ;
1593 _pixel.takeStorage( vshape, wp_p, SHARE ) ;
1594 dircoord_->toPixel( _pixel, _world ) ;
1595 ww_p += 2 ;
1596 wp_p += 2 ;
1597 }
1598 world.putStorage( w_p, bw ) ;
1599 pixel.putStorage( p_p, bp ) ;
1600}
1601
1602void STGrid::boxFunc( Vector<Float> &convFunc, Int &convSize )
1603{
1604 convFunc = 0.0 ;
1605 for ( Int i = 0 ; i < convSize/2 ; i++ )
1606 convFunc(i) = 1.0 ;
1607}
1608
1609#define NEED_UNDERSCORES
1610#if defined(NEED_UNDERSCORES)
1611#define grdsf grdsf_
1612#define grdgauss grdgauss_
1613#define grdjinc1 grdjinc1_
1614#endif
1615#if defined(USE_CASAPY)
1616extern "C" {
1617 void grdsf(Double*, Double*);
1618 void grdgauss(Double*, Double*, Double*);
1619 void grdjinc1(Double*, Double*, Int*, Double*);
1620}
1621#else
1622extern "C" {
1623 void grdsf(Double*, Double*);
1624}
1625void grdgauss(Double *hwhm, Double *val, Double *out)
1626{
1627 *out = exp(-log(2.0) * (*val / *hwhm) * (*val / *hwhm));
1628}
1629void grdjinc1(Double *c, Double *val, Int *normalize, Double *out)
1630{
1631 // Calculate J_1(x) using approximate formula
1632 Double x = C::pi * *val / *c;
1633 Double ax = fabs(x);
1634 Double ans;
1635 if ( ax < 8.0 ) {
1636 Double y = x * x;
1637 Double ans1 = x * (72362614232.0 + y * (-7895059235.0
1638 + y * (242396853.1 + y * (-2972611.439
1639 + y * (15704.48260 + y * (-30.16036606))))));
1640 Double ans2 = 144725228442.0 + y * (2300535178.0
1641 + y * (18583304.74 + y * (99447.43394
1642 + y * (376.9991397 + y * 1.0))));
1643 ans = ans1 / ans2;
1644 }
1645 else {
1646 Double z = 8.0 / ax;
1647 Double y = z * z;
1648 Double xx = ax - 2.356194491;
1649 Double ans1 = 1.0 + y * (0.183105e-2 + y * (-0.3516396496e-4
1650 + y * (0.2457520174e-5 + y * (-0.240337019e-6))));
1651 Double ans2 = 0.04687499995 + y * (-0.2002690873e-3
1652 + y * (0.8449199096e-5 + y * (-0.88228987e-6
1653 + y * (0.105787412e-6))));
1654 ans = sqrt(0.636619772 / ax) * (cos(xx) * ans1
1655 - z * sin(xx) * ans2);
1656 if (x < 0.0)
1657 ans = -ans;
1658 }
1659
1660 // Then, calculate Jinc
1661 if (x == 0.0) {
1662 *out = 0.5;
1663 }
1664 else {
1665 *out = ans / x;
1666 }
1667
1668 if (*normalize == 1)
1669 *out = *out / 0.5;
1670}
1671#endif
1672void STGrid::spheroidalFunc( Vector<Float> &convFunc )
1673{
1674 convFunc = 0.0 ;
1675 for ( Int i = 0 ; i < convSampling_*convSupport_ ; i++ ) {
1676 Double nu = Double(i) / Double(convSupport_*convSampling_) ;
1677 Double val ;
1678 grdsf( &nu, &val ) ;
1679 convFunc(i) = ( 1.0 - nu * nu ) * val ;
1680 }
1681}
1682
1683void STGrid::gaussFunc( Vector<Float> &convFunc, Double hwhm, Double truncate )
1684{
1685 convFunc = 0.0 ;
1686 Int len = (Int)(truncate*Double(convSampling_)+0.5);
1687 Double out, val;
1688 for ( Int i = 0 ; i < len ; i++ ) {
1689 val = Double(i) / Double(convSampling_) ;
1690 grdgauss(&hwhm, &val, &out);
1691 convFunc(i) = out;
1692 }
1693}
1694
1695void STGrid::gjincFunc( Vector<Float> &convFunc, Double hwhm, Double c, Double truncate )
1696{
1697 convFunc = 0.0;
1698 Double out1, out2, val;
1699 Int normalize = 1;
1700 if (truncate >= 0.0) {
1701 Int len = (Int)(truncate*Double(convSampling_)+0.5);
1702 for (Int i = 0 ; i < len ; i++) {
1703 val = Double(i) / Double(convSampling_);
1704 grdgauss(&hwhm, &val, &out1);
1705 grdjinc1(&c, &val, &normalize, &out2);
1706 convFunc(i) = out1 * out2;
1707 }
1708 }
1709 else {
1710 Int len = convFunc.nelements();
1711 for (Int i = 0 ; i < len ; i++) {
1712 val = Double(i) / Double(convSampling_);
1713 grdjinc1(&c, &val, &normalize, &out2);
1714 if (out2 <= 0.0) {
1715 LogIO os(LogOrigin("STGrid","gjincFunc",WHERE));
1716 os << LogIO::DEBUG1 << "convFunc is automatically truncated at radius " << val << LogIO::POST;
1717 break;
1718 }
1719 grdgauss(&hwhm, &val, &out1);
1720 convFunc(i) = out1 * out2;
1721 }
1722 }
1723}
1724
1725void STGrid::pbFunc( Vector<Float> &convFunc )
1726{
1727 convFunc = 0.0 ;
1728}
1729
1730vector<float> STGrid::getConvFunc()
1731{
1732 LogIO os(LogOrigin("STGrid","getConvFunc",WHERE));
1733 Vector<Float> convFunc;
1734 vector<float> out;
1735
1736 if (cellx_ <= 0.0 || celly_ <= 0.0) {
1737 selectData();
1738 setupGrid();
1739 }
1740
1741 if (convType_ == "BOX" || convType_ == "SF") {
1742 setConvFunc(convFunc);
1743 }
1744 else if (convType_ == "GAUSS") {
1745 Quantum<Double> q1,q2;
1746 readQuantity(q1,gwidth_);
1747 readQuantity(q2,truncate_);
1748// if (celly_ <= 0.0
1749// && ((!q1.getUnit().empty()&&q1.getUnit()!="pixel") ||
1750// (!q2.getUnit().empty()&&q2.getUnit()!="pixel"))) {
1751// throw AipsError("You have to call defineImage to get correct convFunc");
1752// }
1753 setConvFunc(convFunc);
1754 }
1755 else if (convType_ == "GJINC") {
1756 Quantum<Double> q1,q2,q3;
1757 readQuantity(q1,gwidth_);
1758 readQuantity(q2,truncate_);
1759 readQuantity(q3,jwidth_);
1760// if (celly_ <= 0.0
1761// && ((!q1.getUnit().empty()&&q1.getUnit()!="pixel") ||
1762// (!q2.getUnit().empty()&&q2.getUnit()!="pixel") ||
1763// (!q3.getUnit().empty()&&q3.getUnit()!="pixel"))) {
1764// throw AipsError("You have to call defineImage to get correct convFunc");
1765// }
1766 setConvFunc(convFunc);
1767 }
1768 else if (convType_ == "PB") {
1769 throw AipsError("Grid function PB is not available");
1770 }
1771 else {
1772 throw AipsError("Unknown grid function: "+convType_);
1773 }
1774
1775 convFunc.tovector(out);
1776 return out;
1777}
1778
1779void STGrid::setConvFunc( Vector<Float> &convFunc )
1780{
1781 LogIO os(LogOrigin("STGrid","setConvFunc",WHERE));
1782 convSupport_ = userSupport_ ;
1783 if ( convType_ == "BOX" ) {
1784 if ( convSupport_ < 0 )
1785 convSupport_ = 0 ;
1786 Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
1787 convFunc.resize( convSize ) ;
1788 boxFunc( convFunc, convSize ) ;
1789 os << LogIO::DEBUGGING
1790 << "convType_ = " << convType_ << endl
1791 << "convSupport_ = " << convSupport_ << LogIO::POST;
1792 }
1793 else if ( convType_ == "SF" ) {
1794 if ( convSupport_ < 0 )
1795 convSupport_ = 3 ;
1796 Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
1797 convFunc.resize( convSize ) ;
1798 spheroidalFunc( convFunc ) ;
1799 os << LogIO::DEBUGGING
1800 << "convType_ = " << convType_ << endl
1801 << "convSupport_ = " << convSupport_ << LogIO::POST;
1802 }
1803 else if ( convType_ == "GAUSS" ) {
1804 // determine pixel gwidth
1805 // default is HWHM corresponding to b = 1.0 (Mangum et al. 2007)
1806 Double pixelGW = -1.0;
1807 Quantum<Double> q ;
1808 if (!gwidth_.empty()) {
1809 readQuantity( q, gwidth_ );
1810 if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
1811 pixelGW = q.getValue();
1812 }
1813 else {
1814 pixelGW = q.getValue("rad")/celly_;
1815 }
1816 }
1817 pixelGW = (pixelGW >= 0.0) ? pixelGW : sqrt(log(2.0));
1818 if (pixelGW < 0.0) {
1819 os << LogIO::SEVERE
1820 << "Negative width is specified for gaussian" << LogIO::EXCEPTION;
1821 }
1822 // determine truncation radius
1823 // default is 3 * HWHM
1824 Double truncate = -1.0;
1825 if (!truncate_.empty()) {
1826 readQuantity( q, truncate_ );
1827 if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
1828 truncate = q.getValue();
1829 }
1830 else {
1831 truncate = q.getValue("rad")/celly_;
1832 }
1833 }
1834 //convSupport_ = (Int)(truncate+0.5);
1835 truncate = (truncate >= 0.0) ? truncate : 3.0 * pixelGW;
1836 convSupport_ = Int(truncate);
1837 convSupport_ += (((truncate-(Double)convSupport_) > 0.0) ? 1 : 0);
1838 Int convSize = convSampling_ * ( 2*convSupport_ + 2 ) ;
1839 convFunc.resize( convSize ) ;
1840 gaussFunc( convFunc, pixelGW, truncate ) ;
1841 os << LogIO::DEBUGGING
1842 << "convType_ = " << convType_ << endl
1843 << "convSupport_ = " << convSupport_ << endl
1844 << "truncate_ = " << truncate << "pixel" << endl
1845 << "gwidth_ = " << pixelGW << "pixel" << LogIO::POST;
1846 }
1847 else if ( convType_ == "GJINC" ) {
1848 // determine pixel gwidth
1849 // default is HWHM corresponding to b = 2.52 (Mangum et al. 2007)
1850 Double pixelGW = -1.0;
1851 Quantum<Double> q ;
1852 if (!gwidth_.empty()) {
1853 readQuantity( q, gwidth_ );
1854 if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
1855 pixelGW = q.getValue();
1856 }
1857 else {
1858 pixelGW = q.getValue("rad")/celly_;
1859 }
1860 }
1861 pixelGW = (pixelGW >= 0.0) ? pixelGW : sqrt(log(2.0)) * 2.52;
1862 if (pixelGW < 0.0) {
1863 os << LogIO::SEVERE
1864 << "Negative width is specified for gaussian" << LogIO::EXCEPTION;
1865 }
1866 // determine pixel c
1867 // default is c = 1.55 (Mangum et al. 2007)
1868 Double pixelJW = -1.0;
1869 if (!jwidth_.empty()) {
1870 readQuantity( q, jwidth_ );
1871 if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
1872 pixelJW = q.getValue();
1873 }
1874 else {
1875 pixelJW = q.getValue("rad")/celly_;
1876 }
1877 }
1878 pixelJW = (pixelJW >= 0.0) ? pixelJW : 1.55;
1879 if (pixelJW < 0.0) {
1880 os << LogIO::SEVERE
1881 << "Negative width is specified for jinc" << LogIO::EXCEPTION;
1882 }
1883 // determine truncation radius
1884 // default is -1.0 (truncate at first null)
1885 Double truncate = -1.0;
1886 if (!truncate_.empty()) {
1887 readQuantity( q, truncate_ );
1888 if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
1889 truncate = q.getValue();
1890 }
1891 else {
1892 truncate = q.getValue("rad")/celly_;
1893 }
1894 }
1895 //convSupport_ = (truncate >= 0.0) ? (Int)(truncate+0.5) : (Int)(2*pixelJW+0.5);
1896 Double convSupportF = (truncate >= 0.0) ? truncate : (2*pixelJW);
1897 convSupport_ = (Int)convSupportF;
1898 convSupport_ += (((convSupportF-(Double)convSupport_) > 0.0) ? 1 : 0);
1899 Int convSize = convSampling_ * ( 2*convSupport_ + 2 ) ;
1900 convFunc.resize( convSize ) ;
1901 gjincFunc( convFunc, pixelGW, pixelJW, truncate ) ;
1902 os << LogIO::DEBUGGING
1903 << "convType_ = " << convType_ << endl
1904 << "convSupport_ = " << convSupport_ << endl
1905 << "truncate_ = " << truncate << "pixel" << endl
1906 << "gwidth_ = " << pixelGW << "pixel" << endl
1907 << "jwidth_ = " << pixelJW << "pixel" << LogIO::POST;
1908 }
1909 else if ( convType_ == "PB" ) {
1910 if ( convSupport_ < 0 )
1911 convSupport_ = 0 ;
1912 pbFunc( convFunc ) ;
1913 }
1914 else {
1915 throw AipsError( "Unsupported convolution function" ) ;
1916 }
1917}
1918
1919string STGrid::saveData( string outfile )
1920{
1921 LogIO os( LogOrigin("STGrid", "saveData", WHERE) ) ;
1922 double t0, t1 ;
1923 t0 = mathutil::gettimeofday_sec() ;
1924
1925 //Int polno = 0 ;
1926 String outfile_ ;
1927 if ( outfile.size() == 0 ) {
1928 if ( infileList_[0].lastchar() == '/' ) {
1929 outfile_ = infileList_[0].substr( 0, infileList_[0].size()-1 ) ;
1930 }
1931 else {
1932 outfile_ = infileList_[0] ;
1933 }
1934 outfile_ += ".grid" ;
1935 }
1936 else {
1937 outfile_ = outfile ;
1938 }
1939 Table tab ;
1940 prepareTable( tab, outfile_ ) ;
1941 fillTable( tab ) ;
1942
1943 t1 = mathutil::gettimeofday_sec() ;
1944 os << LogIO::DEBUGGING << "saveData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
1945
1946 return outfile_ ;
1947}
1948
1949void STGrid::prepareTable( Table &tab, String &name )
1950{
1951 Table t( infileList_[0], Table::Old ) ;
1952 t.deepCopy( name, Table::New, False, t.endianFormat(), True ) ;
1953 tab = Table( name, Table::Update ) ;
1954 // 2012/02/13 TN
1955 // explicitly copy subtables since no rows including subtables are
1956 // copied by Table::deepCopy with noRows=True
1957 //TableCopy::copySubTables( tab, t ) ;
1958 const TableRecord &inrec = t.keywordSet();
1959 TableRecord &outrec = tab.rwKeywordSet();
1960 for (uInt i = 0 ; i < inrec.nfields() ; i++) {
1961 if (inrec.type(i) == TpTable) {
1962 String name = inrec.name(i);
1963 Table intable = inrec.asTable(name);
1964 Table outtable = outrec.asTable(name);
1965 TableCopy::copyRows(outtable, intable);
1966 }
1967 }
1968}
1969
1970void STGrid::fillTable( Table &tab )
1971{
1972 //IPosition dshape = data_.shape() ;
1973 Int nrow = nx_ * ny_ * npol_ ;
1974 tab.rwKeywordSet().define( "nPol", npol_ ) ;
1975 tab.addRow( nrow ) ;
1976 Vector<Double> cpix( 2 ) ;
1977 cpix(0) = Double( nx_ - 1 ) * 0.5 ;
1978 cpix(1) = Double( ny_ - 1 ) * 0.5 ;
1979 Vector<Double> dir( 2 ) ;
1980 Vector<Double> pix( 2 );
1981 ArrayColumn<Double> directionCol( tab, "DIRECTION" ) ;
1982 ArrayColumn<Float> spectraCol( tab, "SPECTRA" ) ;
1983 ArrayColumn<uChar> flagtraCol( tab, "FLAGTRA" ) ;
1984 ScalarColumn<uInt> flagRowCol( tab, "FLAGROW" );
1985 ScalarColumn<uInt> polnoCol( tab, "POLNO" ) ;
1986 ScalarColumn<uInt> scannoCol( tab, "SCANNO" ) ;
1987 Int irow = 0 ;
1988 Vector<Float> sp( nchan_ ) ;
1989 Vector<uChar> flag( nchan_ ) ;
1990 Bool bsp, bdata, bflag ;
1991 const Float *data_p = data_.getStorage( bdata ) ;
1992 const uChar *flag_p = flag_.getStorage( bflag ) ;
1993 Float *wsp_p, *sp_p ;
1994 const Float *wdata_p = data_p ;
1995 const uChar *wflag_p = flag_p ;
1996 long step = nx_ * ny_ * npol_ ;
1997 long offset ;
1998 uInt scanno = 0 ;
1999 uChar rflag;
2000 for ( Int iy = 0 ; iy < ny_ ; iy++ ) {
2001 pix(1) = (Double)(iy);
2002 for ( Int ix = 0 ; ix < nx_ ; ix++ ) {
2003 pix(0) = (Double)(ix);
2004 dircoord_->toWorld(dir,pix);
2005 //os << "dir[" << ix << "," << iy << "]=" << dir << LogIO::POST;
2006 for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
2007 offset = ix + nx_ * (iy + ipol * ny_) ;
2008 //os << "offset = " << offset << LogIO::POST ;
2009 sp_p = sp.getStorage( bsp ) ;
2010 wsp_p = sp_p ;
2011 wdata_p = data_p + offset ;
2012 wflag_p = flag_p + offset ;
2013 rflag = ~0 ; //11111111
2014 for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
2015 *wsp_p = *wdata_p ;
2016 wsp_p++ ;
2017 wdata_p += step ;
2018 flag[ichan] = *wflag_p ;
2019 rflag &= flag[ichan] ;
2020 wflag_p += step ;
2021 }
2022 sp.putStorage( sp_p, bsp ) ;
2023 spectraCol.put( irow, sp ) ;
2024 flagtraCol.put( irow, flag ) ;
2025 flagRowCol.put( irow, ((rflag>0) ? 1 : 0) ) ;
2026 directionCol.put( irow, dir ) ;
2027 polnoCol.put( irow, pollist_[ipol] ) ;
2028 scannoCol.put( irow, scanno ) ;
2029 irow++ ;
2030 }
2031 scanno++ ;
2032 }
2033 }
2034 data_.freeStorage( data_p, bdata ) ;
2035 flag_.freeStorage( flag_p, bflag ) ;
2036
2037 fillMainColumns( tab ) ;
2038}
2039
2040void STGrid::fillMainColumns( Table &tab )
2041{
2042 // values for fill
2043 //Table t( infileList_[0], Table::Old ) ;
2044 Table t ;
2045 table( t, 0 ) ;
2046 Table tsel = t( t.col( "IFNO" ) == (uInt)ifno_, 1 ) ;
2047 ROTableRow row( tsel ) ;
2048 row.get( 0 ) ;
2049 const TableRecord &rec = row.record() ;
2050 uInt freqId = rec.asuInt( "FREQ_ID" ) ;
2051 uInt molId = rec.asuInt( "MOLECULE_ID" ) ;
2052 uInt tcalId = rec.asuInt( "TCAL_ID" ) ;
2053 uInt focusId = rec.asuInt( "FOCUS_ID" ) ;
2054 uInt weatherId = rec.asuInt( "WEATHER_ID" ) ;
2055 String srcname = rec.asString( "SRCNAME" ) ;
2056 String fieldname = rec.asString( "FIELDNAME" ) ;
2057 Vector<Float> defaultTsys( 1, 1.0 ) ;
2058 // @todo how to set flagtra for gridded spectra?
2059 Vector<uChar> flagtra = rec.asArrayuChar( "FLAGTRA" ) ;
2060 flagtra = (uChar)0 ;
2061 Float opacity = rec.asFloat( "OPACITY" ) ;
2062 Double srcvel = rec.asDouble( "SRCVELOCITY" ) ;
2063 Vector<Double> srcpm = rec.asArrayDouble( "SRCPROPERMOTION" ) ;
2064 Vector<Double> srcdir = rec.asArrayDouble( "SRCDIRECTION" ) ;
2065 Vector<Double> scanrate = rec.asArrayDouble( "SCANRATE" ) ;
2066 Double time = rec.asDouble( "TIME" ) ;
2067 Double interval = rec.asDouble( "INTERVAL" ) ;
2068
2069 // fill columns
2070 Int nrow = tab.nrow() ;
2071 ScalarColumn<uInt> ifnoCol( tab, "IFNO" ) ;
2072 ScalarColumn<uInt> beamnoCol(tab, "BEAMNO");
2073 ScalarColumn<uInt> freqIdCol( tab, "FREQ_ID" ) ;
2074 ScalarColumn<uInt> molIdCol( tab, "MOLECULE_ID" ) ;
2075 ScalarColumn<uInt> tcalidCol( tab, "TCAL_ID" ) ;
2076 ScalarColumn<Int> fitidCol( tab, "FIT_ID" ) ;
2077 ScalarColumn<uInt> focusidCol( tab, "FOCUS_ID" ) ;
2078 ScalarColumn<uInt> weatheridCol( tab, "WEATHER_ID" ) ;
2079 ArrayColumn<uChar> flagtraCol( tab, "FLAGTRA" ) ;
2080 ScalarColumn<uInt> rflagCol( tab, "FLAGROW" ) ;
2081 ArrayColumn<Float> tsysCol( tab, "TSYS" ) ;
2082 ScalarColumn<String> srcnameCol( tab, "SRCNAME" ) ;
2083 ScalarColumn<String> fieldnameCol( tab, "FIELDNAME" ) ;
2084 ScalarColumn<Int> srctypeCol( tab, "SRCTYPE" ) ;
2085 ScalarColumn<Float> opacityCol( tab, "OPACITY" ) ;
2086 ScalarColumn<Double> srcvelCol( tab, "SRCVELOCITY" ) ;
2087 ArrayColumn<Double> srcpmCol( tab, "SRCPROPERMOTION" ) ;
2088 ArrayColumn<Double> srcdirCol( tab, "SRCDIRECTION" ) ;
2089 ArrayColumn<Double> scanrateCol( tab, "SCANRATE" ) ;
2090 ScalarColumn<Double> timeCol( tab, "TIME" ) ;
2091 ScalarColumn<Double> intervalCol( tab, "INTERVAL" ) ;
2092 for ( Int i = 0 ; i < nrow ; i++ ) {
2093 ifnoCol.put( i, (uInt)ifno_ ) ;
2094 beamnoCol.put(i, 0);
2095 freqIdCol.put( i, freqId ) ;
2096 molIdCol.put( i, molId ) ;
2097 tcalidCol.put( i, tcalId ) ;
2098 fitidCol.put( i, -1 ) ;
2099 focusidCol.put( i, focusId ) ;
2100 weatheridCol.put( i, weatherId ) ;
2101 //flagtraCol.put( i, flagtra ) ;
2102 //rflagCol.put( i, 0 ) ;
2103 tsysCol.put( i, defaultTsys ) ;
2104 srcnameCol.put( i, srcname ) ;
2105 fieldnameCol.put( i, fieldname ) ;
2106 srctypeCol.put( i, (Int)SrcType::PSON ) ;
2107 opacityCol.put( i, opacity ) ;
2108 srcvelCol.put( i, srcvel ) ;
2109 srcpmCol.put( i, srcpm ) ;
2110 srcdirCol.put( i, srcdir ) ;
2111 scanrateCol.put( i, scanrate ) ;
2112 timeCol.put( i, time ) ;
2113 intervalCol.put( i, interval ) ;
2114 if ( (i + 1) % npol_ == 0 ) {
2115 time += interval / 86400.0;
2116 }
2117 }
2118}
2119
2120vector<int> STGrid::getResultantMapSize()
2121{
2122 vector<int> r(2);
2123 r[0] = nx_;
2124 r[1] = ny_;
2125 return r;
2126}
2127
2128vector<double> STGrid::getResultantCellSize()
2129{
2130 vector<double> r(2);
2131 r[0] = cellx_;
2132 r[1] = celly_;
2133 return r;
2134}
2135
2136// STGrid2
2137STGrid2::STGrid2()
2138 : STGrid()
2139{
2140}
2141
2142STGrid2::STGrid2( const ScantableWrapper &s )
2143 : STGrid()
2144{
2145 setScantable( s ) ;
2146}
2147
2148STGrid2::STGrid2( const vector<ScantableWrapper> &v )
2149 : STGrid()
2150{
2151 setScantableList( v ) ;
2152}
2153
2154void STGrid2::setScantable( const ScantableWrapper &s )
2155{
2156 nfile_ = 1 ;
2157 dataList_.resize( nfile_ ) ;
2158 dataList_[0] = s ;
2159 infileList_.resize( nfile_ ) ;
2160 infileList_[0] = s.getCP()->table().tableName() ;
2161}
2162
2163void STGrid2::setScantableList( const vector<ScantableWrapper> &v )
2164{
2165 nfile_ = v.size() ;
2166 dataList_.resize( nfile_ ) ;
2167 infileList_.resize( nfile_ ) ;
2168 for ( uInt i = 0 ; i < nfile_ ; i++ ) {
2169 dataList_[i] = v[i] ;
2170 infileList_[i] = v[i].getCP()->table().tableName() ;
2171 }
2172}
2173
2174ScantableWrapper STGrid2::getResultAsScantable( int tp )
2175{
2176 ScantableWrapper sw( tp ) ;
2177 CountedPtr<Scantable> s = sw.getCP() ;
2178 s->setHeader( dataList_[0].getCP()->getHeader() ) ;
2179 Table tout, tin ;
2180 String subt[] = { "FREQUENCIES", "FOCUS", "WEATHER",
2181 "TCAL", "MOLECULES", "HISTORY", "FIT" } ;
2182 for ( uInt i = 0 ; i < 7 ; i++ ) {
2183 tout = s->table().rwKeywordSet().asTable(subt[i]) ;
2184 tin = dataList_[0].getCP()->table().rwKeywordSet().asTable(subt[i]) ;
2185 TableCopy::copyRows( tout, tin ) ;
2186 tout.rwKeywordSet() = tin.rwKeywordSet();
2187 }
2188 fillTable( s->table() ) ;
2189 return sw ;
2190}
2191
2192void STGrid2::table( Table &tab, uInt i )
2193{
2194 if ( i < nfile_ )
2195 tab = dataList_[i].getCP()->table() ;
2196}
2197
2198}
Note: See TracBrowser for help on using the repository browser.