source: trunk/src/STGrid.cpp@ 2390

Last change on this file since 2390 was 2390, checked in by Takeshi Nakazato, 13 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...

More than one input data are supported.
File size: 36.7 KB
Line 
1#include <iostream>
2#include <fstream>
3
4#include <casa/BasicSL/String.h>
5#include <casa/Arrays/Vector.h>
6#include <casa/Arrays/Matrix.h>
7#include <casa/Arrays/Cube.h>
8#include <casa/Arrays/ArrayMath.h>
9#include <casa/Arrays/ArrayIter.h>
10#include <casa/Quanta/Quantum.h>
11#include <casa/Quanta/QuantumHolder.h>
12#include <casa/Utilities/CountedPtr.h>
13#include <casa/Logging/LogIO.h>
14
15#include <tables/Tables/Table.h>
16#include <tables/Tables/TableRecord.h>
17#include <tables/Tables/ExprNode.h>
18#include <tables/Tables/ScalarColumn.h>
19#include <tables/Tables/ArrayColumn.h>
20#include <tables/Tables/TableIter.h>
21
22#include <measures/Measures/MDirection.h>
23
24#include <MathUtils.h>
25
26#include "STGrid.h"
27
28using namespace std ;
29using namespace casa ;
30using namespace asap ;
31
32namespace asap {
33
34// for performance check
35double eToInt = 0.0 ;
36double eGetWeight = 0.0 ;
37
38// constructor
39STGrid::STGrid()
40{
41 init() ;
42}
43
44STGrid::STGrid( const string infile )
45{
46 init() ;
47
48 setFileIn( infile ) ;
49}
50
51STGrid::STGrid( const vector<string> infile )
52{
53 init() ;
54
55 setFileList( infile ) ;
56}
57
58void STGrid::init()
59{
60 ifno_ = -1 ;
61 nx_ = -1 ;
62 ny_ = -1 ;
63 npol_ = 0 ;
64 nchan_ = 0 ;
65 nrow_ = 0 ;
66 cellx_ = 0.0 ;
67 celly_ = 0.0 ;
68 center_ = Vector<Double> ( 2, 0.0 ) ;
69 convType_ = "BOX" ;
70 wtype_ = "UNIFORM" ;
71 convSupport_ = -1 ;
72 userSupport_ = -1 ;
73 convSampling_ = 100 ;
74 nprocessed_ = 0 ;
75 nchunk_ = 0 ;
76
77 // initialize user input
78 nxUI_ = -1 ;
79 nyUI_ = -1 ;
80 cellxUI_ = "" ;
81 cellyUI_ = "" ;
82 centerUI_ = "" ;
83}
84
85void STGrid::setFileIn( const string infile )
86{
87 String name( infile ) ;
88 if ( infileList_.size() == 0 || infileList_[0].compare( name ) != 0 ) {
89 infileList_.resize( 1 ) ;
90 infileList_[0] = String(infile) ;
91 nfile_ = 1 ;
92 }
93}
94
95void STGrid::setFileList( const vector<string> infile )
96{
97 nfile_ = infile.size() ;
98 infileList_.resize( nfile_ ) ;
99 for ( uInt i = 0 ; i < nfile_ ; i++ ) {
100 infileList_[i] = infile[i] ;
101 }
102}
103
104void STGrid::setPolList( vector<unsigned int> pols )
105{
106 pollist_.assign( Vector<uInt>( pols ) ) ;
107}
108
109void STGrid::setScanList( vector<unsigned int> scans )
110{
111 scanlist_.assign( Vector<uInt>( scans ) ) ;
112}
113
114void STGrid::setWeight( const string wType )
115{
116 wtype_ = String( wType ) ;
117 wtype_.upcase() ;
118}
119
120void STGrid::defineImage( int nx,
121 int ny,
122 string scellx,
123 string scelly,
124 string scenter )
125{
126 nxUI_ = (Int)nx ;
127 nyUI_ = (Int)ny ;
128 cellxUI_ = String( scellx ) ;
129 cellyUI_ = String( scelly ) ;
130 centerUI_ = String( scenter ) ;
131}
132
133void STGrid::setFunc( string convType,
134 int convSupport )
135{
136 convType_ = String( convType ) ;
137 convType_.upcase() ;
138 userSupport_ = (Int)convSupport ;
139}
140
141#define NEED_UNDERSCORES
142#if defined(NEED_UNDERSCORES)
143#define ggridsd ggridsd_
144#endif
145extern "C" {
146 void ggridsd(Double*,
147 const Complex*,
148 Int*,
149 Int*,
150 Int*,
151 const Int*,
152 const Int*,
153 const Float*,
154 Int*,
155 Int*,
156 Complex*,
157 Float*,
158 Int*,
159 Int*,
160 Int *,
161 Int *,
162 Int*,
163 Int*,
164 Float*,
165 Int*,
166 Int*,
167 Double*);
168}
169void STGrid::call_ggridsd( Array<Double> &xypos,
170 Array<Complex> &spectra,
171 Int &nvispol,
172 Int &nvischan,
173 Array<Int> &flagtra,
174 Array<Int> &flagrow,
175 Array<Float> &weight,
176 Int &nrow,
177 Int &irow,
178 Array<Complex> &gdata,
179 Array<Float> &gwgt,
180 Int &nx,
181 Int &ny,
182 Int &npol,
183 Int &nchan,
184 Int &support,
185 Int &sampling,
186 Vector<Float> &convFunc,
187 Int *chanMap,
188 Int *polMap )
189{
190 // parameters for gridding
191 Int idopsf = 0 ;
192 Int len = npol*nchan ;
193 Double *sumw_p = new Double[len] ;
194 {
195 Double *work_p = sumw_p ;
196 for ( Int i = 0 ; i < len ; i++ ) {
197 *work_p = 0.0 ;
198 work_p++ ;
199 }
200 }
201
202 // prepare pointer
203 Bool deletePos, deleteData, deleteWgt, deleteFlag, deleteFlagR, deleteConv, deleteDataG, deleteWgtG ;
204 Double *xy_p = xypos.getStorage( deletePos ) ;
205 const Complex *values_p = spectra.getStorage( deleteData ) ;
206 const Int *flag_p = flagtra.getStorage( deleteFlag ) ;
207 const Int *rflag_p = flagrow.getStorage( deleteFlagR ) ;
208 const Float *wgt_p = weight.getStorage( deleteWgt ) ;
209 Complex *grid_p = gdata.getStorage( deleteDataG ) ;
210 Float *wgrid_p = gwgt.getStorage( deleteWgtG ) ;
211 Float *conv_p = convFunc.getStorage( deleteConv ) ;
212
213 // pass copy of irow to ggridsd since it will be modified in theroutine
214 Int irowCopy = irow ;
215
216 // call ggridsd
217 ggridsd( xy_p,
218 values_p,
219 &nvispol,
220 &nvischan,
221 &idopsf,
222 flag_p,
223 rflag_p,
224 wgt_p,
225 &nrow,
226 &irowCopy,
227 grid_p,
228 wgrid_p,
229 &nx,
230 &ny,
231 &npol,
232 &nchan,
233 &support,
234 &sampling,
235 conv_p,
236 chanMap,
237 polMap,
238 sumw_p ) ;
239
240 // finalization
241 xypos.putStorage( xy_p, deletePos ) ;
242 spectra.freeStorage( values_p, deleteData ) ;
243 flagtra.freeStorage( flag_p, deleteFlag ) ;
244 flagrow.freeStorage( rflag_p, deleteFlagR ) ;
245 weight.freeStorage( wgt_p, deleteWgt ) ;
246 gdata.putStorage( grid_p, deleteDataG ) ;
247 gwgt.putStorage( wgrid_p, deleteWgtG ) ;
248 convFunc.putStorage( conv_p, deleteConv ) ;
249 delete sumw_p ;
250}
251
252Bool STGrid::pastEnd()
253{
254 LogIO os( LogOrigin("STGrid","pastEnd",WHERE) ) ;
255 Bool b = nprocessed_ >= nrow_ ;
256 return b ;
257}
258
259void STGrid::grid()
260{
261 LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
262 double t0,t1 ;
263
264 // data selection
265 t0 = mathutil::gettimeofday_sec() ;
266 selectData() ;
267 t1 = mathutil::gettimeofday_sec() ;
268 os << "selectData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
269
270 setupGrid() ;
271 setupArray() ;
272
273 if ( wtype_.compare("UNIFORM") != 0 &&
274 wtype_.compare("TINT") != 0 &&
275 wtype_.compare("TSYS") != 0 &&
276 wtype_.compare("TINTSYS") != 0 ) {
277 LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
278 os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
279 wtype_ = "UNIFORM" ;
280 }
281
282 // grid parameter
283 os << LogIO::DEBUGGING ;
284 os << "----------" << endl ;
285 os << "Data selection summary" << endl ;
286 os << " ifno = " << ifno_ << endl ;
287 os << " pollist = " << pollist_ << endl ;
288 os << " scanlist = " << scanlist_ << endl ;
289 os << "----------" << endl ;
290 os << "Grid parameter summary" << endl ;
291 os << " (nx,ny) = (" << nx_ << "," << ny_ << ")" << endl ;
292 os << " (cellx,celly) = (" << cellx_ << "," << celly_ << ")" << endl ;
293 os << " center = " << center_ << endl ;
294 os << " weighting = " << wtype_ << endl ;
295 os << " convfunc = " << convType_ << " with support " << convSupport_ << endl ;
296 os << "----------" << LogIO::POST ;
297 os << LogIO::NORMAL ;
298
299 Bool doAll = examine() ;
300
301 if ( doAll )
302 gridPerPol() ;
303 else
304 gridPerRow() ;
305}
306
307Bool STGrid::examine()
308{
309 // TODO: nchunk_ must be determined from nchan_, npol_, and (nx_,ny_)
310 // by considering data size to be allocated for ggridsd input/output
311 nchunk_ = 400 ;
312 Bool b = nchunk_ >= nrow_ ;
313 nchunk_ = min( nchunk_, nrow_ ) ;
314 return b ;
315}
316
317void STGrid::gridPerRow()
318{
319 LogIO os( LogOrigin("STGrid", "gridPerRow", WHERE) ) ;
320 double t0, t1 ;
321
322 // convolution kernel
323 Vector<Float> convFunc ;
324 t0 = mathutil::gettimeofday_sec() ;
325 setConvFunc( convFunc ) ;
326 t1 = mathutil::gettimeofday_sec() ;
327 os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
328
329 // grid data
330 Int gnx = nx_ ;
331 Int gny = ny_ ;
332 // Extend grid plane with convSupport_
333// Int gnx = nx_+convSupport_*2 ;
334// Int gny = ny_+convSupport_*2 ;
335 IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
336 Array<Complex> gdataArrC( gshape, 0.0 ) ;
337 // 2011/12/20 TN
338 // data_ and gwgtArr share storage
339 data_.resize( gshape ) ;
340 data_ = 0.0 ;
341 Array<Float> gwgtArr( data_ ) ;
342
343 // parameters for gridding
344 Int *chanMap = new Int[nchan_] ;
345 {
346 Int *work_p = chanMap ;
347 for ( Int i = 0 ; i < nchan_ ; i++ ) {
348 *work_p = i ;
349 work_p++ ;
350 }
351 }
352 Int *polMap = new Int[1] ;
353 Int nvispol = 1 ;
354 Int irow = -1 ;
355
356 // for performance check
357 double eGetData = 0.0 ;
358 double eToPixel = 0.0 ;
359 double eGGridSD = 0.0 ;
360 double eInitPol = 0.0 ;
361
362 for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) {
363 initTable( ifile ) ;
364
365 os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ;
366
367 for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
368 t0 = mathutil::gettimeofday_sec() ;
369 initPol( ipol ) ;
370 t1 = mathutil::gettimeofday_sec() ;
371 eInitPol += t1-t0 ;
372
373 polMap[0] = ipol ;
374
375 os << "start pol " << ipol << LogIO::POST ;
376
377 while( !pastEnd() ) {
378 // data storage
379 IPosition cshape( 3, npol_, nchan_, nchunk_ ) ;
380 IPosition mshape( 2, npol_, nchunk_ ) ;
381 IPosition vshape( 1, nchunk_ ) ;
382 IPosition dshape( 2, 2, nchunk_ ) ;
383 IPosition wshape( 2, nchan_, nchunk_ ) ;
384 Array<Complex> spectra( wshape ) ;
385 Array<Int> flagtra( wshape ) ;
386 Array<Int> rflag( vshape ) ;
387 Array<Float> weight( wshape ) ;
388 Array<Double> direction( dshape ) ;
389 Array<Double> xypos( dshape ) ;
390
391 spectraF_.resize( wshape ) ;
392 flagtraUC_.resize( wshape ) ;
393 rflagUI_.resize( vshape ) ;
394
395 // retrieve data
396 t0 = mathutil::gettimeofday_sec() ;
397 Int nrow = getDataChunk( spectra, direction, flagtra, rflag, weight ) ;
398 t1 = mathutil::gettimeofday_sec() ;
399 eGetData += t1-t0 ;
400
401 // world -> pixel
402 t0 = mathutil::gettimeofday_sec() ;
403 toPixel( direction, xypos ) ;
404 t1 = mathutil::gettimeofday_sec() ;
405 eToPixel += t1-t0 ;
406
407 // call ggridsd
408 t0 = mathutil::gettimeofday_sec() ;
409 call_ggridsd( xypos,
410 spectra,
411 nvispol,
412 nchan_,
413 flagtra,
414 rflag,
415 weight,
416 nrow,
417 irow,
418 gdataArrC,
419 gwgtArr,
420 gnx,
421 gny,
422 npol_,
423 nchan_,
424 convSupport_,
425 convSampling_,
426 convFunc,
427 chanMap,
428 polMap ) ;
429 t1 = mathutil::gettimeofday_sec() ;
430 eGGridSD += t1-t0 ;
431
432 }
433
434 os << "end pol " << ipol << LogIO::POST ;
435
436 nprocessed_ = 0 ;
437 }
438 }
439 os << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ;
440 os << "getData: elapsed time is " << eGetData-eToInt-eGetWeight << " sec." << LogIO::POST ;
441 os << "toPixel: elapsed time is " << eToPixel << " sec." << LogIO::POST ;
442 os << "ggridsd: elapsed time is " << eGGridSD << " sec." << LogIO::POST ;
443 os << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
444 os << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
445
446 delete polMap ;
447 delete chanMap ;
448
449 // set data
450 setData( gdataArrC, gwgtArr ) ;
451
452}
453
454void STGrid::gridPerPol()
455{
456 LogIO os( LogOrigin("STGrid", "gridPerPol", WHERE) ) ;
457 double t0, t1 ;
458
459 // convolution kernel
460 Vector<Float> convFunc ;
461 t0 = mathutil::gettimeofday_sec() ;
462 setConvFunc( convFunc ) ;
463 t1 = mathutil::gettimeofday_sec() ;
464 os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
465
466 // prepare grid data storage
467 Int gnx = nx_ ;
468 Int gny = ny_ ;
469// // Extend grid plane with convSupport_
470// Int gnx = nx_+convSupport_*2 ;
471// Int gny = ny_+convSupport_*2 ;
472 IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
473 Array<Complex> gdataArrC( gshape, 0.0 ) ;
474 //Array<Float> gwgtArr( gshape, 0.0 ) ;
475 // 2011/12/20 TN
476 // data_ and weight array shares storage
477 data_.resize( gshape ) ;
478 data_ = 0.0 ;
479 Array<Float> gwgtArr( data_ ) ;
480
481 // maps
482 Int *chanMap = new Int[nchan_] ;
483 {
484 Int *work_p = chanMap ;
485 for ( Int i = 0 ; i < nchan_ ; i++ ) {
486 *work_p = i ;
487 work_p++ ;
488 }
489 }
490 Int *polMap = new Int[1] ;
491
492 // some parameters for ggridsd
493 Int nvispol = 1 ;
494 Int irow = -1 ;
495
496 // for performance check
497 double eInitPol = 0.0 ;
498 double eGetData = 0.0 ;
499 double eToPixel = 0.0 ;
500 double eGGridSD = 0.0 ;
501 double eToInt = 0.0 ;
502
503 for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) {
504 initTable( ifile ) ;
505
506 os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ;
507 // to read data from the table
508 IPosition mshape( 2, nchan_, nrow_ ) ;
509 IPosition dshape( 2, 2, nrow_ ) ;
510 Array<Complex> spectra( mshape ) ;
511 Array<Double> direction( dshape ) ;
512 Array<Int> flagtra( mshape ) ;
513 Array<Int> rflag( IPosition(1,nrow_) ) ;
514 Array<Float> weight( mshape ) ;
515 Array<Double> xypos( dshape ) ;
516
517 for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
518 t0 = mathutil::gettimeofday_sec() ;
519 initPol( ipol ) ;
520 t1 = mathutil::gettimeofday_sec() ;
521 eInitPol += t1-t0 ;
522
523 // retrieve data
524 t0 = mathutil::gettimeofday_sec() ;
525 getDataPerPol( spectra, direction, flagtra, rflag, weight ) ;
526 t1 = mathutil::gettimeofday_sec() ;
527 eGetData += t1-t0 ;
528
529 // world -> pixel
530 t0 = mathutil::gettimeofday_sec() ;
531 toPixel( direction, xypos ) ;
532 t1 = mathutil::gettimeofday_sec() ;
533 eToPixel += t1-t0 ;
534
535 // call ggridsd
536 polMap[0] = ipol ;
537 t0 = mathutil::gettimeofday_sec() ;
538 call_ggridsd( xypos,
539 spectra,
540 nvispol,
541 nchan_,
542 flagtra,
543 rflag,
544 weight,
545 nrow_,
546 irow,
547 gdataArrC,
548 gwgtArr,
549 gnx,
550 gny,
551 npol_,
552 nchan_,
553 convSupport_,
554 convSampling_,
555 convFunc,
556 chanMap,
557 polMap ) ;
558 t1 = mathutil::gettimeofday_sec() ;
559 eGGridSD += t1-t0 ;
560 }
561 }
562 os << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ;
563 os << "getData: elapsed time is " << eGetData-eToInt-eGetWeight << " sec." << LogIO::POST ;
564 os << "toPixel: elapsed time is " << eToPixel << " sec." << LogIO::POST ;
565 os << "ggridsd: elapsed time is " << eGGridSD << " sec." << LogIO::POST ;
566 os << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
567 os << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
568
569 // delete maps
570 delete polMap ;
571 delete chanMap ;
572
573 setData( gdataArrC, gwgtArr ) ;
574// os << "gdataArr = " << gdataArr << LogIO::POST ;
575// os << "gwgtArr = " << gwgtArr << LogIO::POST ;
576// os << "data_ " << data_ << LogIO::POST ;
577}
578
579void STGrid::initPol( Int ipol )
580{
581 LogIO os( LogOrigin("STGrid","initPol",WHERE) ) ;
582 if ( npolOrg_ == 1 ) {
583 os << "single polarization data." << LogIO::POST ;
584 ptab_ = tab_ ;
585 }
586 else
587 ptab_ = tab_( tab_.col("POLNO") == (uInt)ipol ) ;
588
589 attach( ptab_ ) ;
590}
591
592void STGrid::initTable( uInt idx )
593{
594 tab_ = tableList_[idx] ;
595 nrow_ = rows_[idx] ;
596}
597
598void STGrid::setData( Array<Complex> &gdata,
599 Array<Float> &gwgt )
600{
601 // 2011/12/20 TN
602 // gwgt and data_ share storage
603 LogIO os( LogOrigin("STGrid","setData",WHERE) ) ;
604 double t0, t1 ;
605 t0 = mathutil::gettimeofday_sec() ;
606 uInt len = data_.nelements() ;
607 const Complex *w1_p ;
608 Float *w2_p ;
609 Bool b1, b2 ;
610 const Complex *gdata_p = gdata.getStorage( b1 ) ;
611 Float *gwgt_p = data_.getStorage( b2 ) ;
612 w1_p = gdata_p ;
613 w2_p = gwgt_p ;
614 for ( uInt i = 0 ; i < len ; i++ ) {
615 if ( *w2_p > 0.0 ) *w2_p = (*w1_p).real() / *w2_p ;
616 w1_p++ ;
617 w2_p++ ;
618 }
619 gdata.freeStorage( gdata_p, b1 ) ;
620 data_.putStorage( gwgt_p, b2 ) ;
621 t1 = mathutil::gettimeofday_sec() ;
622 os << "setData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
623}
624
625void STGrid::setupGrid()
626{
627 Double xmin,xmax,ymin,ymax ;
628 mapExtent( xmin, xmax, ymin, ymax ) ;
629
630 setupGrid( nxUI_, nyUI_, cellxUI_, cellyUI_,
631 xmin, xmax, ymin, ymax, centerUI_ ) ;
632}
633
634void STGrid::setupGrid( Int &nx,
635 Int &ny,
636 String &cellx,
637 String &celly,
638 Double &xmin,
639 Double &xmax,
640 Double &ymin,
641 Double &ymax,
642 String &center )
643{
644 LogIO os( LogOrigin("STGrid","setupGrid",WHERE) ) ;
645 //cout << "nx=" << nx << ", ny=" << ny << endl ;
646
647 // center position
648 if ( center.size() == 0 ) {
649 center_(0) = 0.5 * ( xmin + xmax ) ;
650 center_(1) = 0.5 * ( ymin + ymax ) ;
651 }
652 else {
653 String::size_type pos0 = center.find( " " ) ;
654 if ( pos0 == String::npos ) {
655 throw AipsError( "bad string format in parameter center" ) ;
656 }
657 String::size_type pos1 = center.find( " ", pos0+1 ) ;
658 String typestr, xstr, ystr ;
659 if ( pos1 != String::npos ) {
660 typestr = center.substr( 0, pos0 ) ;
661 xstr = center.substr( pos0+1, pos1-pos0 ) ;
662 ystr = center.substr( pos1+1 ) ;
663 // todo: convert to J2000 (or direction ref for DIRECTION column)
664 }
665 else {
666 typestr = "J2000" ;
667 xstr = center.substr( 0, pos0 ) ;
668 ystr = center.substr( pos0+1 ) ;
669 }
670 QuantumHolder qh ;
671 String err ;
672 qh.fromString( err, xstr ) ;
673 Quantum<Double> xcen = qh.asQuantumDouble() ;
674 qh.fromString( err, ystr ) ;
675 Quantum<Double> ycen = qh.asQuantumDouble() ;
676 center_(0) = xcen.getValue( "rad" ) ;
677 center_(1) = ycen.getValue( "rad" ) ;
678 }
679
680
681 nx_ = nx ;
682 ny_ = ny ;
683 if ( nx < 0 && ny > 0 ) {
684 nx_ = ny ;
685 ny_ = ny ;
686 }
687 if ( ny < 0 && nx > 0 ) {
688 nx_ = nx ;
689 ny_ = nx ;
690 }
691
692 //Double wx = xmax - xmin ;
693 //Double wy = ymax - ymin ;
694 Double wx = max( abs(xmax-center_(0)), abs(xmin-center_(0)) ) * 2 ;
695 Double wy = max( abs(ymax-center_(1)), abs(ymin-center_(1)) ) * 2 ;
696 // take 10% margin
697 wx *= 1.10 ;
698 wy *= 1.10 ;
699
700 Quantum<Double> qcellx ;
701 Quantum<Double> qcelly ;
702 //cout << "nx_ = " << nx_ << ", ny_ = " << ny_ << endl ;
703 if ( cellx.size() != 0 && celly.size() != 0 ) {
704 readQuantity( qcellx, cellx ) ;
705 readQuantity( qcelly, celly ) ;
706 }
707 else if ( celly.size() != 0 ) {
708 os << "Using celly to x-axis..." << LogIO::POST ;
709 readQuantity( qcelly, celly ) ;
710 qcellx = qcelly ;
711 }
712 else if ( cellx.size() != 0 ) {
713 os << "Using cellx to y-axis..." << LogIO::POST ;
714 readQuantity( qcellx, cellx ) ;
715 qcelly = qcellx ;
716 }
717 else {
718 if ( nx_ < 0 ) {
719 os << "No user preference in grid setting. Using default..." << LogIO::POST ;
720 readQuantity( qcellx, "1.0arcmin" ) ;
721 qcelly = qcellx ;
722 }
723 else {
724 if ( wx == 0.0 ) {
725 os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
726 wx = 0.00290888 ;
727 }
728 if ( wy == 0.0 ) {
729 os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
730 wy = 0.00290888 ;
731 }
732 qcellx = Quantum<Double>( wx/nx_, "rad" ) ;
733 qcelly = Quantum<Double>( wy/ny_, "rad" ) ;
734 }
735 }
736 cellx_ = qcellx.getValue( "rad" ) ;
737 celly_ = qcelly.getValue( "rad" ) ;
738 if ( nx_ < 0 ) {
739 if ( wx == 0.0 ) {
740 os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
741 wx = 0.00290888 ;
742 }
743 if ( wy == 0.0 ) {
744 os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
745 wy = 0.00290888 ;
746 }
747 nx_ = Int( ceil( wx/cellx_ ) ) ;
748 ny_ = Int( ceil( wy/celly_ ) ) ;
749 }
750}
751
752void STGrid::mapExtent( Double &xmin, Double &xmax,
753 Double &ymin, Double &ymax )
754{
755 //LogIO os( LogOrigin("STGrid","mapExtent",WHERE) ) ;
756 directionCol_.attach( tableList_[0], "DIRECTION" ) ;
757 Matrix<Double> direction = directionCol_.getColumn() ;
758 //os << "dirCol.nrow() = " << dirCol.nrow() << LogIO::POST ;
759 minMax( xmin, xmax, direction.row( 0 ) ) ;
760 minMax( ymin, ymax, direction.row( 1 ) ) ;
761 Double amin, amax, bmin, bmax ;
762 for ( uInt i = 1 ; i < nfile_ ; i++ ) {
763 directionCol_.attach( tableList_[i], "DIRECTION" ) ;
764 direction.assign( directionCol_.getColumn() ) ;
765 //os << "dirCol.nrow() = " << dirCol.nrow() << LogIO::POST ;
766 minMax( amin, amax, direction.row( 0 ) ) ;
767 minMax( bmin, bmax, direction.row( 1 ) ) ;
768 xmin = min( xmin, amin ) ;
769 xmax = max( xmax, amax ) ;
770 ymin = min( ymin, bmin ) ;
771 ymax = max( ymax, bmax ) ;
772 }
773 //os << "(xmin,xmax)=(" << xmin << "," << xmax << ")" << LogIO::POST ;
774 //os << "(ymin,ymax)=(" << ymin << "," << ymax << ")" << LogIO::POST ;
775}
776
777void STGrid::selectData()
778{
779 LogIO os( LogOrigin("STGrid","selectData",WHERE) ) ;
780 Int ifno = ifno_ ;
781 tableList_.resize( nfile_ ) ;
782 for ( uInt i = 0 ; i < nfile_ ; i++ ) {
783 Table taborg( infileList_[i] ) ;
784 if ( ifno == -1 ) {
785 ROScalarColumn<uInt> ifnoCol( taborg, "IFNO" ) ;
786 ifno = ifnoCol( 0 ) ;
787 os << LogIO::WARN
788 << "IFNO is not given. Using default IFNO: " << ifno << LogIO::POST ;
789 }
790 TableExprNode node ;
791 if ( isMultiIF( taborg ) ) {
792 os << "apply selection on IFNO" << LogIO::POST ;
793 node = taborg.col("IFNO") == ifno ;
794 }
795 if ( scanlist_.size() > 0 ) {
796 os << "apply selection on SCANNO" << LogIO::POST ;
797 node = node && taborg.col("SCANNO").in( scanlist_ ) ;
798 }
799 if ( node.isNull() ) {
800 tableList_[i] = taborg ;
801 }
802 else {
803 tableList_[i] = taborg( node ) ;
804 }
805 if ( tableList_[i].nrow() == 0 ) {
806 os << LogIO::SEVERE
807 << "No corresponding rows for given selection: IFNO " << ifno ;
808 if ( scanlist_.size() > 0 )
809 os << " SCANNO " << scanlist_ ;
810 os << LogIO::EXCEPTION ;
811 }
812 }
813}
814
815Bool STGrid::isMultiIF( Table &tab )
816{
817 ROScalarColumn<uInt> ifnoCol( tab, "IFNO" ) ;
818 Vector<uInt> ifnos = ifnoCol.getColumn() ;
819 return anyNE( ifnos, ifnos[0] ) ;
820}
821
822void STGrid::attach( Table &tab )
823{
824 // attach to table
825 spectraCol_.attach( tab, "SPECTRA" ) ;
826 flagtraCol_.attach( tab, "FLAGTRA" ) ;
827 directionCol_.attach( tab, "DIRECTION" ) ;
828 flagRowCol_.attach( tab, "FLAGROW" ) ;
829 tsysCol_.attach( tab, "TSYS" ) ;
830 intervalCol_.attach( tab, "INTERVAL" ) ;
831}
832
833void STGrid::getDataPerPol( Array<Float> &spectra,
834 Array<Double> &direction,
835 Array<uChar> &flagtra,
836 Array<uInt> &rflag,
837 Array<Float> &weight )
838{
839 LogIO os( LogOrigin("STGrid","getDataPerPol",WHERE) ) ;
840
841 // 2011/12/22 TN
842 // weight and tsys shares its storage
843 Array<Float> tsys( weight ) ;
844 Array<Double> tint( rflag.shape() ) ;
845
846 Vector<uInt> rflagVec( rflag ) ;
847 Vector<Double> tintVec( tint ) ;
848
849 spectraCol_.getColumn( spectra ) ;
850 flagtraCol_.getColumn( flagtra ) ;
851 flagRowCol_.getColumn( rflagVec ) ;
852 directionCol_.getColumn( direction ) ;
853 intervalCol_.getColumn( tintVec ) ;
854 Vector<Float> tsysTemp = tsysCol_( 0 ) ;
855 if ( tsysTemp.nelements() == (uInt)nchan_ ) {
856 tsysCol_.getColumn( tsys ) ;
857 }
858 else {
859 tsys = tsysTemp[0] ;
860 }
861
862 double t0,t1 ;
863 t0 = mathutil::gettimeofday_sec() ;
864 getWeight( weight, tsys, tint ) ;
865 t1 = mathutil::gettimeofday_sec() ;
866 eGetWeight += t1-t0 ;
867}
868
869void STGrid::getDataPerPol( Array<Complex> &spectra,
870 Array<Double> &direction,
871 Array<Int> &flagtra,
872 Array<Int> &rflag,
873 Array<Float> &weight )
874{
875 LogIO os( LogOrigin("STGrid","getDataPerPol",WHERE) ) ;
876 double t0, t1 ;
877
878 Array<uChar> flagUC( flagtra.shape() ) ;
879 Array<uInt> rflagUI( rflag.shape() ) ;
880 Array<Float> spectraF( spectra.shape() ) ;
881 getDataPerPol( spectraF, direction, flagUC, rflagUI, weight ) ;
882
883 convertArray( spectra, spectraF ) ;
884 t0 = mathutil::gettimeofday_sec() ;
885 toInt( flagUC, flagtra ) ;
886 toInt( rflagUI, rflag ) ;
887 t1 = mathutil::gettimeofday_sec() ;
888 eToInt += t1-t0 ;
889 //os << "toInt: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
890}
891
892Int STGrid::getDataChunk( Array<Complex> &spectra,
893 Array<Double> &direction,
894 Array<Int> &flagtra,
895 Array<Int> &rflag,
896 Array<Float> &weight )
897{
898 LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
899 Int nrow = getDataChunk( spectraF_, direction, flagtraUC_, rflagUI_, weight ) ;
900 if ( nrow < nchunk_ ) {
901 spectra.resize( spectraF_.shape() ) ;
902 flagtra.resize( flagtraUC_.shape() ) ;
903 rflag.resize( rflagUI_.shape() ) ;
904 }
905 double t0, t1 ;
906 t0 = mathutil::gettimeofday_sec() ;
907 convertArray( spectra, spectraF_ ) ;
908 toInt( flagtraUC_, flagtra ) ;
909 toInt( rflagUI_, rflag ) ;
910 t1 = mathutil::gettimeofday_sec() ;
911 eToInt = t1 - t0 ;
912
913 return nrow ;
914}
915
916Int STGrid::getDataChunk( Array<Float> &spectra,
917 Array<Double> &direction,
918 Array<uChar> &flagtra,
919 Array<uInt> &rflag,
920 Array<Float> &weight )
921{
922 LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
923 Int nrow = spectra.shape()[1] ;
924 Int remainingRow = nrow_ - nprocessed_ ;
925 if ( remainingRow < nrow ) {
926 nrow = remainingRow ;
927 IPosition mshape( 2, nchan_, nrow ) ;
928 IPosition vshape( 1, nrow ) ;
929 spectra.resize( mshape ) ;
930 flagtra.resize( mshape ) ;
931 direction.resize( IPosition(2,2,nrow) ) ;
932 rflag.resize( vshape ) ;
933 weight.resize( mshape ) ;
934 }
935 // 2011/12/22 TN
936 // tsys shares its storage with weight
937 Array<Float> tsys( weight ) ;
938 Array<Double> tint( rflag.shape() ) ;
939
940 Vector<uInt> rflagVec( rflag ) ;
941 Vector<Double> tintVec( tint ) ;
942
943 RefRows rows( nprocessed_, nprocessed_+nrow-1, 1 ) ;
944 //os<<LogIO::DEBUGGING<<"nprocessed_="<<nprocessed_<<": rows.nrows()="<<rows.nrows()<<LogIO::POST ;
945 spectraCol_.getColumnCells( rows, spectra ) ;
946 flagtraCol_.getColumnCells( rows, flagtra ) ;
947 directionCol_.getColumnCells( rows, direction ) ;
948 flagRowCol_.getColumnCells( rows, rflagVec ) ;
949 intervalCol_.getColumnCells( rows, tintVec ) ;
950 Vector<Float> tsysTemp = tsysCol_( nprocessed_ ) ;
951 if ( tsysTemp.nelements() == (uInt)nchan_ )
952 tsysCol_.getColumnCells( rows, tsys ) ;
953 else
954 tsys = tsysTemp[0] ;
955
956 double t0,t1 ;
957 t0 = mathutil::gettimeofday_sec() ;
958 getWeight( weight, tsys, tint ) ;
959 t1 = mathutil::gettimeofday_sec() ;
960 eGetWeight += t1-t0 ;
961
962 nprocessed_ += nrow ;
963
964 return nrow ;
965}
966
967void STGrid::setupArray()
968{
969 LogIO os( LogOrigin("STGrid","setupArray",WHERE) ) ;
970 ROScalarColumn<uInt> polnoCol( tableList_[0], "POLNO" ) ;
971 Vector<uInt> pols = polnoCol.getColumn() ;
972 //os << pols << LogIO::POST ;
973 Vector<uInt> pollistOrg ;
974 npolOrg_ = 0 ;
975 uInt polno ;
976 for ( uInt i = 0 ; i < polnoCol.nrow() ; i++ ) {
977 polno = pols( i ) ;
978 if ( allNE( pollistOrg, polno ) ) {
979 pollistOrg.resize( npolOrg_+1, True ) ;
980 pollistOrg[npolOrg_] = polno ;
981 npolOrg_++ ;
982 }
983 }
984 if ( pollist_.size() == 0 )
985 pollist_ = pollistOrg ;
986 else {
987 Vector<uInt> newlist ;
988 uInt newsize = 0 ;
989 for ( uInt i = 0 ; i < pollist_.size() ; i++ ) {
990 if ( anyEQ( pollistOrg, pollist_[i] ) ) {
991 newlist.resize( newsize+1, True ) ;
992 newlist[newsize] = pollist_[i] ;
993 newsize++ ;
994 }
995 }
996 pollist_.assign( newlist ) ;
997 }
998 npol_ = pollist_.size() ;
999 if ( npol_ == 0 ) {
1000 os << LogIO::SEVERE << "Empty pollist" << LogIO::EXCEPTION ;
1001 }
1002 rows_.resize( nfile_ ) ;
1003 for ( uInt i = 0 ; i < nfile_ ; i++ ) {
1004 rows_[i] = tableList_[i].nrow() / npolOrg_ ;
1005 }
1006 flagtraCol_.attach( tableList_[0], "FLAGTRA" ) ;
1007 nchan_ = flagtraCol_( 0 ).nelements() ;
1008// os << "npol_ = " << npol_ << "(" << pollist_ << ")" << endl
1009// << "nchan_ = " << nchan_ << endl
1010// << "nrow_ = " << nrow_ << LogIO::POST ;
1011}
1012
1013void STGrid::getWeight( Array<Float> &w,
1014 Array<Float> &tsys,
1015 Array<Double> &tint )
1016{
1017 LogIO os( LogOrigin("STGrid","getWeight",WHERE) ) ;
1018
1019 // 2011/12/22 TN
1020 // w (weight) and tsys share storage
1021 IPosition refShape = tsys.shape() ;
1022 Int nchan = refShape[0] ;
1023 Int nrow = refShape[1] ;
1024// os << "nchan=" << nchan << ", nrow=" << nrow << LogIO::POST ;
1025// os << "w.shape()=" << w.shape() << endl
1026// << "tsys.shape()=" << tsys.shape() << endl
1027// << "tint.shape()=" << tint.shape() << LogIO::POST ;
1028
1029 // set weight
1030 if ( wtype_.compare( "UNIFORM" ) == 0 ) {
1031 w = 1.0 ;
1032 }
1033 else if ( wtype_.compare( "TINT" ) == 0 ) {
1034 Bool b0, b1 ;
1035 Float *w_p = w.getStorage( b0 ) ;
1036 Float *w0_p = w_p ;
1037 const Double *ti_p = tint.getStorage( b1 ) ;
1038 const Double *w1_p = ti_p ;
1039 for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1040 for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1041 *w0_p = *w1_p ;
1042 w0_p++ ;
1043 }
1044 w1_p++ ;
1045 }
1046 w.putStorage( w_p, b0 ) ;
1047 tint.freeStorage( ti_p, b1 ) ;
1048 }
1049 else if ( wtype_.compare( "TSYS" ) == 0 ) {
1050 Bool b0 ;
1051 Float *w_p = w.getStorage( b0 ) ;
1052 Float *w0_p = w_p ;
1053 for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1054 for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1055 Float temp = *w0_p ;
1056 *w0_p = 1.0 / ( temp * temp ) ;
1057 w0_p++ ;
1058 }
1059 }
1060 w.putStorage( w_p, b0 ) ;
1061 }
1062 else if ( wtype_.compare( "TINTSYS" ) == 0 ) {
1063 Bool b0, b1 ;
1064 Float *w_p = w.getStorage( b0 ) ;
1065 Float *w0_p = w_p ;
1066 const Double *ti_p = tint.getStorage( b1 ) ;
1067 const Double *w1_p = ti_p ;
1068 for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1069 Float interval = *w1_p ;
1070 for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1071 Float temp = *w0_p ;
1072 *w0_p = interval / ( temp * temp ) ;
1073 w0_p++ ;
1074 }
1075 w1_p++ ;
1076 }
1077 w.putStorage( w_p, b0 ) ;
1078 tint.freeStorage( ti_p, b1 ) ;
1079 }
1080 else {
1081 //LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
1082 //os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
1083 w = 1.0 ;
1084 }
1085}
1086
1087void STGrid::toInt( Array<uChar> &u, Array<Int> &v )
1088{
1089 uInt len = u.nelements() ;
1090 Int *int_p = new Int[len] ;
1091 Bool deleteIt ;
1092 const uChar *data_p = u.getStorage( deleteIt ) ;
1093 Int *i_p = int_p ;
1094 const uChar *u_p = data_p ;
1095 for ( uInt i = 0 ; i < len ; i++ ) {
1096 *i_p = ( *u_p == 0 ) ? 0 : 1 ;
1097 i_p++ ;
1098 u_p++ ;
1099 }
1100 u.freeStorage( data_p, deleteIt ) ;
1101 v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
1102}
1103
1104void STGrid::toInt( Array<uInt> &u, Array<Int> &v )
1105{
1106 uInt len = u.nelements() ;
1107 Int *int_p = new Int[len] ;
1108 Bool deleteIt ;
1109 const uInt *data_p = u.getStorage( deleteIt ) ;
1110 Int *i_p = int_p ;
1111 const uInt *u_p = data_p ;
1112 for ( uInt i = 0 ; i < len ; i++ ) {
1113 *i_p = ( *u_p == 0 ) ? 0 : 1 ;
1114 i_p++ ;
1115 u_p++ ;
1116 }
1117 u.freeStorage( data_p, deleteIt ) ;
1118 v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
1119}
1120
1121void STGrid::toPixel( Array<Double> &world, Array<Double> &pixel )
1122{
1123 // gridding will be done on (nx_+2*convSupport_) x (ny_+2*convSupport_)
1124 // grid plane to avoid unexpected behavior on grid edge
1125 Block<Double> pixc( 2 ) ;
1126 pixc[0] = Double( nx_-1 ) * 0.5 ;
1127 pixc[1] = Double( ny_-1 ) * 0.5 ;
1128// pixc[0] = Double( nx_+2*convSupport_-1 ) * 0.5 ;
1129// pixc[1] = Double( ny_+2*convSupport_-1 ) * 0.5 ;
1130 uInt nrow = world.shape()[1] ;
1131 Bool bw, bp ;
1132 const Double *w_p = world.getStorage( bw ) ;
1133 Double *p_p = pixel.getStorage( bp ) ;
1134 const Double *ww_p = w_p ;
1135 Double *wp_p = p_p ;
1136 for ( uInt i = 0 ; i < nrow ; i++ ) {
1137 *wp_p = pixc[0] + ( *ww_p - center_[0] ) / cellx_ ;
1138 wp_p++ ;
1139 ww_p++ ;
1140 *wp_p = pixc[1] + ( *ww_p - center_[1] ) / celly_ ;
1141 wp_p++ ;
1142 ww_p++ ;
1143 }
1144 world.freeStorage( w_p, bw ) ;
1145 pixel.putStorage( p_p, bp ) ;
1146// String gridfile = "grid."+convType_+"."+String::toString(convSupport_)+".dat" ;
1147// ofstream ofs( gridfile.c_str(), ios::out ) ;
1148// ofs << "center " << center_(0) << " " << pixc(0)
1149// << " " << center_(1) << " " << pixc(1) << endl ;
1150// for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
1151// ofs << irow ;
1152// for ( uInt i = 0 ; i < 2 ; i++ ) {
1153// ofs << " " << world(i, irow) << " " << pixel(i, irow) ;
1154// }
1155// ofs << endl ;
1156// }
1157// ofs.close() ;
1158}
1159
1160void STGrid::boxFunc( Vector<Float> &convFunc, Int &convSize )
1161{
1162 convFunc = 0.0 ;
1163 for ( Int i = 0 ; i < convSize/2 ; i++ )
1164 convFunc(i) = 1.0 ;
1165}
1166
1167#define NEED_UNDERSCORES
1168#if defined(NEED_UNDERSCORES)
1169#define grdsf grdsf_
1170#endif
1171extern "C" {
1172 void grdsf(Double*, Double*);
1173}
1174void STGrid::spheroidalFunc( Vector<Float> &convFunc )
1175{
1176 convFunc = 0.0 ;
1177 for ( Int i = 0 ; i < convSampling_*convSupport_ ; i++ ) {
1178 Double nu = Double(i) / Double(convSupport_*convSampling_) ;
1179 Double val ;
1180 grdsf( &nu, &val ) ;
1181 convFunc(i) = ( 1.0 - nu * nu ) * val ;
1182 }
1183}
1184
1185void STGrid::gaussFunc( Vector<Float> &convFunc )
1186{
1187 convFunc = 0.0 ;
1188 // HWHM of the Gaussian is convSupport_ / 4
1189 // To take into account Gaussian tail, kernel cutoff is set to 4 * HWHM
1190 Int len = convSampling_ * convSupport_ ;
1191 Double hwhm = len * 0.25 ;
1192 for ( Int i = 0 ; i < len ; i++ ) {
1193 Double val = Double(i) / hwhm ;
1194 convFunc(i) = exp( -log(2)*val*val ) ;
1195 }
1196}
1197
1198void STGrid::pbFunc( Vector<Float> &convFunc )
1199{
1200 convFunc = 0.0 ;
1201}
1202
1203void STGrid::setConvFunc( Vector<Float> &convFunc )
1204{
1205 convSupport_ = userSupport_ ;
1206 if ( convType_ == "BOX" ) {
1207 if ( convSupport_ < 0 )
1208 convSupport_ = 0 ;
1209 Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
1210 convFunc.resize( convSize ) ;
1211 boxFunc( convFunc, convSize ) ;
1212 }
1213 else if ( convType_ == "SF" ) {
1214 if ( convSupport_ < 0 )
1215 convSupport_ = 3 ;
1216 Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
1217 convFunc.resize( convSize ) ;
1218 spheroidalFunc( convFunc ) ;
1219 }
1220 else if ( convType_ == "GAUSS" ) {
1221 // to take into account Gaussian tail
1222 if ( convSupport_ < 0 )
1223 convSupport_ = 12 ; // 3 * 4
1224 else {
1225 convSupport_ = userSupport_ * 4 ;
1226 }
1227 Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
1228 convFunc.resize( convSize ) ;
1229 gaussFunc( convFunc ) ;
1230 }
1231 else if ( convType_ == "PB" ) {
1232 if ( convSupport_ < 0 )
1233 convSupport_ = 0 ;
1234 pbFunc( convFunc ) ;
1235 }
1236 else {
1237 throw AipsError( "Unsupported convolution function" ) ;
1238 }
1239}
1240
1241string STGrid::saveData( string outfile )
1242{
1243 LogIO os( LogOrigin("STGrid", "saveData", WHERE) ) ;
1244 double t0, t1 ;
1245 t0 = mathutil::gettimeofday_sec() ;
1246
1247 String outfile_ ;
1248 if ( outfile.size() == 0 ) {
1249 if ( infileList_[0].lastchar() == '/' ) {
1250 outfile_ = infileList_[0].substr( 0, infileList_[0].size()-1 ) ;
1251 }
1252 else {
1253 outfile_ = infileList_[0] ;
1254 }
1255 outfile_ += ".grid" ;
1256 }
1257 else {
1258 outfile_ = outfile ;
1259 }
1260 Table tab ;
1261 prepareTable( tab, outfile_ ) ;
1262 IPosition dshape = data_.shape() ;
1263 Int nrow = nx_ * ny_ * npol_ ;
1264 tab.rwKeywordSet().define( "nPol", npol_ ) ;
1265 tab.addRow( nrow ) ;
1266 Vector<Double> cpix( 2 ) ;
1267 cpix(0) = Double( nx_ - 1 ) * 0.5 ;
1268 cpix(1) = Double( ny_ - 1 ) * 0.5 ;
1269 Vector<Double> dir( 2 ) ;
1270 ArrayColumn<Double> directionCol( tab, "DIRECTION" ) ;
1271 ArrayColumn<Float> spectraCol( tab, "SPECTRA" ) ;
1272 ScalarColumn<uInt> polnoCol( tab, "POLNO" ) ;
1273 Int irow = 0 ;
1274 Vector<Float> sp( nchan_ ) ;
1275 Bool bsp, bdata ;
1276 const Float *data_p = data_.getStorage( bdata ) ;
1277 Float *wsp_p, *sp_p ;
1278 const Float *wdata_p = data_p ;
1279 long step = nx_ * ny_ * npol_ ;
1280 long offset ;
1281 for ( Int iy = 0 ; iy < ny_ ; iy++ ) {
1282 dir(1) = center_(1) - ( cpix(1) - (Double)iy ) * celly_ ;
1283 for ( Int ix = 0 ; ix < nx_ ; ix++ ) {
1284 dir(0) = center_(0) - ( cpix(0) - (Double)ix ) * cellx_ ;
1285 for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
1286 offset = ix + iy * nx_ + ipol * nx_ * ny_ ;
1287 //os << "offset = " << offset << LogIO::POST ;
1288 sp_p = sp.getStorage( bsp ) ;
1289 wsp_p = sp_p ;
1290 wdata_p = data_p + offset ;
1291 for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
1292 *wsp_p = *wdata_p ;
1293 wsp_p++ ;
1294 wdata_p += step ;
1295 }
1296 sp.putStorage( sp_p, bsp ) ;
1297 spectraCol.put( irow, sp ) ;
1298 directionCol.put( irow, dir ) ;
1299 polnoCol.put( irow, pollist_[ipol] ) ;
1300 irow++ ;
1301 }
1302 }
1303 }
1304 data_.freeStorage( data_p, bdata ) ;
1305
1306 t1 = mathutil::gettimeofday_sec() ;
1307 os << "saveData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
1308
1309 return outfile_ ;
1310}
1311
1312void STGrid::prepareTable( Table &tab, String &name )
1313{
1314 Table t( infileList_[0], Table::Old ) ;
1315 t.deepCopy( name, Table::New, False, t.endianFormat(), True ) ;
1316 tab = Table( name, Table::Update ) ;
1317}
1318
1319}
Note: See TracBrowser for help on using the repository browser.