source: trunk/src/MSFiller.cpp@ 2023

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

New Development: No

JIRA Issue: Yes CAS-2718

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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