source: trunk/src/MSFiller.cpp@ 2290

Last change on this file since 2290 was 2289, checked in by ShinnosukeKawakami, 13 years ago

merged parallel branch to trunk
File size: 68.5 KB
Line 
1//
2// C++ Interface: MSFiller
3//
4// Description:
5//
6// This class is specific filler for MS format
7//
8// Takeshi Nakazato <takeshi.nakazato@nao.ac.jp>, (C) 2010
9//
10// Copyright: See COPYING file that comes with this distribution
11//
12//
13
14#include <iostream>
15#include <map>
16
17#include <tables/Tables/ExprNode.h>
18#include <tables/Tables/TableIter.h>
19#include <tables/Tables/TableColumn.h>
20#include <tables/Tables/ScalarColumn.h>
21#include <tables/Tables/ArrayColumn.h>
22#include <tables/Tables/TableParse.h>
23#include <tables/Tables/TableRow.h>
24
25#include <casa/Containers/RecordField.h>
26#include <casa/Logging/LogIO.h>
27#include <casa/Arrays/Slicer.h>
28#include <casa/Quanta/MVTime.h>
29#include <casa/OS/Path.h>
30
31#include <measures/Measures/Stokes.h>
32#include <measures/Measures/MEpoch.h>
33#include <measures/Measures/MCEpoch.h>
34#include <measures/Measures/MFrequency.h>
35#include <measures/Measures/MCFrequency.h>
36#include <measures/Measures/MPosition.h>
37#include <measures/Measures/MCPosition.h>
38#include <measures/Measures/MDirection.h>
39#include <measures/Measures/MCDirection.h>
40#include <measures/Measures/MeasConvert.h>
41#include <measures/TableMeasures/ScalarMeasColumn.h>
42#include <measures/TableMeasures/ArrayMeasColumn.h>
43#include <measures/TableMeasures/ScalarQuantColumn.h>
44#include <measures/TableMeasures/ArrayQuantColumn.h>
45
46#include <ms/MeasurementSets/MSAntennaIndex.h>
47
48#include <atnf/PKSIO/SrcType.h>
49
50#include "MSFiller.h"
51#include "STHeader.h"
52
53// #include <ctime>
54// #include <sys/time.h>
55
56#include "MathUtils.h"
57
58using namespace casa ;
59using namespace std ;
60
61namespace asap {
62// double MSFiller::gettimeofday_sec()
63// {
64// struct timeval tv ;
65// gettimeofday( &tv, NULL ) ;
66// return tv.tv_sec + (double)tv.tv_usec*1.0e-6 ;
67// }
68
69MSFiller::MSFiller( casa::CountedPtr<Scantable> stable )
70 : table_( stable ),
71 tablename_( "" ),
72 antenna_( -1 ),
73 antennaStr_(""),
74 getPt_( True ),
75 isFloatData_( False ),
76 isData_( False ),
77 isDoppler_( False ),
78 isFlagCmd_( False ),
79 isFreqOffset_( False ),
80 isHistory_( False ),
81 isProcessor_( False ),
82 isSysCal_( False ),
83 isWeather_( False ),
84 colTsys_( "TSYS_SPECTRUM" ),
85 colTcal_( "TCAL_SPECTRUM" )
86{
87 os_ = LogIO() ;
88 os_.origin( LogOrigin( "MSFiller", "MSFiller()", WHERE ) ) ;
89}
90
91MSFiller::~MSFiller()
92{
93 os_.origin( LogOrigin( "MSFiller", "~MSFiller()", WHERE ) ) ;
94}
95
96bool MSFiller::open( const std::string &filename, const casa::Record &rec )
97{
98 os_.origin( LogOrigin( "MSFiller", "open()", WHERE ) ) ;
99 //double startSec = mathutil::gettimeofday_sec() ;
100 //os_ << "start MSFiller::open() startsec=" << startSec << LogIO::POST ;
101 //os_ << " filename = " << filename << endl ;
102
103 // parsing MS options
104 if ( rec.isDefined( "ms" ) ) {
105 Record msrec = rec.asRecord( "ms" ) ;
106 if ( msrec.isDefined( "getpt" ) ) {
107 getPt_ = msrec.asBool( "getpt" ) ;
108 }
109 if ( msrec.isDefined( "antenna" ) ) {
110 if ( msrec.type( msrec.fieldNumber( "antenna" ) ) == TpInt ) {
111 antenna_ = msrec.asInt( "antenna" ) ;
112 }
113 else {
114 //antenna_ = atoi( msrec.asString( "antenna" ).c_str() ) ;
115 antennaStr_ = msrec.asString( "antenna" ) ;
116 }
117 }
118 else {
119 antenna_ = 0 ;
120 }
121 }
122
123 MeasurementSet *tmpMS = new MeasurementSet( filename, Table::Old ) ;
124 //mstable_ = (*tmpMS)( tmpMS->col("ANTENNA1") == antenna_
125 // && tmpMS->col("ANTENNA1") == tmpMS->col("ANTENNA2") ) ;
126 tablename_ = tmpMS->tableName() ;
127 if ( antenna_ == -1 && antennaStr_.size() > 0 ) {
128 MSAntennaIndex msAntIdx( tmpMS->antenna() ) ;
129 Vector<Int> id = msAntIdx.matchAntennaName( antennaStr_ ) ;
130 if ( id.size() > 0 )
131 antenna_ = id[0] ;
132 }
133
134 os_ << "Parsing MS options" << endl ;
135 os_ << " getPt = " << getPt_ << endl ;
136 os_ << " antenna = " << antenna_ << endl ;
137 os_ << " antennaStr = " << antennaStr_ << LogIO::POST ;
138
139 mstable_ = MeasurementSet( (*tmpMS)( tmpMS->col("ANTENNA1") == antenna_
140 && tmpMS->col("ANTENNA1") == tmpMS->col("ANTENNA2") ) ) ;
141// stringstream ss ;
142// ss << "SELECT FROM $1 WHERE ANTENNA1 == ANTENNA2 && ANTENNA1 == " << antenna_ ;
143// String taql( ss.str() ) ;
144// mstable_ = MeasurementSet( tableCommand( taql, *tmpMS ) ) ;
145 delete tmpMS ;
146
147 // check which data column exists
148 isFloatData_ = mstable_.tableDesc().isColumn( "FLOAT_DATA" ) ;
149 isData_ = mstable_.tableDesc().isColumn( "DATA" ) ;
150
151 //double endSec = mathutil::gettimeofday_sec() ;
152 //os_ << "end MSFiller::open() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
153 return true ;
154}
155
156void MSFiller::fill()
157{
158 os_.origin( LogOrigin( "MSFiller", "fill()", WHERE ) ) ;
159 //double startSec = mathutil::gettimeofday_sec() ;
160 //os_ << "start MSFiller::fill() startSec=" << startSec << LogIO::POST ;
161
162 //double time0 = mathutil::gettimeofday_sec() ;
163 //os_ << "start init fill: " << time0 << LogIO::POST ;
164
165 // Initialize header
166 STHeader sdh ;
167 initHeader( sdh ) ;
168
169 // check if optional table exists
170 //const TableRecord msrec = tablesel_.keywordSet() ;
171 const TableRecord msrec = mstable_.keywordSet() ;
172 isDoppler_ = msrec.isDefined( "DOPPLER" ) ;
173 if ( isDoppler_ )
174 if ( mstable_.doppler().nrow() == 0 )
175 isDoppler_ = False ;
176 isFlagCmd_ = msrec.isDefined( "FLAG_CMD" ) ;
177 if ( isFlagCmd_ )
178 if ( mstable_.flagCmd().nrow() == 0 )
179 isFlagCmd_ = False ;
180 isFreqOffset_ = msrec.isDefined( "FREQ_OFFSET" ) ;
181 if ( isFreqOffset_ )
182 if ( mstable_.freqOffset().nrow() == 0 )
183 isFreqOffset_ = False ;
184 isHistory_ = msrec.isDefined( "HISTORY" ) ;
185 if ( isHistory_ )
186 if ( mstable_.history().nrow() == 0 )
187 isHistory_ = False ;
188 isProcessor_ = msrec.isDefined( "PROCESSOR" ) ;
189 if ( isProcessor_ )
190 if ( mstable_.processor().nrow() == 0 )
191 isProcessor_ = False ;
192 isSysCal_ = msrec.isDefined( "SYSCAL" ) ;
193 if ( isSysCal_ )
194 if ( mstable_.sysCal().nrow() == 0 )
195 isSysCal_ = False ;
196 isWeather_ = msrec.isDefined( "WEATHER" ) ;
197 if ( isWeather_ )
198 if ( mstable_.weather().nrow() == 0 )
199 isWeather_ = False ;
200
201 // Access to MS subtables
202 MSField fieldtab = mstable_.field() ;
203 MSPolarization poltab = mstable_.polarization() ;
204 MSDataDescription ddtab = mstable_.dataDescription() ;
205 MSObservation obstab = mstable_.observation() ;
206 MSSource srctab = mstable_.source() ;
207 MSSpectralWindow spwtab = mstable_.spectralWindow() ;
208 MSSysCal caltab = mstable_.sysCal() ;
209 if ( caltab.nrow() == 0 )
210 isSysCal_ = False ;
211 else {
212 if ( !caltab.tableDesc().isColumn( colTcal_ ) ) {
213 colTcal_ = "TCAL" ;
214 if ( !caltab.tableDesc().isColumn( colTcal_ ) )
215 colTcal_ = "NONE" ;
216 }
217 if ( !caltab.tableDesc().isColumn( colTsys_ ) ) {
218 colTsys_ = "TSYS" ;
219 if ( !caltab.tableDesc().isColumn( colTcal_ ) )
220 colTsys_ = "NONE" ;
221 }
222 }
223// colTcal_ = "TCAL" ;
224// colTsys_ = "TSYS" ;
225 MSPointing pointtab = mstable_.pointing() ;
226 if ( mstable_.weather().nrow() == 0 )
227 isWeather_ = False ;
228 MSState stattab = mstable_.state() ;
229 MSAntenna anttab = mstable_.antenna() ;
230
231 // TEST
232 // memory allocation by boost::object_pool
233 boost::object_pool<ROTableColumn> *tpoolr = new boost::object_pool<ROTableColumn> ;
234 //
235
236 // SUBTABLES: FREQUENCIES
237 //string freqFrame = getFrame() ;
238 string freqFrame = "LSRK" ;
239 table_->frequencies().setFrame( freqFrame ) ;
240 table_->frequencies().setFrame( freqFrame, True ) ;
241
242 // SUBTABLES: WEATHER
243 fillWeather() ;
244
245 // SUBTABLES: FOCUS
246 fillFocus() ;
247
248 // SUBTABLES: TCAL
249 fillTcal( tpoolr ) ;
250
251 // SUBTABLES: FIT
252 //fillFit() ;
253
254 // SUBTABLES: HISTORY
255 //fillHistory() ;
256
257 // MAIN
258 // Iterate over several ids
259 map<Int, uInt> ifmap ; // (IFNO, FREQ_ID) pair
260 ROArrayQuantColumn<Double> *sharedQDArrCol = new ROArrayQuantColumn<Double>( anttab, "POSITION" ) ;
261 Vector< Quantum<Double> > antpos = (*sharedQDArrCol)( antenna_ ) ;
262 delete sharedQDArrCol ;
263 MPosition mp( MVPosition( antpos ), MPosition::ITRF ) ;
264 Vector<Double> pt ;
265 ROArrayColumn<Double> pdcol ;
266 Vector<Double> defaultScanrate( 2, 0.0 ) ;
267 if ( getPt_ ) {
268 pointtab = MSPointing( pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort("TIME") ) ;
269 ROScalarColumn<Double> ptcol( pointtab, "TIME" ) ;
270 ptcol.getColumn( pt ) ;
271 TableRecord trec = ptcol.keywordSet() ;
272 String tUnit = trec.asArrayString( "QuantumUnits" ).data()[0] ;
273 if ( tUnit == "d" )
274 pt *= 86400.0 ;
275 pdcol.attach( pointtab, "DIRECTION" ) ;
276 }
277 String stationName = asString( "STATION", antenna_, anttab, tpoolr ) ;
278 String antennaName = asString( "NAME", antenna_, anttab, tpoolr ) ;
279 sdh.antennaposition.resize( 3 ) ;
280 for ( int i = 0 ; i < 3 ; i++ )
281 sdh.antennaposition[i] = antpos[i].getValue( "m" ) ;
282 String telescopeName = "" ;
283
284 //double time1 = mathutil::gettimeofday_sec() ;
285 //os_ << "end init fill: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
286
287 // row based
288 Table &stab = table_->table() ;
289 TableRow row( stab ) ;
290 TableRecord &trec = row.record() ;
291 RecordFieldPtr< Array<Float> > spRF( trec, "SPECTRA" ) ;
292 RecordFieldPtr< Array<uChar> > ucarrRF( trec, "FLAGTRA" ) ;
293 RecordFieldPtr<Double> timeRF( trec, "TIME" ) ;
294 RecordFieldPtr< Array<Float> > tsysRF( trec, "TSYS" ) ;
295 RecordFieldPtr<Double> intervalRF( trec, "INTERVAL" ) ;
296 RecordFieldPtr< Array<Double> > dirRF( trec, "DIRECTION" ) ;
297 RecordFieldPtr<Float> azRF( trec, "AZIMUTH" ) ;
298 RecordFieldPtr<Float> elRF( trec, "ELEVATION" ) ;
299 RecordFieldPtr< Array<Double> > scrRF( trec, "SCANRATE" ) ;
300 RecordFieldPtr<uInt> cycleRF( trec, "CYCLENO" ) ;
301 RecordFieldPtr<uInt> flrRF( trec, "FLAGROW" ) ;
302 RecordFieldPtr<uInt> tcalidRF( trec, "TCAL_ID" ) ;
303 RecordFieldPtr<uInt> widRF( trec, "WEATHER_ID" ) ;
304 RecordFieldPtr<uInt> polnoRF( trec, "POLNO" ) ;
305 RecordFieldPtr<Int> refbRF( trec, "REFBEAMNO" ) ;
306 RecordFieldPtr<Int> fitidRF( trec, "FIT_ID" ) ;
307 RecordFieldPtr<Float> tauRF( trec, "OPACITY" ) ;
308 RecordFieldPtr<uInt> beamRF( trec, "BEAMNO" ) ;
309 RecordFieldPtr<uInt> focusidRF( trec, "FOCUS_ID" ) ;
310 RecordFieldPtr<uInt> ifnoRF( trec, "IFNO" ) ;
311 RecordFieldPtr<uInt> molidRF( trec, "MOLECULE_ID" ) ;
312 RecordFieldPtr<uInt> freqidRF( trec, "FREQ_ID" ) ;
313 RecordFieldPtr<uInt> scanRF( trec, "SCANNO" ) ;
314 RecordFieldPtr<Int> srctypeRF( trec, "SRCTYPE" ) ;
315 RecordFieldPtr<String> srcnameRF( trec, "SRCNAME" ) ;
316 RecordFieldPtr<String> fieldnameRF( trec, "FIELDNAME" ) ;
317 RecordFieldPtr< Array<Double> > srcpmRF( trec, "SRCPROPERMOTION" ) ;
318 RecordFieldPtr< Array<Double> > srcdirRF( trec, "SRCDIRECTION" ) ;
319 RecordFieldPtr<Double> sysvelRF( trec, "SRCVELOCITY" ) ;
320
321 // REFBEAMNO
322 *refbRF = -1 ;
323
324 // FIT_ID
325 *fitidRF = -1 ;
326
327 // OPACITY
328 *tauRF = 0.0 ;
329
330 //
331 // ITERATION: OBSERVATION_ID
332 //
333 TableIterator iter0( mstable_, "OBSERVATION_ID" ) ;
334 while( !iter0.pastEnd() ) {
335 //time0 = mathutil::gettimeofday_sec() ;
336 //os_ << "start 0th iteration: " << time0 << LogIO::POST ;
337 Table t0 = iter0.table() ;
338 Int obsId = asInt( "OBSERVATION_ID", 0, t0, tpoolr ) ;
339 if ( sdh.observer == "" ) {
340 sdh.observer = asString( "OBSERVER", obsId, obstab, tpoolr ) ;
341 }
342 if ( sdh.project == "" ) {
343 sdh.project = asString( "PROJECT", obsId, obstab, tpoolr ) ;
344 }
345 ROArrayMeasColumn<MEpoch> *tmpMeasCol = new ROArrayMeasColumn<MEpoch>( obstab, "TIME_RANGE" ) ;
346 MEpoch me = (*tmpMeasCol)( obsId )( IPosition(1,0) ) ;
347 delete tmpMeasCol ;
348 if ( sdh.utc == 0.0 ) {
349 sdh.utc = me.get( "d" ).getValue() ;
350 }
351 if ( telescopeName == "" ) {
352 telescopeName = asString( "TELESCOPE_NAME", obsId, obstab, tpoolr ) ;
353 }
354 Int nbeam = 0 ;
355 //time1 = mathutil::gettimeofday_sec() ;
356 //os_ << "end 0th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
357 //
358 // ITERATION: FEED1
359 //
360 TableIterator iter1( t0, "FEED1" ) ;
361 while( !iter1.pastEnd() ) {
362 //time0 = mathutil::gettimeofday_sec() ;
363 //os_ << "start 1st iteration: " << time0 << LogIO::POST ;
364 Table t1 = iter1.table() ;
365 // assume FEED1 == FEED2
366 Int feedId = asInt( "FEED1", 0, t1, tpoolr ) ;
367 nbeam++ ;
368
369 // BEAMNO
370 *beamRF = feedId ;
371
372 // FOCUS_ID
373 *focusidRF = 0 ;
374
375 //time1 = mathutil::gettimeofday_sec() ;
376 //os_ << "end 1st iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
377 //
378 // ITERATION: FIELD_ID
379 //
380 TableIterator iter2( t1, "FIELD_ID" ) ;
381 while( !iter2.pastEnd() ) {
382 //time0 = mathutil::gettimeofday_sec() ;
383 //os_ << "start 2nd iteration: " << time0 << LogIO::POST ;
384 Table t2 = iter2.table() ;
385 Int fieldId = asInt( "FIELD_ID", 0, t2, tpoolr ) ;
386 Int srcId = asInt( "SOURCE_ID", fieldId, fieldtab, tpoolr ) ;
387 String fieldName = asString( "NAME", fieldId, fieldtab, tpoolr ) ;
388 fieldName += "__" + String::toString(fieldId) ;
389 ROArrayMeasColumn<MDirection> *delayDirCol = new ROArrayMeasColumn<MDirection>( fieldtab, "DELAY_DIR" ) ;
390 Vector<MDirection> delayDir = (*delayDirCol)( fieldId ) ;
391 delete delayDirCol ;
392
393 // FIELDNAME
394 *fieldnameRF = fieldName ;
395
396
397 //time1 = mathutil::gettimeofday_sec() ;
398 //os_ << "end 2nd iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
399 //
400 // ITERATION: DATA_DESC_ID
401 //
402 TableIterator iter3( t2, "DATA_DESC_ID" ) ;
403 while( !iter3.pastEnd() ) {
404 //time0 = mathutil::gettimeofday_sec() ;
405 //os_ << "start 3rd iteration: " << time0 << LogIO::POST ;
406 Table t3 = iter3.table() ;
407 Int ddId = asInt( "DATA_DESC_ID", 0, t3, tpoolr ) ;
408 Int polId = asInt( "POLARIZATION_ID", ddId, ddtab, tpoolr ) ;
409 Int spwId = asInt( "SPECTRAL_WINDOW_ID", ddId, ddtab, tpoolr ) ;
410
411 // IFNO
412 *ifnoRF = (uInt)spwId ;
413
414 // polarization information
415 Int npol = asInt( "NUM_CORR", polId, poltab, tpoolr ) ;
416 ROArrayColumn<Int> *roArrICol = new ROArrayColumn<Int>( poltab, "CORR_TYPE" ) ;
417 Vector<Int> corrtype = (*roArrICol)( polId ) ;
418 delete roArrICol ;
419 String srcName( "" ) ;
420 Vector<Double> srcPM( 2, 0.0 ) ;
421 Vector<Double> srcDir( 2, 0.0 ) ;
422 MDirection md ;
423 Vector<Double> restFreqs ;
424 Vector<String> transitionName ;
425 Vector<Double> sysVels ;
426// os_ << "npol = " << npol << LogIO::POST ;
427// os_ << "corrtype = " << corrtype << LogIO::POST ;
428
429 // source and molecular transition
430 sourceInfo( srcId, spwId, srcName, md, srcPM, restFreqs, transitionName, sysVels, tpoolr ) ;
431// os_ << "srcId = " << srcId << ", spwId = " << spwId << LogIO::POST ;
432
433 // SRCNAME
434 *srcnameRF = srcName ;
435
436// os_ << "srcName = " << srcName << LogIO::POST ;
437
438 // SRCPROPERMOTION
439 *srcpmRF = srcPM ;
440
441 //os_ << "srcPM = " << srcPM << LogIO::POST ;
442
443 // SRCDIRECTION
444 *srcdirRF = md.getAngle().getValue( "rad" ) ;
445
446 //os_ << "srcDir = " << srcDir << LogIO::POST ;
447
448 // SRCVELOCITY
449 Double sysVel = 0.0 ;
450 if ( !sysVels.empty() )
451 sysVel = sysVels[0] ;
452 *sysvelRF = sysVel ;
453
454// os_ << "sysVel = " << sysVel << LogIO::POST ;
455
456 uInt molId = table_->molecules().addEntry( restFreqs, transitionName, transitionName ) ;
457
458 // MOLECULE_ID
459 *molidRF = molId ;
460
461 // spectral setup
462 uInt freqId ;
463 Int nchan = asInt( "NUM_CHAN", spwId, spwtab, tpoolr ) ;
464 Bool iswvr = False ;
465 if ( nchan == 4 ) iswvr = True ;
466 sdh.nchan = max( sdh.nchan, nchan ) ;
467 map<Int,uInt>::iterator iter = ifmap.find( spwId ) ;
468 if ( iter == ifmap.end() ) {
469 ROScalarMeasColumn<MEpoch> *tmpMeasCol = new ROScalarMeasColumn<MEpoch>( t3, "TIME" ) ;
470 me = (*tmpMeasCol)( 0 ) ;
471 delete tmpMeasCol ;
472 Double refpix ;
473 Double refval ;
474 Double increment ;
475 spectralSetup( spwId,
476 me,
477 mp,
478 md,
479 refpix,
480 refval,
481 increment,
482 nchan,
483 sdh.freqref,
484 sdh.reffreq,
485 sdh.bandwidth,
486 tpoolr ) ;
487 freqId = table_->frequencies().addEntry( refpix, refval, increment ) ;
488 ifmap.insert( pair<Int, uInt>(spwId,freqId) ) ;
489 //os_ << "added to ifmap: (" << spwId << "," << freqId << ")" << LogIO::POST ;
490 }
491 else {
492 freqId = iter->second ;
493 }
494
495 // FREQ_ID
496 *freqidRF = freqId ;
497
498 // for TSYS and TCAL
499 Vector<MEpoch> scTime ;
500 Vector<Double> scInterval ;
501 ROArrayColumn<Float> scTsysCol ;
502 MSSysCal caltabsel ;
503 if ( isSysCal_ ) {
504 caltabsel = caltab( caltab.col("ANTENNA_ID") == antenna_ && caltab.col("FEED_ID") == feedId && caltab.col("SPECTRAL_WINDOW_ID") == spwId ).sort("TIME") ;
505 ROScalarMeasColumn<MEpoch> scTimeCol( caltabsel, "TIME" ) ;
506 scTime.resize( caltabsel.nrow() ) ;
507 for ( uInt irow = 0 ; irow < caltabsel.nrow() ; irow++ )
508 scTime[irow] = scTimeCol( irow ) ;
509 ROScalarColumn<Double> scIntervalCol( caltabsel, "INTERVAL" ) ;
510 scIntervalCol.getColumn( scInterval ) ;
511 if ( colTsys_ != "NONE" )
512 scTsysCol.attach( caltabsel, colTsys_ ) ;
513 }
514
515 sdh.npol = max( sdh.npol, npol ) ;
516 if ( !iswvr && sdh.poltype == "" ) sdh.poltype = getPolType( corrtype[0] ) ;
517
518 //time1 = mathutil::gettimeofday_sec() ;
519 //os_ << "end 3rd iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
520 //
521 // ITERATION: SCAN_NUMBER
522 //
523 TableIterator iter4( t3, "SCAN_NUMBER" ) ;
524 while( !iter4.pastEnd() ) {
525 //time0 = mathutil::gettimeofday_sec() ;
526 //os_ << "start 4th iteration: " << time0 << LogIO::POST ;
527 Table t4 = iter4.table() ;
528 Int scanNum = asInt( "SCAN_NUMBER", 0, t4, tpoolr ) ;
529
530 // SCANNO
531 *scanRF = scanNum - 1 ;
532
533 uInt cycle = 0 ;
534
535 //time1 = mathutil::gettimeofday_sec() ;
536 //os_ << "end 4th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
537 //
538 // ITERATION: STATE_ID
539 //
540 TableIterator iter5( t4, "STATE_ID" ) ;
541 while( !iter5.pastEnd() ) {
542 //time0 = mathutil::gettimeofday_sec() ;
543 //os_ << "start 5th iteration: " << time0 << LogIO::POST ;
544 Table t5 = iter5.table() ;
545 Int stateId = asInt( "STATE_ID", 0, t5, tpoolr ) ;
546 String obstype = asString( "OBS_MODE", 0, stattab, tpoolr ) ;
547 if ( sdh.obstype == "" ) sdh.obstype = obstype ;
548
549 Int nrow = t5.nrow() ;
550 //time1 = mathutil::gettimeofday_sec() ;
551 //os_ << "end 5th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
552
553 Cube<Float> spArr ;
554 Cube<Bool> flArr ;
555 reshapeSpectraAndFlagtra( spArr,
556 flArr,
557 t5,
558 npol,
559 nchan,
560 nrow,
561 corrtype ) ;
562 if ( sdh.fluxunit == "" ) {
563 String colName = "FLOAT_DATA" ;
564 if ( isData_ ) colName = "DATA" ;
565 ROTableColumn dataCol( t5, colName ) ;
566 const TableRecord &dataColKeys = dataCol.keywordSet() ;
567 if ( dataColKeys.isDefined( "UNIT" ) )
568 sdh.fluxunit = dataColKeys.asString( "UNIT" ) ;
569 else if ( dataColKeys.isDefined( "QuantumUnits" ) )
570 sdh.fluxunit = dataColKeys.asString( "QuantumUnits" ) ;
571 }
572
573 ROScalarMeasColumn<MEpoch> *mTimeCol = new ROScalarMeasColumn<MEpoch>( t5, "TIME" ) ;
574 Block<MEpoch> mTimeB( nrow ) ;
575 for ( Int irow = 0 ; irow < nrow ; irow++ )
576 mTimeB[irow] = (*mTimeCol)( irow ) ;
577 Block<Int> sysCalIdx( nrow, -1 ) ;
578 if ( isSysCal_ ) {
579 getSysCalTime( scTime, scInterval, mTimeB, sysCalIdx ) ;
580 }
581 delete mTimeCol ;
582 Matrix<Float> defaulttsys( npol, 1, 1.0 ) ;
583 Int srcType = getSrcType( stateId, tpoolr ) ;
584 uInt diridx = 0 ;
585 uInt wid = 0 ;
586 Int pidx = 0 ;
587 Bool crossOK = False ;
588 Block<uInt> polnos( npol, 99 ) ;
589 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
590 Block<uInt> p = getPolNo( corrtype[ipol] ) ;
591 if ( p.size() > 1 ) {
592 if ( crossOK ) continue ;
593 else {
594 polnos[pidx] = p[0] ;
595 pidx++ ;
596 polnos[pidx] = p[1] ;
597 pidx++ ;
598 crossOK = True ;
599 }
600 }
601 else {
602 polnos[pidx] = p[0] ;
603 pidx++ ;
604 }
605 }
606
607 // SRCTYPE
608 *srctypeRF = srcType ;
609
610 for ( Int irow = 0 ; irow < nrow ; irow++ ) {
611 // CYCLENO
612 *cycleRF = cycle ;
613
614 // FLAGROW
615 *flrRF = (uInt)asBool( "FLAG_ROW", irow, t5, tpoolr ) ;
616
617 // SPECTRA, FLAG
618 Matrix<Float> sp = spArr.xyPlane( irow ) ;
619 Matrix<Bool> flb = flArr.xyPlane( irow ) ;
620 Matrix<uChar> fl( flb.shape() ) ;
621 convertArray( fl, flb ) ;
622
623 // TIME
624 *timeRF = mTimeB[irow].get("d").getValue() ;
625
626 // INTERVAL
627 *intervalRF = asDouble( "INTERVAL", irow, t5, tpoolr ) ;
628
629 // TSYS
630 Matrix<Float> tsys ;
631 if ( sysCalIdx[irow] != -1 && colTsys_ != "NONE" )
632 tsys = scTsysCol( sysCalIdx[irow] ) ;
633 else
634 tsys = defaulttsys ;
635
636 // TCAL_ID
637 Block<uInt> tcalids( npol, 0 ) ;
638 if ( sysCalIdx[irow] != -1 && colTcal_ != "NONE" ) {
639 tcalids = getTcalId( feedId, spwId, scTime[sysCalIdx[irow]] ) ;
640 }
641
642 // WEATHER_ID
643 if ( isWeather_ ) {
644 wid = getWeatherId( wid, mTimeB[irow].get("s").getValue() ) ;
645 *widRF = mwIndex_[wid] ;
646 }
647 else {
648 *widRF = wid ;
649 }
650
651
652 // DIRECTION, AZEL, SCANRATE
653 Vector<Double> dir ;
654 Vector<Double> azel ;
655 Vector<Double> scanrate = defaultScanrate ;
656 String refString ;
657 if ( getPt_ )
658 diridx = getDirection( diridx, dir, azel, scanrate, pt, pdcol, mTimeB[irow], mp ) ;
659 else
660 getSourceDirection( dir, azel, scanrate, mTimeB[irow], mp, delayDir ) ;
661 *dirRF = dir ;
662 *azRF = azel[0] ;
663 *elRF = azel[1] ;
664 *scrRF = scanrate ;
665
666 // Polarization dependent things
667 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
668 // POLNO
669 *polnoRF = polnos[ipol] ;
670
671 spRF.define( sp.row( ipol ) ) ;
672 ucarrRF.define( fl.row( ipol ) ) ;
673 tsysRF.define( tsys.row( ipol ) ) ;
674 *tcalidRF = tcalids[ipol] ;
675
676 // Commit row
677 stab.addRow() ;
678 row.put( stab.nrow()-1 ) ;
679 }
680
681 cycle++ ;
682 }
683
684 //time1 = mathutil::gettimeofday_sec() ;
685 //os_ << "end 5th iteration: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
686
687 iter5.next() ;
688 }
689 iter4.next() ;
690 }
691 iter3.next() ;
692 }
693 iter2.next() ;
694 }
695 iter1.next() ;
696 }
697 if ( sdh.nbeam < nbeam ) sdh.nbeam = nbeam ;
698
699 iter0.next() ;
700 }
701
702
703 delete tpoolr ;
704
705
706 // Table Keywords
707 sdh.nif = ifmap.size() ;
708 if ( ( telescopeName == "" ) || ( antennaName == telescopeName ) ) {
709 sdh.antennaname = antennaName ;
710 }
711 else {
712 sdh.antennaname = telescopeName + "//" + antennaName ;
713 }
714 if ( stationName != "" && stationName != antennaName ) {
715 sdh.antennaname += "@" + stationName ;
716 }
717 ROArrayColumn<Double> pdirCol( pointtab, "DIRECTION" ) ;
718 String dirref = pdirCol.keywordSet().asRecord("MEASINFO").asString("Ref") ;
719 if ( dirref == "AZELGEO" || dirref == "AZEL" ) {
720 dirref = "J2000" ;
721 }
722 sscanf( dirref.chars()+1, "%f", &sdh.equinox ) ;
723 sdh.epoch = "UTC" ;
724 if (sdh.freqref == "TOPO") {
725 sdh.freqref = "TOPOCENT";
726 } else if (sdh.freqref == "GEO") {
727 sdh.freqref = "GEOCENTR";
728 } else if (sdh.freqref == "BARY") {
729 sdh.freqref = "BARYCENT";
730 } else if (sdh.freqref == "GALACTO") {
731 sdh.freqref = "GALACTOC";
732 } else if (sdh.freqref == "LGROUP") {
733 sdh.freqref = "LOCALGRP";
734 } else if (sdh.freqref == "CMB") {
735 sdh.freqref = "CMBDIPOL";
736 } else if (sdh.freqref == "REST") {
737 sdh.freqref = "SOURCE";
738 }
739
740 if ( sdh.fluxunit == "" || sdh.fluxunit == "CNTS" )
741 sdh.fluxunit = "K" ;
742 table_->setHeader( sdh ) ;
743
744 // save path to POINTING table
745 // 2011/07/06 TN
746 // Path to POINTING table in original MS will not be written
747 // if getPt_ is True
748 Path datapath( tablename_ ) ;
749 if ( !getPt_ ) {
750 String pTabName = datapath.absoluteName() + "/POINTING" ;
751 stab.rwKeywordSet().define( "POINTING", pTabName ) ;
752 }
753
754 // for GBT
755 if ( antennaName.contains( "GBT" ) ) {
756 String goTabName = datapath.absoluteName() + "/GBT_GO" ;
757 stab.rwKeywordSet().define( "GBT_GO", goTabName ) ;
758 }
759
760 // for MS created from ASDM
761 //mstable_.keywordSet().print(cout) ;
762 const TableRecord &msKeys = mstable_.keywordSet() ;
763 uInt nfields = msKeys.nfields() ;
764 for ( uInt ifield = 0 ; ifield < nfields ; ifield++ ) {
765 String name = msKeys.name( ifield ) ;
766 //os_ << "name = " << name << LogIO::POST ;
767 if ( name.find( "ASDM" ) != String::npos ) {
768 String asdmpath = msKeys.asTable( ifield ).tableName() ;
769 os_ << "ASDM table: " << asdmpath << LogIO::POST ;
770 stab.rwKeywordSet().define( name, asdmpath ) ;
771 }
772 }
773
774 //double endSec = mathutil::gettimeofday_sec() ;
775 //os_ << "end MSFiller::fill() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
776}
777
778void MSFiller::close()
779{
780 //tablesel_.closeSubTables() ;
781 mstable_.closeSubTables() ;
782 //tablesel_.unlock() ;
783 mstable_.unlock() ;
784}
785
786Int MSFiller::getSrcType( Int stateId, boost::object_pool<ROTableColumn> *tpool )
787{
788 //double startSec = mathutil::gettimeofday_sec() ;
789 //os_ << "start MSFiller::getSrcType() startSec=" << startSec << LogIO::POST ;
790
791 MSState statetab = mstable_.state() ;
792 String obsMode = asString( "OBS_MODE", stateId, statetab, tpool ) ;
793 Bool sig = asBool( "SIG", stateId, statetab, tpool ) ;
794 Bool ref = asBool( "REF", stateId, statetab, tpool ) ;
795 Double cal = asDouble( "CAL", stateId, statetab, tpool ) ;
796 //os_ << "OBS_MODE = " << obsMode << LogIO::POST ;
797
798 // determine separator
799 String sep = "" ;
800 String tmpStr = obsMode.substr( 0, obsMode.find_first_of( "," ) ) ;
801 //os_ << "tmpStr = " << tmpStr << LogIO::POST ;
802 //if ( obsMode.find( ":" ) != String::npos ) {
803 if ( tmpStr.find( ":" ) != String::npos ) {
804 sep = ":" ;
805 }
806 //else if ( obsMode.find( "." ) != String::npos ) {
807 else if ( tmpStr.find( "." ) != String::npos ) {
808 sep = "." ;
809 }
810 else if ( tmpStr.find( "#" ) != String::npos ) {
811 sep = "#" ;
812 }
813 //else if ( obsMode.find( "_" ) != String::npos ) {
814 else if ( tmpStr.find( "_" ) != String::npos ) {
815 sep = "_" ;
816 }
817 //os_ << "separator = " << sep << LogIO::POST ;
818
819 // determine SRCTYPE
820 Int srcType = SrcType::NOTYPE ;
821 if ( sep == ":" ) {
822 // sep == ":"
823 //
824 // GBT case
825 //
826 // obsMode1=Nod
827 // NOD
828 // obsMode1=OffOn
829 // obsMode2=PSWITCHON: PSON
830 // obsMode2=PSWITCHOFF: PSOFF
831 // obsMode1=??
832 // obsMode2=FSWITCH:
833 // SIG=1: FSON
834 // REF=1: FSOFF
835 // Calibration scan if CAL != 0
836 Int epos = obsMode.find_first_of( sep ) ;
837 Int nextpos = obsMode.find_first_of( sep, epos+1 ) ;
838 String obsMode1 = obsMode.substr( 0, epos ) ;
839 String obsMode2 = obsMode.substr( epos+1, nextpos-epos-1 ) ;
840 if ( obsMode1 == "Nod" ) {
841 srcType = SrcType::NOD ;
842 }
843 else if ( obsMode1 == "OffOn" ) {
844 if ( obsMode2 == "PSWITCHON" ) srcType = SrcType::PSON ;
845 if ( obsMode2 == "PSWITCHOFF" ) srcType = SrcType::PSOFF ;
846 }
847 else {
848 if ( obsMode2 == "FSWITCH" ) {
849 if ( sig ) srcType = SrcType::FSON ;
850 if ( ref ) srcType = SrcType::FSOFF ;
851 }
852 }
853 if ( cal > 0.0 ) {
854 if ( srcType == SrcType::NOD )
855 srcType = SrcType::NODCAL ;
856 else if ( srcType == SrcType::PSON )
857 srcType = SrcType::PONCAL ;
858 else if ( srcType == SrcType::PSOFF )
859 srcType = SrcType::POFFCAL ;
860 else if ( srcType == SrcType::FSON )
861 srcType = SrcType::FONCAL ;
862 else if ( srcType == SrcType::FSOFF )
863 srcType = SrcType::FOFFCAL ;
864 else
865 srcType = SrcType::CAL ;
866 }
867 }
868 else if ( sep == "." || sep == "#" ) {
869 // sep == "." or "#"
870 //
871 // ALMA & EVLA case (MS via ASDM) before3.1
872 //
873 // obsMode1=CALIBRATE_*
874 // obsMode2=ON_SOURCE: PONCAL
875 // obsMode2=OFF_SOURCE: POFFCAL
876 // obsMode1=OBSERVE_TARGET
877 // obsMode2=ON_SOURCE: PON
878 // obsMode2=OFF_SOURCE: POFF
879 string substr[2] ;
880 int numSubstr = split( obsMode, substr, 2, "," ) ;
881 //os_ << "numSubstr = " << numSubstr << LogIO::POST ;
882 //for ( int i = 0 ; i < numSubstr ; i++ )
883 //os_ << "substr[" << i << "] = " << substr[i] << LogIO::POST ;
884 String obsType( substr[0] ) ;
885 //os_ << "obsType = " << obsType << LogIO::POST ;
886 Int epos = obsType.find_first_of( sep ) ;
887 Int nextpos = obsType.find_first_of( sep, epos+1 ) ;
888 String obsMode1 = obsType.substr( 0, epos ) ;
889 String obsMode2 = obsType.substr( epos+1, nextpos-epos-1 ) ;
890 //os_ << "obsMode1 = " << obsMode1 << LogIO::POST ;
891 //os_ << "obsMode2 = " << obsMode2 << LogIO::POST ;
892 if ( obsMode1.find( "CALIBRATE_" ) == 0 ) {
893 if ( obsMode2 == "ON_SOURCE" ) srcType = SrcType::PONCAL ;
894 if ( obsMode2 == "OFF_SOURCE" ) srcType = SrcType::POFFCAL ;
895 }
896 else if ( obsMode1 == "OBSERVE_TARGET" ) {
897 if ( obsMode2 == "ON_SOURCE" ) srcType = SrcType::PSON ;
898 if ( obsMode2 == "OFF_SOURCE" ) srcType = SrcType::PSOFF ;
899 }
900 }
901 else if ( sep == "_" ) {
902 // sep == "_"
903 //
904 // ALMA & EVLA case (MS via ASDM) after 3.2
905 //
906 // obsMode1=CALIBRATE_*
907 // obsMode2=ON_SOURCE: PONCAL
908 // obsMode2=OFF_SOURCE: POFFCAL
909 // obsMode1=OBSERVE_TARGET
910 // obsMode2=ON_SOURCE: PON
911 // obsMode2=OFF_SOURCE: POFF
912 string substr[2] ;
913 int numSubstr = split( obsMode, substr, 2, "," ) ;
914 //os_ << "numSubstr = " << numSubstr << LogIO::POST ;
915 //for ( int i = 0 ; i < numSubstr ; i++ )
916 //os_ << "substr[" << i << "] = " << substr[i] << LogIO::POST ;
917 String obsType( substr[0] ) ;
918 //os_ << "obsType = " << obsType << LogIO::POST ;
919 string substr2[4] ;
920 int numSubstr2 = split( obsType, substr2, 4, sep ) ;
921 //Int epos = obsType.find_first_of( sep ) ;
922 //Int nextpos = obsType.find_first_of( sep, epos+1 ) ;
923 //String obsMode1 = obsType.substr( 0, epos ) ;
924 //String obsMode2 = obsType.substr( epos+1, nextpos-epos-1 ) ;
925 String obsMode1( substr2[0] ) ;
926 String obsMode2( substr2[2] ) ;
927 //os_ << "obsMode1 = " << obsMode1 << LogIO::POST ;
928 //os_ << "obsMode2 = " << obsMode2 << LogIO::POST ;
929 if ( obsMode1.find( "CALIBRATE" ) == 0 ) {
930 if ( obsMode2 == "ON" ) srcType = SrcType::PONCAL ;
931 if ( obsMode2 == "OFF" ) srcType = SrcType::POFFCAL ;
932 }
933 else if ( obsMode1 == "OBSERVE" ) {
934 if ( obsMode2 == "ON" ) srcType = SrcType::PSON ;
935 if ( obsMode2 == "OFF" ) srcType = SrcType::PSOFF ;
936 }
937 }
938 else {
939 if ( sig ) srcType = SrcType::SIG ;
940 if ( ref ) srcType = SrcType::REF ;
941 }
942
943 //os_ << "srcType = " << srcType << LogIO::POST ;
944 //double endSec = mathutil::gettimeofday_sec() ;
945 //os_ << "end MSFiller::getSrcType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
946 return srcType ;
947}
948
949//Vector<uInt> MSFiller::getPolNo( Int corrType )
950Block<uInt> MSFiller::getPolNo( Int corrType )
951{
952 //double startSec = mathutil::gettimeofday_sec() ;
953 //os_ << "start MSFiller::getPolNo() startSec=" << startSec << LogIO::POST ;
954 Block<uInt> polno( 1 ) ;
955
956 if ( corrType == Stokes::I || corrType == Stokes::RR || corrType == Stokes::XX ) {
957 polno = 0 ;
958 }
959 else if ( corrType == Stokes::Q || corrType == Stokes::LL || corrType == Stokes::YY ) {
960 polno = 1 ;
961 }
962 else if ( corrType == Stokes::U ) {
963 polno = 2 ;
964 }
965 else if ( corrType == Stokes::V ) {
966 polno = 3 ;
967 }
968 else if ( corrType == Stokes::RL || corrType == Stokes::XY || corrType == Stokes::LR || corrType == Stokes::RL ) {
969 polno.resize( 2 ) ;
970 polno[0] = 2 ;
971 polno[1] = 3 ;
972 }
973 else if ( corrType == Stokes::Plinear ) {
974 polno[0] = 1 ;
975 }
976 else if ( corrType == Stokes::Pangle ) {
977 polno[0] = 2 ;
978 }
979 else {
980 polno = 99 ;
981 }
982 //os_ << "polno = " << polno << LogIO::POST ;
983 //double endSec = mathutil::gettimeofday_sec() ;
984 //os_ << "end MSFiller::getPolNo() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
985
986 return polno ;
987}
988
989String MSFiller::getPolType( Int corrType )
990{
991 //double startSec = mathutil::gettimeofday_sec() ;
992 //os_ << "start MSFiller::getPolType() startSec=" << startSec << LogIO::POST ;
993 String poltype = "" ;
994
995 if ( corrType == Stokes::I || corrType == Stokes::Q || corrType == Stokes::U || corrType == Stokes::V )
996 poltype = "stokes" ;
997 else if ( corrType == Stokes::XX || corrType == Stokes::YY || corrType == Stokes::XY || corrType == Stokes::YX )
998 poltype = "linear" ;
999 else if ( corrType == Stokes::RR || corrType == Stokes::LL || corrType == Stokes::RL || corrType == Stokes::LR )
1000 poltype = "circular" ;
1001 else if ( corrType == Stokes::Plinear || corrType == Stokes::Pangle )
1002 poltype = "linpol" ;
1003
1004 //double endSec = mathutil::gettimeofday_sec() ;
1005 //os_ << "end MSFiller::getPolType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1006 return poltype ;
1007}
1008
1009void MSFiller::fillWeather()
1010{
1011 //double startSec = mathutil::gettimeofday_sec() ;
1012 //os_ << "start MSFiller::fillWeather() startSec=" << startSec << LogIO::POST ;
1013
1014 if ( !isWeather_ ) {
1015 // add dummy row
1016 table_->weather().table().addRow(1,True) ;
1017 return ;
1018 }
1019
1020 Table mWeather = mstable_.weather() ;
1021 //Table mWeatherSel = mWeather( mWeather.col("ANTENNA_ID") == antenna_ ).sort("TIME") ;
1022 Table mWeatherSel( mWeather( mWeather.col("ANTENNA_ID") == antenna_ ).sort("TIME") ) ;
1023 //os_ << "mWeatherSel.nrow() = " << mWeatherSel.nrow() << LogIO::POST ;
1024 if ( mWeatherSel.nrow() == 0 ) {
1025 os_ << "No rows with ANTENNA_ID = " << antenna_ << " in WEATHER table, Try -1..." << LogIO::POST ;
1026 mWeatherSel = Table( MSWeather( mWeather( mWeather.col("ANTENNA_ID") == -1 ) ) ) ;
1027 if ( mWeatherSel.nrow() == 0 ) {
1028 os_ << "No rows in WEATHER table" << LogIO::POST ;
1029 }
1030 }
1031 uInt wnrow = mWeatherSel.nrow() ;
1032 //os_ << "wnrow = " << wnrow << LogIO::POST ;
1033
1034 if ( wnrow == 0 )
1035 return ;
1036
1037 Table wtab = table_->weather().table() ;
1038 wtab.addRow( wnrow ) ;
1039
1040 Bool stationInfoExists = mWeatherSel.tableDesc().isColumn( "NS_WX_STATION_ID" ) ;
1041 Int stationId = -1 ;
1042 if ( stationInfoExists ) {
1043 // determine which station is closer
1044 ROScalarColumn<Int> stationCol( mWeatherSel, "NS_WX_STATION_ID" ) ;
1045 ROArrayColumn<Double> stationPosCol( mWeatherSel, "NS_WX_STATION_POSITION" ) ;
1046 Vector<Int> stationIds = stationCol.getColumn() ;
1047 Vector<Int> stationIdList( 0 ) ;
1048 Matrix<Double> stationPosList( 0, 3, 0.0 ) ;
1049 uInt numStation = 0 ;
1050 for ( uInt i = 0 ; i < stationIds.size() ; i++ ) {
1051 if ( !anyEQ( stationIdList, stationIds[i] ) ) {
1052 numStation++ ;
1053 stationIdList.resize( numStation, True ) ;
1054 stationIdList[numStation-1] = stationIds[i] ;
1055 stationPosList.resize( numStation, 3, True ) ;
1056 stationPosList.row( numStation-1 ) = stationPosCol( i ) ;
1057 }
1058 }
1059 //os_ << "staionIdList = " << stationIdList << endl ;
1060 Table mAntenna = mstable_.antenna() ;
1061 ROArrayColumn<Double> antposCol( mAntenna, "POSITION" ) ;
1062 Vector<Double> antpos = antposCol( antenna_ ) ;
1063 Double minDiff = -1.0 ;
1064 for ( uInt i = 0 ; i < stationIdList.size() ; i++ ) {
1065 Double diff = sum( square( antpos - stationPosList.row( i ) ) ) ;
1066 if ( minDiff < 0.0 || minDiff > diff ) {
1067 minDiff = diff ;
1068 stationId = stationIdList[i] ;
1069 }
1070 }
1071 }
1072 //os_ << "stationId = " << stationId << endl ;
1073
1074 ScalarColumn<Float> *fCol ;
1075 ROScalarColumn<Float> *sharedFloatCol ;
1076 if ( mWeatherSel.tableDesc().isColumn( "TEMPERATURE" ) ) {
1077 fCol = new ScalarColumn<Float>( wtab, "TEMPERATURE" ) ;
1078 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "TEMPERATURE" ) ;
1079 fCol->putColumn( *sharedFloatCol ) ;
1080 delete sharedFloatCol ;
1081 delete fCol ;
1082 }
1083 if ( mWeatherSel.tableDesc().isColumn( "PRESSURE" ) ) {
1084 fCol = new ScalarColumn<Float>( wtab, "PRESSURE" ) ;
1085 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "PRESSURE" ) ;
1086 fCol->putColumn( *sharedFloatCol ) ;
1087 delete sharedFloatCol ;
1088 delete fCol ;
1089 }
1090 if ( mWeatherSel.tableDesc().isColumn( "REL_HUMIDITY" ) ) {
1091 fCol = new ScalarColumn<Float>( wtab, "HUMIDITY" ) ;
1092 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "REL_HUMIDITY" ) ;
1093 fCol->putColumn( *sharedFloatCol ) ;
1094 delete sharedFloatCol ;
1095 delete fCol ;
1096 }
1097 if ( mWeatherSel.tableDesc().isColumn( "WIND_SPEED" ) ) {
1098 fCol = new ScalarColumn<Float>( wtab, "WINDSPEED" ) ;
1099 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "WIND_SPEED" ) ;
1100 fCol->putColumn( *sharedFloatCol ) ;
1101 delete sharedFloatCol ;
1102 delete fCol ;
1103 }
1104 if ( mWeatherSel.tableDesc().isColumn( "WIND_DIRECTION" ) ) {
1105 fCol = new ScalarColumn<Float>( wtab, "WINDAZ" ) ;
1106 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "WIND_DIRECTION" ) ;
1107 fCol->putColumn( *sharedFloatCol ) ;
1108 delete sharedFloatCol ;
1109 delete fCol ;
1110 }
1111 ScalarColumn<uInt> idCol( wtab, "ID" ) ;
1112 for ( uInt irow = 0 ; irow < wnrow ; irow++ )
1113 idCol.put( irow, irow ) ;
1114
1115 ROScalarQuantColumn<Double> tqCol( mWeatherSel, "TIME" ) ;
1116 ROScalarColumn<Double> tCol( mWeatherSel, "TIME" ) ;
1117 String tUnit = tqCol.getUnits() ;
1118 Vector<Double> mwTime = tCol.getColumn() ;
1119 if ( tUnit == "d" )
1120 mwTime *= 86400.0 ;
1121 tqCol.attach( mWeatherSel, "INTERVAL" ) ;
1122 tCol.attach( mWeatherSel, "INTERVAL" ) ;
1123 String iUnit = tqCol.getUnits() ;
1124 Vector<Double> mwInterval = tCol.getColumn() ;
1125 if ( iUnit == "d" )
1126 mwInterval *= 86400.0 ;
1127
1128 if ( stationId > 0 ) {
1129 ROScalarColumn<Int> stationCol( mWeatherSel, "NS_WX_STATION_ID" ) ;
1130 Vector<Int> stationVec = stationCol.getColumn() ;
1131 uInt wsnrow = ntrue( stationVec == stationId ) ;
1132 mwTime_.resize( wsnrow ) ;
1133 mwInterval_.resize( wsnrow ) ;
1134 mwIndex_.resize( wsnrow ) ;
1135 uInt wsidx = 0 ;
1136 for ( uInt irow = 0 ; irow < wnrow ; irow++ ) {
1137 if ( stationId == stationVec[irow] ) {
1138 mwTime_[wsidx] = mwTime[irow] ;
1139 mwInterval_[wsidx] = mwInterval[irow] ;
1140 mwIndex_[wsidx] = irow ;
1141 wsidx++ ;
1142 }
1143 }
1144 }
1145 else {
1146 mwTime_ = mwTime ;
1147 mwInterval_ = mwInterval ;
1148 mwIndex_.resize( mwTime_.size() ) ;
1149 indgen( mwIndex_ ) ;
1150 }
1151 //os_ << "mwTime[0] = " << mwTime_[0] << " mwInterval[0] = " << mwInterval_[0] << LogIO::POST ;
1152 //double endSec = mathutil::gettimeofday_sec() ;
1153 //os_ << "end MSFiller::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1154}
1155
1156void MSFiller::fillFocus()
1157{
1158 //double startSec = mathutil::gettimeofday_sec() ;
1159 //os_ << "start MSFiller::fillFocus() startSec=" << startSec << LogIO::POST ;
1160 // tentative
1161 table_->focus().addEntry( 0.0, 0.0, 0.0, 0.0 ) ;
1162 //double endSec = mathutil::gettimeofday_sec() ;
1163 //os_ << "end MSFiller::fillFocus() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1164}
1165
1166void MSFiller::fillTcal( boost::object_pool<ROTableColumn> *tpoolr )
1167{
1168 //double startSec = mathutil::gettimeofday_sec() ;
1169 //os_ << "start MSFiller::fillTcal() startSec=" << startSec << LogIO::POST ;
1170
1171 if ( !isSysCal_ ) {
1172 // add dummy row
1173 os_ << "No SYSCAL rows" << LogIO::POST ;
1174 table_->tcal().table().addRow(1,True) ;
1175 Vector<Float> defaultTcal( 1, 1.0 ) ;
1176 ArrayColumn<Float> tcalCol( table_->tcal().table(), "TCAL" ) ;
1177 tcalCol.put( 0, defaultTcal ) ;
1178 return ;
1179 }
1180
1181 if ( colTcal_ == "NONE" ) {
1182 // add dummy row
1183 os_ << "No TCAL column" << LogIO::POST ;
1184 table_->tcal().table().addRow(1,True) ;
1185 Vector<Float> defaultTcal( 1, 1.0 ) ;
1186 ArrayColumn<Float> tcalCol( table_->tcal().table(), "TCAL" ) ;
1187 tcalCol.put( 0, defaultTcal ) ;
1188 return ;
1189 }
1190
1191 Table sctab = mstable_.sysCal() ;
1192 if ( sctab.nrow() == 0 ) {
1193 os_ << "No SYSCAL rows" << LogIO::POST ;
1194 return ;
1195 }
1196 Table sctabsel( sctab( sctab.col("ANTENNA_ID") == antenna_ ) ) ;
1197 if ( sctabsel.nrow() == 0 ) {
1198 os_ << "No SYSCAL rows" << LogIO::POST ;
1199 return ;
1200 }
1201 ROArrayColumn<Float> *tmpTcalCol = new ROArrayColumn<Float>( sctabsel, colTcal_ ) ;
1202 // return if any rows without Tcal value exists
1203 Bool notDefined = False ;
1204 for ( uInt irow = 0 ; irow < sctabsel.nrow() ; irow++ ) {
1205 if ( !tmpTcalCol->isDefined( irow ) ) {
1206 notDefined = True ;
1207 break ;
1208 }
1209 }
1210 if ( notDefined ) {
1211 os_ << "No TCAL value" << LogIO::POST ;
1212 delete tmpTcalCol ;
1213 table_->tcal().table().addRow(1,True) ;
1214 Vector<Float> defaultTcal( 1, 1.0 ) ;
1215 ArrayColumn<Float> tcalCol( table_->tcal().table(), "TCAL" ) ;
1216 tcalCol.put( 0, defaultTcal ) ;
1217 return ;
1218 }
1219 uInt npol = tmpTcalCol->shape( 0 )(0) ;
1220 delete tmpTcalCol ;
1221 //os_ << "fillTcal(): npol = " << npol << LogIO::POST ;
1222 Table tab = table_->tcal().table() ;
1223 ArrayColumn<Float> tcalCol( tab, "TCAL" ) ;
1224 uInt oldnr = 0 ;
1225 uInt newnr = 0 ;
1226 TableRow row( tab ) ;
1227 TableRecord &trec = row.record() ;
1228 RecordFieldPtr<uInt> idRF( trec, "ID" ) ;
1229 RecordFieldPtr<String> timeRF( trec, "TIME" ) ;
1230 RecordFieldPtr< Array<Float> > tcalRF( trec, "TCAL" ) ;
1231 TableIterator iter0( sctabsel, "FEED_ID" ) ;
1232 while( !iter0.pastEnd() ) {
1233 Table t0 = iter0.table() ;
1234 Int feedId = asInt( "FEED_ID", 0, t0, tpoolr ) ;
1235 TableIterator iter1( t0, "SPECTRAL_WINDOW_ID" ) ;
1236 while( !iter1.pastEnd() ) {
1237 Table t1 = iter1.table() ;
1238 Int spwId = asInt( "SPECTRAL_WINDOW_ID", 0, t1, tpoolr ) ;
1239 tmpTcalCol = new ROArrayColumn<Float>( t1, colTcal_ ) ;
1240 ROScalarQuantColumn<Double> scTimeCol( t1, "TIME" ) ;
1241 Vector<uInt> idminmax( 2, oldnr ) ;
1242 for ( uInt irow = 0 ; irow < t1.nrow() ; irow++ ) {
1243 String sTime = MVTime( scTimeCol(irow) ).string( MVTime::YMD ) ;
1244 *timeRF = sTime ;
1245 uInt idx = oldnr ;
1246 Matrix<Float> subtcal = (*tmpTcalCol)( irow ) ;
1247 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
1248 *idRF = idx++ ;
1249 //*tcalRF = subtcal.row( ipol ) ;
1250 tcalRF.define( subtcal.row( ipol ) ) ;
1251
1252 // commit row
1253 tab.addRow() ;
1254 row.put( tab.nrow()-1 ) ;
1255
1256 newnr++ ;
1257 }
1258 idminmax[0] = oldnr ;
1259 idminmax[1] = newnr - 1 ;
1260 oldnr = newnr ;
1261
1262 String key = keyTcal( feedId, spwId, sTime ) ;
1263 tcalrec_.define( key, idminmax ) ;
1264 }
1265 delete tmpTcalCol ;
1266 iter1++ ;
1267 }
1268 iter0++ ;
1269 }
1270
1271 //tcalrec_.print( std::cout ) ;
1272 //double endSec = mathutil::gettimeofday_sec() ;
1273 //os_ << "end MSFiller::fillTcal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1274}
1275
1276uInt MSFiller::getWeatherId( uInt idx, Double wtime )
1277{
1278 //double startSec = mathutil::gettimeofday_sec() ;
1279 //os_ << "start MSFiller::getWeatherId() startSec=" << startSec << LogIO::POST ;
1280 uInt nrow = mwTime_.size() ;
1281 if ( nrow < 2 )
1282 return 0 ;
1283 uInt wid = nrow ;
1284 if ( idx == 0 ) {
1285 wid = 0 ;
1286 Double tStart = mwTime_[wid]-0.5*mwInterval_[wid] ;
1287 if ( wtime < tStart )
1288 return wid ;
1289 }
1290 for ( uInt i = idx ; i < nrow-1 ; i++ ) {
1291 Double tStart = mwTime_[i]-0.5*mwInterval_[i] ;
1292 // use of INTERVAL column is problematic
1293 // since there are "blank" time of weather monitoring
1294 //Double tEnd = tStart + mwInterval_[i] ;
1295 Double tEnd = mwTime_[i+1]-0.5*mwInterval_[i+1] ;
1296 //os_ << "tStart = " << tStart << " dtEnd = " << tEnd-tStart << " dwtime = " << wtime-tStart << LogIO::POST ;
1297 if ( wtime >= tStart && wtime <= tEnd ) {
1298 wid = i ;
1299 break ;
1300 }
1301 }
1302 if ( wid == nrow ) {
1303 uInt i = nrow - 1 ;
1304 Double tStart = mwTime_[i-1]+0.5*mwInterval_[i-1] ;
1305 Double tEnd = mwTime_[i]+0.5*mwInterval_[i] ;
1306 //os_ << "tStart = " << tStart << " dtEnd = " << tEnd-tStart << " dwtime = " << wtime-tStart << LogIO::POST ;
1307 if ( wtime >= tStart && wtime <= tEnd )
1308 wid = i-1 ;
1309 else
1310 wid = i ;
1311 }
1312
1313 //if ( wid == nrow )
1314 //os_ << LogIO::WARN << "Couldn't find correct WEATHER_ID for time " << wtime << LogIO::POST ;
1315
1316 //double endSec = mathutil::gettimeofday_sec() ;
1317 //os_ << "end MSFiller::getWeatherId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1318 return wid ;
1319}
1320
1321void MSFiller::getSysCalTime( Vector<MEpoch> &scTime, Vector<Double> &scInterval, Block<MEpoch> &tcol, Block<Int> &tidx )
1322{
1323 //double startSec = mathutil::gettimeofday_sec() ;
1324 //os_ << "start MSFiller::getSysCalTime() startSec=" << startSec << LogIO::POST ;
1325
1326 if ( !isSysCal_ )
1327 return ;
1328
1329 uInt nrow = tidx.nelements() ;
1330 if ( scTime.nelements() == 0 )
1331 return ;
1332 else if ( scTime.nelements() == 1 ) {
1333 tidx[0] = 0 ;
1334 return ;
1335 }
1336 uInt scnrow = scTime.nelements() ;
1337 uInt idx = 0 ;
1338 const Double half = 0.5e0 ;
1339 // execute binary search
1340 idx = binarySearch( scTime, tcol[0].get( "s" ).getValue() ) ;
1341 if ( idx != 0 )
1342 idx -= 1 ;
1343 for ( uInt i = 0 ; i < nrow ; i++ ) {
1344 Double t = tcol[i].get( "s" ).getValue() ;
1345 Double tsc = scTime[0].get( "s" ).getValue() ;
1346 if ( t < tsc ) {
1347 tidx[i] = 0 ;
1348 continue ;
1349 }
1350 for ( uInt j = idx ; j < scnrow-1 ; j++ ) {
1351 Double tsc1 = scTime[j].get( "s" ).getValue() ;
1352 Double dt1 = scInterval[j] ;
1353 Double tsc2 = scTime[j+1].get( "s" ).getValue() ;
1354 Double dt2 = scInterval[j+1] ;
1355 if ( t > tsc1-half*dt1 && t <= tsc2-half*dt2 ) {
1356 tidx[i] = j ;
1357 idx = j ;
1358 break ;
1359 }
1360 }
1361 if ( tidx[i] == -1 ) {
1362// Double tsc = scTime[scnrow-1].get( "s" ).getValue() ;
1363// Double dt = scInterval[scnrow-1] ;
1364// if ( t <= tsc+0.5*dt ) {
1365// tidx[i] = scnrow-1 ;
1366// }
1367 tidx[i] = scnrow-1 ;
1368 }
1369 }
1370 //double endSec = mathutil::gettimeofday_sec() ;
1371 //os_ << "end MSFiller::getSysCalTime() endSec=" << endSec << " (" << endSec-startSec << "sec) scnrow = " << scnrow << " tcol.nelements = " << tcol.nelements() << LogIO::POST ;
1372 return ;
1373}
1374
1375Block<uInt> MSFiller::getTcalId( Int fid, Int spwid, MEpoch &t )
1376{
1377 //double startSec = mathutil::gettimeofday_sec() ;
1378 //os_ << "start MSFiller::getTcalId() startSec=" << startSec << LogIO::POST ;
1379 //if ( table_->tcal().table().nrow() == 0 ) {
1380 if ( !isSysCal_ ) {
1381 os_ << "No TCAL rows" << LogIO::POST ;
1382 Block<uInt> tcalids( 4, 0 ) ;
1383 return tcalids ;
1384 }
1385 //String sctime = MVTime( Quantum<Double>(t,"s") ).string(MVTime::YMD) ;
1386 String sctime = MVTime( t.getValue() ).string(MVTime::YMD) ;
1387 String key = keyTcal( fid, spwid, sctime ) ;
1388 if ( !tcalrec_.isDefined( key ) ) {
1389 os_ << "No TCAL rows" << LogIO::POST ;
1390 Block<uInt> tcalids( 4, 0 ) ;
1391 return tcalids ;
1392 }
1393 Vector<uInt> ids = tcalrec_.asArrayuInt( key ) ;
1394 uInt npol = ids[1] - ids[0] + 1 ;
1395 Block<uInt> tcalids( npol ) ;
1396 tcalids[0] = ids[0] ;
1397 tcalids[1] = ids[1] ;
1398 for ( uInt ipol = 2 ; ipol < npol ; ipol++ )
1399 tcalids[ipol] = ids[0] + ipol - 1 ;
1400
1401 //double endSec = mathutil::gettimeofday_sec() ;
1402 //os_ << "end MSFiller::getTcalId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1403 return tcalids ;
1404}
1405
1406uInt MSFiller::getDirection( uInt idx,
1407 Vector<Double> &dir,
1408 Vector<Double> &srate,
1409 String &ref,
1410 Vector<Double> &tcol,
1411 ROArrayColumn<Double> &dcol,
1412 Double t )
1413{
1414 //double startSec = mathutil::gettimeofday_sec() ;
1415 //os_ << "start MSFiller::getDirection1() startSec=" << startSec << LogIO::POST ;
1416 //double time0 = mathutil::gettimeofday_sec() ;
1417 //os_ << "start getDirection 1st stage startSec=" << time0 << LogIO::POST ;
1418 // assume that cols is sorted by TIME
1419 Bool doInterp = False ;
1420 uInt nrow = dcol.nrow() ;
1421 if ( nrow == 0 )
1422 return 0 ;
1423 if ( idx == 0 ) {
1424 uInt nrowb = 1000 ;
1425 if ( nrow > nrowb ) {
1426 uInt nblock = nrow / nrowb + 1 ;
1427 for ( uInt iblock = 0 ; iblock < nblock ; iblock++ ) {
1428 uInt high = min( nrowb, nrow-iblock*nrowb ) ;
1429
1430 if ( tcol( high-1 ) < t ) {
1431 idx = iblock * nrowb ;
1432 continue ;
1433 }
1434
1435 Slice slice( iblock*nrowb, nrowb ) ;
1436 Vector<Double> tarr = tcol( slice ) ;
1437
1438 uInt bidx = binarySearch( tarr, t ) ;
1439
1440 idx = iblock * nrowb + bidx ;
1441 break ;
1442 }
1443 }
1444 else {
1445 idx = binarySearch( tcol, t ) ;
1446 }
1447 }
1448 //double time1 = mathutil::gettimeofday_sec() ;
1449 //os_ << "end getDirection 1st stage endSec=" << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
1450 // ensure that tcol(idx) < t
1451 //os_ << "tcol(idx) = " << tcol(idx).get("s").getValue() << " t = " << t << " diff = " << tcol(idx).get("s").getValue()-t << endl ;
1452 //time0 = mathutil::gettimeofday_sec() ;
1453 //os_ << "start getDirection 2nd stage startSec=" << time0 << LogIO::POST ;
1454 //while( tcol( idx ) * factor > t && idx > 0 )
1455 while( tcol[idx] > t && idx > 0 )
1456 idx-- ;
1457 //os_ << "idx = " << idx << LogIO::POST ;
1458
1459 // index search
1460 for ( uInt i = idx ; i < nrow ; i++ ) {
1461 Double tref = tcol[i] ;
1462 if ( tref == t ) {
1463 idx = i ;
1464 break ;
1465 }
1466 else if ( tref > t ) {
1467 if ( i == 0 ) {
1468 idx = i ;
1469 }
1470 else {
1471 idx = i-1 ;
1472 doInterp = True ;
1473 }
1474 break ;
1475 }
1476 else {
1477 idx = nrow - 1 ;
1478 }
1479 }
1480 //time1 = mathutil::gettimeofday_sec() ;
1481 //os_ << "end getDirection 2nd stage endSec=" << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
1482 //os_ << "searched idx = " << idx << LogIO::POST ;
1483
1484 //time0 = mathutil::gettimeofday_sec() ;
1485 //os_ << "start getDirection 3rd stage startSec=" << time0 << LogIO::POST ;
1486 //os_ << "dmcol(idx).shape() = " << dmcol(idx).shape() << LogIO::POST ;
1487 //IPosition ip( dmcol(idx).shape().nelements(), 0 ) ;
1488 IPosition ip( dcol(idx).shape().nelements(), 0 ) ;
1489 //os_ << "ip = " << ip << LogIO::POST ;
1490 //ref = dmcol(idx)(ip).getRefString() ;
1491 TableRecord trec = dcol.keywordSet() ;
1492 Record rec = trec.asRecord( "MEASINFO" ) ;
1493 ref = rec.asString( "Ref" ) ;
1494 //os_ << "ref = " << ref << LogIO::POST ;
1495 if ( doInterp ) {
1496 //os_ << "do interpolation" << LogIO::POST ;
1497 //os_ << "dcol(idx).shape() = " << dcol(idx).shape() << LogIO::POST ;
1498 Double tref0 = tcol[idx] ;
1499 Double tref1 = tcol[idx+1] ;
1500 Matrix<Double> mdir0 = dcol( idx ) ;
1501 Matrix<Double> mdir1 = dcol( idx+1 ) ;
1502 Vector<Double> dir0 = mdir0.column( 0 ) ;
1503 //os_ << "dir0 = " << dir0 << LogIO::POST ;
1504 Vector<Double> dir1 = mdir1.column( 0 ) ;
1505 //os_ << "dir1 = " << dir1 << LogIO::POST ;
1506 Double dt0 = t - tref0 ;
1507 Double dt1 = tref1 - t ;
1508 dir.reference( (dt0*dir1+dt1*dir0)/(dt0+dt1) ) ;
1509 if ( mdir0.ncolumn() > 1 ) {
1510 if ( dt0 >= dt1 )
1511 srate.reference( mdir0.column( 1 ) ) ;
1512 else
1513 srate.reference( mdir1.column( 1 ) ) ;
1514 }
1515 //os_ << "dir = " << dir << LogIO::POST ;
1516 }
1517 else {
1518 //os_ << "no interpolation" << LogIO::POST ;
1519 Matrix<Double> mdir0 = dcol( idx ) ;
1520 dir.reference( mdir0.column( 0 ) ) ;
1521 if ( mdir0.ncolumn() > 1 )
1522 srate.reference( mdir0.column( 1 ) ) ;
1523 }
1524
1525 //time1 = mathutil::gettimeofday_sec() ;
1526 //os_ << "end getDirection 3rd stage endSec=" << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
1527 //double endSec = mathutil::gettimeofday_sec() ;
1528 //os_ << "end MSFiller::getDirection1() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1529 return idx ;
1530}
1531
1532String MSFiller::keyTcal( Int feedid, Int spwid, String stime )
1533{
1534 String sfeed = "FEED" + String::toString( feedid ) ;
1535 String sspw = "SPW" + String::toString( spwid ) ;
1536 return sfeed+":"+sspw+":"+stime ;
1537}
1538
1539uInt MSFiller::binarySearch( Vector<MEpoch> &timeList, Double target )
1540{
1541 Int low = 0 ;
1542 Int high = timeList.nelements() ;
1543 uInt idx = 0 ;
1544
1545 while ( low <= high ) {
1546 idx = (Int)( 0.5 * ( low + high ) ) ;
1547 Double t = timeList[idx].get( "s" ).getValue() ;
1548 if ( t < target )
1549 low = idx + 1 ;
1550 else if ( t > target )
1551 high = idx - 1 ;
1552 else {
1553 return idx ;
1554 }
1555 }
1556
1557 idx = max( 0, min( low, high ) ) ;
1558
1559 return idx ;
1560}
1561
1562uInt MSFiller::binarySearch( Vector<Double> &timeList, Double target )
1563{
1564 Int low = 0 ;
1565 Int high = timeList.nelements() ;
1566 uInt idx = 0 ;
1567
1568 while ( low <= high ) {
1569 idx = (Int)( 0.5 * ( low + high ) ) ;
1570 Double t = timeList[idx] ;
1571 if ( t < target )
1572 low = idx + 1 ;
1573 else if ( t > target )
1574 high = idx - 1 ;
1575 else {
1576 return idx ;
1577 }
1578 }
1579
1580 idx = max( 0, min( low, high ) ) ;
1581
1582 return idx ;
1583}
1584
1585string MSFiller::getFrame()
1586{
1587 MFrequency::Types frame = MFrequency::DEFAULT ;
1588 ROTableColumn numChanCol( mstable_.spectralWindow(), "NUM_CHAN" ) ;
1589 ROTableColumn measFreqRefCol( mstable_.spectralWindow(), "MEAS_FREQ_REF" ) ;
1590 uInt nrow = numChanCol.nrow() ;
1591 Vector<Int> measFreqRef( nrow, MFrequency::DEFAULT ) ;
1592 uInt nref = 0 ;
1593 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
1594 if ( numChanCol.asInt( irow ) != 4 ) { // exclude WVR
1595 measFreqRef[nref] = measFreqRefCol.asInt( irow ) ;
1596 nref++ ;
1597 }
1598 }
1599 if ( nref > 0 )
1600 frame = (MFrequency::Types)measFreqRef[0] ;
1601
1602 return MFrequency::showType( frame ) ;
1603}
1604
1605void MSFiller::reshapeSpectraAndFlagtra( Cube<Float> &sp,
1606 Cube<Bool> &fl,
1607 Table &tab,
1608 Int &npol,
1609 Int &nchan,
1610 Int &nrow,
1611 Vector<Int> &corrtype )
1612{
1613 //double startSec = mathutil::gettimeofday_sec() ;
1614 //os_ << "start MSFiller::reshapeSpectraAndFlagtra() startSec=" << startSec << LogIO::POST ;
1615 if ( isFloatData_ ) {
1616 ROArrayColumn<Bool> mFlagCol( tab, "FLAG" ) ;
1617 ROArrayColumn<Float> mFloatDataCol( tab, "FLOAT_DATA" ) ;
1618 mFloatDataCol.getColumn( sp ) ;
1619 mFlagCol.getColumn( fl ) ;
1620 }
1621 else if ( isData_ ) {
1622 sp.resize( npol, nchan, nrow ) ;
1623 fl.resize( npol, nchan, nrow ) ;
1624 ROArrayColumn<Bool> mFlagCol( tab, "FLAG" ) ;
1625 ROArrayColumn<Complex> mDataCol( tab, "DATA" ) ;
1626 if ( npol < 3 ) {
1627 Cube<Float> tmp = ComplexToReal( mDataCol.getColumn() ) ;
1628 IPosition start( 3, 0, 0, 0 ) ;
1629 IPosition end( 3, 2*npol-1, nchan-1, nrow-1 ) ;
1630 IPosition inc( 3, 2, 1, 1 ) ;
1631 sp = tmp( start, end, inc ) ;
1632 fl = mFlagCol.getColumn() ;
1633 }
1634 else {
1635 for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1636 Bool crossOK = False ;
1637 Matrix<Complex> mSp = mDataCol( irow ) ;
1638 Matrix<Bool> mFl = mFlagCol( irow ) ;
1639 Matrix<Float> spxy = sp.xyPlane( irow ) ;
1640 Matrix<Bool> flxy = fl.xyPlane( irow ) ;
1641 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
1642 if ( corrtype[ipol] == Stokes::XY || corrtype[ipol] == Stokes::YX
1643 || corrtype[ipol] == Stokes::RL || corrtype[ipol] == Stokes::LR ) {
1644 if ( !crossOK ) {
1645 Vector<Float> tmp = ComplexToReal( mSp.row( ipol ) ) ;
1646 IPosition start( 1, 0 ) ;
1647 IPosition end( 1, 2*nchan-1 ) ;
1648 IPosition inc( 1, 2 ) ;
1649 spxy.row( ipol ) = tmp( start, end, inc ) ;
1650 flxy.row( ipol ) = mFl.row( ipol ) ;
1651 start = IPosition( 1, 1 ) ;
1652 spxy.row( ipol+1 ) = tmp( start, end, inc ) ;
1653 flxy.row( ipol+1 ) = mFl.row( ipol ) ;
1654 if ( corrtype[ipol] == Stokes::YX || corrtype[ipol] == Stokes::LR ) {
1655 spxy.row( ipol+1 ) = spxy.row( ipol+1 ) * (Float)-1.0 ;
1656 }
1657 crossOK = True ;
1658 }
1659 }
1660 else {
1661 Vector<Float> tmp = ComplexToReal( mSp.row( ipol ) ) ;
1662 IPosition start( 1, 0 ) ;
1663 IPosition end( 1, 2*nchan-1 ) ;
1664 IPosition inc( 1, 2 ) ;
1665 spxy.row( ipol ) = tmp( start, end, inc ) ;
1666 flxy.row( ipol ) = mFl.row( ipol ) ;
1667 }
1668 }
1669 }
1670 }
1671 }
1672 //double endSec = mathutil::gettimeofday_sec() ;
1673 //os_ << "end MSFiller::reshapeSpectraAndFlagtra() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1674}
1675
1676uInt MSFiller::getDirection( uInt idx,
1677 Vector<Double> &dir,
1678 Vector<Double> &azel,
1679 Vector<Double> &srate,
1680 Vector<Double> &ptcol,
1681 ROArrayColumn<Double> &pdcol,
1682 MEpoch &t,
1683 MPosition &antpos )
1684{
1685 //double startSec = mathutil::gettimeofday_sec() ;
1686 //os_ << "start MSFiller::getDirection2() startSec=" << startSec << LogIO::POST ;
1687 String refString ;
1688 MDirection::Types dirType ;
1689 uInt diridx = getDirection( idx, dir, srate, refString, ptcol, pdcol, t.get("s").getValue() ) ;
1690 MDirection::getType( dirType, refString ) ;
1691 MeasFrame mf( t, antpos ) ;
1692 if ( refString == "J2000" ) {
1693 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
1694 azel = toazel( dir ).getAngle("rad").getValue() ;
1695 }
1696 else if ( refString(0,4) == "AZEL" ) {
1697 azel = dir.copy() ;
1698 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
1699 dir = toj2000( dir ).getAngle("rad").getValue() ;
1700 }
1701 else {
1702 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
1703 azel = toazel( dir ).getAngle("rad").getValue() ;
1704 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
1705 dir = toj2000( dir ).getAngle("rad").getValue() ;
1706 }
1707 if ( srate.size() == 0 ) {
1708 srate.resize( 2 ) ;
1709 srate = 0.0 ;
1710 }
1711 //double endSec = mathutil::gettimeofday_sec() ;
1712 //os_ << "end MSFiller::getDirection2() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1713 return diridx ;
1714}
1715
1716void MSFiller::getSourceDirection( Vector<Double> &dir,
1717 Vector<Double> &azel,
1718 Vector<Double> &srate,
1719 MEpoch &t,
1720 MPosition &antpos,
1721 Vector<MDirection> &srcdir )
1722{
1723 //double startSec = mathutil::gettimeofday_sec() ;
1724 //os_ << "start MSFiller::getSourceDirection() startSec=" << startSec << LogIO::POST ;
1725 Vector<Double> defaultDir = srcdir[0].getAngle( "rad" ).getValue() ;
1726 if ( srcdir.nelements() > 1 )
1727 srate = srcdir[1].getAngle( "rad" ).getValue() ;
1728 String ref = srcdir[0].getRefString() ;
1729 MDirection::Types dirType ;
1730 MDirection::getType( dirType, ref ) ;
1731 MeasFrame mf( t, antpos ) ;
1732 if ( ref != "J2000" ) {
1733 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
1734 dir = toj2000( defaultDir ).getAngle("rad").getValue() ;
1735 }
1736 else
1737 dir = defaultDir ;
1738 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ;
1739 azel = toazel( defaultDir ).getAngle("rad").getValue() ;
1740 //double endSec = mathutil::gettimeofday_sec() ;
1741 //os_ << "end MSFiller::getSourceDirection() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1742}
1743
1744void MSFiller::initHeader( STHeader &header )
1745{
1746 header.nchan = 0 ;
1747 header.npol = 0 ;
1748 header.nif = 0 ;
1749 header.nbeam = 0 ;
1750 header.observer = "" ;
1751 header.project = "" ;
1752 header.obstype = "" ;
1753 header.antennaname = "" ;
1754 header.antennaposition.resize( 0 ) ;
1755 header.equinox = 0.0 ;
1756 header.freqref = "" ;
1757 header.reffreq = -1.0 ;
1758 header.bandwidth = 0.0 ;
1759 header.utc = 0.0 ;
1760 header.fluxunit = "" ;
1761 header.epoch = "" ;
1762 header.poltype = "" ;
1763}
1764
1765String MSFiller::asString( String name,
1766 uInt idx,
1767 Table tab,
1768 boost::object_pool<ROTableColumn> *pool )
1769{
1770 ROTableColumn *col = pool->construct( tab, name ) ;
1771 String v = col->asString( idx ) ;
1772 pool->destroy( col ) ;
1773 return v ;
1774}
1775
1776Bool MSFiller::asBool( String name,
1777 uInt idx,
1778 Table &tab,
1779 boost::object_pool<ROTableColumn> *pool )
1780{
1781 ROTableColumn *col = pool->construct( tab, name ) ;
1782 Bool v = col->asBool( idx ) ;
1783 pool->destroy( col ) ;
1784 return v ;
1785}
1786
1787uInt MSFiller::asuInt( String name,
1788 uInt idx,
1789 Table &tab,
1790 boost::object_pool<ROTableColumn> *pool )
1791{
1792 ROTableColumn *col = pool->construct( tab, name ) ;
1793 uInt v = col->asuInt( idx ) ;
1794 pool->destroy( col ) ;
1795 return v ;
1796}
1797
1798Int MSFiller::asInt( String name,
1799 uInt idx,
1800 Table &tab,
1801 boost::object_pool<ROTableColumn> *pool )
1802{
1803 ROTableColumn *col = pool->construct( tab, name ) ;
1804 Int v = col->asInt( idx ) ;
1805 pool->destroy( col ) ;
1806 return v ;
1807}
1808
1809Float MSFiller::asFloat( String name,
1810 uInt idx,
1811 Table &tab,
1812 boost::object_pool<ROTableColumn> *pool )
1813{
1814 ROTableColumn *col = pool->construct( tab, name ) ;
1815 Float v = col->asfloat( idx ) ;
1816 pool->destroy( col ) ;
1817 return v ;
1818}
1819
1820Double MSFiller::asDouble( String name,
1821 uInt idx,
1822 Table &tab,
1823 boost::object_pool<ROTableColumn> *pool )
1824{
1825 ROTableColumn *col = pool->construct( tab, name ) ;
1826 Double v = col->asdouble( idx ) ;
1827 pool->destroy( col ) ;
1828 return v ;
1829}
1830
1831void MSFiller::sourceInfo( Int sourceId,
1832 Int spwId,
1833 String &name,
1834 MDirection &direction,
1835 Vector<casa::Double> &properMotion,
1836 Vector<casa::Double> &restFreqs,
1837 Vector<casa::String> &transitions,
1838 Vector<casa::Double> &sysVels,
1839 boost::object_pool<ROTableColumn> *tpoolr )
1840{
1841 //double startSec = mathutil::gettimeofday_sec() ;
1842 //os_ << "start MSFiller::sourceInfo() startSec=" << startSec << LogIO::POST ;
1843
1844 MSSource srctab = mstable_.source() ;
1845 MSSource srctabSel = srctab( srctab.col("SOURCE_ID") == sourceId && srctab.col("SPECTRAL_WINDOW_ID") == spwId ) ;
1846 if ( srctabSel.nrow() == 0 ) {
1847 srctabSel = srctab( srctab.col("SOURCE_ID") == sourceId && srctab.col("SPECTRAL_WINDOW_ID") == -1 ) ;
1848 }
1849 Int numLines = 0 ;
1850 if ( srctabSel.nrow() > 0 ) {
1851 // source name
1852 name = asString( "NAME", 0, srctabSel, tpoolr ) ;
1853
1854 // source proper motion
1855 ROArrayColumn<Double> roArrDCol( srctabSel, "PROPER_MOTION" ) ;
1856 properMotion = roArrDCol( 0 ) ;
1857
1858 // source direction as MDirection object
1859 ROScalarMeasColumn<MDirection> tmpMeasCol( srctabSel, "DIRECTION" ) ;
1860 direction = tmpMeasCol( 0 ) ;
1861
1862 // number of lines
1863 numLines = asInt( "NUM_LINES", 0, srctabSel, tpoolr ) ;
1864 }
1865 else {
1866 name = "" ;
1867 properMotion = Vector<Double>( 2, 0.0 ) ;
1868 direction = MDirection( Quantum<Double>(0.0,Unit("rad")), Quantum<Double>(0.0,Unit("rad")) ) ;
1869 }
1870
1871 restFreqs.resize( numLines ) ;
1872 transitions.resize( numLines ) ;
1873 sysVels.resize( numLines ) ;
1874 if ( numLines > 0 ) {
1875 if ( srctabSel.tableDesc().isColumn( "REST_FREQUENCY" ) ) {
1876 ROArrayQuantColumn<Double> quantArrCol( srctabSel, "REST_FREQUENCY" ) ;
1877 Array< Quantum<Double> > qRestFreqs = quantArrCol( 0 ) ;
1878 for ( int i = 0 ; i < numLines ; i++ ) {
1879 restFreqs[i] = qRestFreqs( IPosition( 1, i ) ).getValue( "Hz" ) ;
1880 }
1881 }
1882 //os_ << "restFreqs = " << restFreqs << LogIO::POST ;
1883 if ( srctabSel.tableDesc().isColumn( "TRANSITION" ) ) {
1884 ROArrayColumn<String> transitionCol( srctabSel, "TRANSITION" ) ;
1885 if ( transitionCol.isDefined( 0 ) )
1886 transitions = transitionCol( 0 ) ;
1887 //os_ << "transitionNameCol.nrow() = " << transitionCol.nrow() << LogIO::POST ;
1888 }
1889 if ( srctabSel.tableDesc().isColumn( "SYSVEL" ) ) {
1890 ROArrayColumn<Double> roArrDCol( srctabSel, "SYSVEL" ) ;
1891 sysVels = roArrDCol( 0 ) ;
1892 }
1893 }
1894
1895 //double endSec = mathutil::gettimeofday_sec() ;
1896 //os_ << "end MSFiller::sourceInfo() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1897}
1898
1899void MSFiller::spectralSetup( Int spwId,
1900 MEpoch &me,
1901 MPosition &mp,
1902 MDirection &md,
1903 Double &refpix,
1904 Double &refval,
1905 Double &increment,
1906 Int &nchan,
1907 String &freqref,
1908 Double &reffreq,
1909 Double &bandwidth,
1910 boost::object_pool<ROTableColumn> *tpoolr )
1911{
1912 //double startSec = mathutil::gettimeofday_sec() ;
1913 //os_ << "start MSFiller::spectralSetup() startSec=" << startSec << LogIO::POST ;
1914
1915 MSSpectralWindow spwtab = mstable_.spectralWindow() ;
1916 MeasFrame mf( me, mp, md ) ;
1917 MFrequency::Types freqRef = MFrequency::castType( (uInt)asInt( "MEAS_FREQ_REF", spwId, spwtab, tpoolr ) ) ;
1918 Bool even = False ;
1919 if ( (nchan/2)*2 == nchan ) even = True ;
1920 ROScalarQuantColumn<Double> tmpQuantCol( spwtab, "TOTAL_BANDWIDTH" ) ;
1921 Double totbw = tmpQuantCol( spwId ).getValue( "Hz" ) ;
1922 if ( nchan != 4 )
1923 bandwidth = max( bandwidth, totbw ) ;
1924 if ( freqref == "" && nchan != 4)
1925 //sdh.freqref = MFrequency::showType( freqRef ) ;
1926 freqref = "LSRK" ;
1927 if ( reffreq == -1.0 && nchan != 4 ) {
1928 tmpQuantCol.attach( spwtab, "REF_FREQUENCY" ) ;
1929 Quantum<Double> qreffreq = tmpQuantCol( spwId ) ;
1930 if ( freqRef == MFrequency::LSRK ) {
1931 reffreq = qreffreq.getValue("Hz") ;
1932 }
1933 else {
1934 MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
1935 reffreq = tolsr( qreffreq ).get("Hz").getValue() ;
1936 }
1937 }
1938 Int refchan = nchan / 2 ;
1939 IPosition refip( 1, refchan ) ;
1940 refpix = 0.5*(nchan-1) ;
1941 refval = 0.0 ;
1942 ROArrayQuantColumn<Double> sharedQDArrCol( spwtab, "CHAN_WIDTH" ) ;
1943 ROTableColumn netSidebandCol( spwtab, "NET_SIDEBAND" ) ;
1944 Int netSideband = netSidebandCol.asInt( spwId ) ;
1945 increment = sharedQDArrCol( spwId )( refip ).getValue( "Hz" ) ;
1946 // os_ << "nchan = " << nchan << " refchan = " << refchan << "(even=" << even << ") refpix = " << refpix << LogIO::POST ;
1947 sharedQDArrCol.attach( spwtab, "CHAN_FREQ" ) ;
1948 Vector< Quantum<Double> > chanFreqs = sharedQDArrCol( spwId ) ;
1949 if ( ( nchan > 1 &&
1950 chanFreqs[0].getValue("Hz") > chanFreqs[1].getValue("Hz") )
1951 || ( nchan == 1 && netSideband == 1 ) )
1952 increment *= -1.0 ;
1953 if ( freqRef == MFrequency::LSRK ) {
1954 if ( even ) {
1955 IPosition refip0( 1, refchan-1 ) ;
1956 Double refval0 = chanFreqs(refip0).getValue("Hz") ;
1957 Double refval1 = chanFreqs(refip).getValue("Hz") ;
1958 refval = 0.5 * ( refval0 + refval1 ) ;
1959 }
1960 else {
1961 refval = chanFreqs(refip).getValue("Hz") ;
1962 }
1963 }
1964 else {
1965 MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
1966 if ( even ) {
1967 IPosition refip0( 1, refchan-1 ) ;
1968 Double refval0 = chanFreqs(refip0).getValue("Hz") ;
1969 Double refval1 = chanFreqs(refip).getValue("Hz") ;
1970 refval = 0.5 * ( refval0 + refval1 ) ;
1971 refval = tolsr( refval ).get("Hz").getValue() ;
1972 }
1973 else {
1974 refval = tolsr( chanFreqs(refip) ).get("Hz").getValue() ;
1975 }
1976 }
1977
1978 //double endSec = mathutil::gettimeofday_sec() ;
1979 //os_ << "end MSFiller::spectralSetup() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1980}
1981
1982} ;
1983
Note: See TracBrowser for help on using the repository browser.