source: trunk/src/STGrid.cpp@ 2397

Last change on this file since 2397 was 2397, checked in by Takeshi Nakazato, 14 years ago

New Development: No

JIRA Issue: Yes CAS-2816

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

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