source: tags/casa3.2.0asap/src/MSFiller.cpp@ 2747

Last change on this file since 2747 was 2191, checked in by Kana Sugimoto, 13 years ago

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