source: branches/parallelCasa3.3/src/MSWriter.cpp@ 2450

Last change on this file since 2450 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.1 KB
Line 
1//
2// C++ Interface: MSWriter
3//
4// Description:
5//
6// This class is specific writer 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 <casa/OS/File.h>
15#include <casa/OS/RegularFile.h>
16#include <casa/OS/Directory.h>
17#include <casa/OS/SymLink.h>
18#include <casa/BasicSL/String.h>
19#include <casa/Containers/RecordField.h>
20#include <casa/Arrays/Cube.h>
21
22#include <tables/Tables/ExprNode.h>
23#include <tables/Tables/TableDesc.h>
24#include <tables/Tables/SetupNewTab.h>
25#include <tables/Tables/TableIter.h>
26#include <tables/Tables/RefRows.h>
27#include <tables/Tables/TableRow.h>
28
29#include <ms/MeasurementSets/MeasurementSet.h>
30#include <ms/MeasurementSets/MSColumns.h>
31#include <ms/MeasurementSets/MSPolIndex.h>
32#include <ms/MeasurementSets/MSDataDescIndex.h>
33#include <ms/MeasurementSets/MSSourceIndex.h>
34
35#include "MSWriter.h"
36#include "STHeader.h"
37#include "STFrequencies.h"
38#include "STMolecules.h"
39#include "STTcal.h"
40
41#include <ctime>
42#include <sys/time.h>
43
44using namespace casa ;
45using namespace std ;
46
47namespace asap {
48double MSWriter::gettimeofday_sec()
49{
50 struct timeval tv ;
51 gettimeofday( &tv, NULL ) ;
52 return tv.tv_sec + (double)tv.tv_usec*1.0e-6 ;
53}
54
55MSWriter::MSWriter(CountedPtr<Scantable> stable)
56 : table_(stable),
57 isTcal_(False),
58 isWeather_(False),
59 tcalSpec_(False),
60 tsysSpec_(False),
61 ptTabName_("")
62{
63 os_ = LogIO() ;
64 os_.origin( LogOrigin( "MSWriter", "MSWriter()", WHERE ) ) ;
65// os_ << "MSWriter::MSWriter()" << LogIO::POST ;
66
67 // initialize writer
68 init() ;
69}
70
71MSWriter::~MSWriter()
72{
73 os_.origin( LogOrigin( "MSWriter", "~MSWriter()", WHERE ) ) ;
74// os_ << "MSWriter::~MSWriter()" << LogIO::POST ;
75
76 if ( mstable_ != 0 )
77 delete mstable_ ;
78}
79
80bool MSWriter::write(const string& filename, const Record& rec)
81{
82 os_.origin( LogOrigin( "MSWriter", "write()", WHERE ) ) ;
83// double startSec = gettimeofday_sec() ;
84// os_ << "start MSWriter::write() startSec=" << startSec << LogIO::POST ;
85
86 filename_ = filename ;
87
88 // parsing MS options
89 Bool overwrite = False ;
90 if ( rec.isDefined( "ms" ) ) {
91 Record msrec = rec.asRecord( "ms" ) ;
92 if ( msrec.isDefined( "overwrite" ) ) {
93 overwrite = msrec.asBool( "overwrite" ) ;
94 }
95 }
96
97 os_ << "Parsing MS options" << endl ;
98 os_ << " overwrite = " << overwrite << LogIO::POST ;
99
100 File file( filename_ ) ;
101 if ( file.exists() ) {
102 if ( overwrite ) {
103 os_ << filename_ << " exists. Overwrite existing data... " << LogIO::POST ;
104 if ( file.isRegular() ) RegularFile(file).remove() ;
105 else if ( file.isDirectory() ) Directory(file).removeRecursive() ;
106 else SymLink(file).remove() ;
107 }
108 else {
109 os_ << LogIO::SEVERE << "ERROR: " << filename_ << " exists..." << LogIO::POST ;
110 return False ;
111 }
112 }
113
114 // set up MS
115 setupMS() ;
116
117 // subtables
118 // OBSERVATION
119 fillObservation() ;
120
121 // ANTENNA
122 fillAntenna() ;
123
124 // PROCESSOR
125 fillProcessor() ;
126
127 // SOURCE
128 fillSource() ;
129
130 // WEATHER
131 if ( isWeather_ )
132 fillWeather() ;
133
134 // MAIN
135 // Iterate over several ids
136 Vector<uInt> processedFreqId( 0 ) ;
137 Int defaultFieldId = 0 ;
138
139 // row based
140 TableRow row( *mstable_ ) ;
141 TableRecord &trec = row.record() ;
142 NoticeTarget *dataRF = 0 ;
143 if ( useFloatData_ )
144 dataRF = new RecordFieldPtr< Array<Float> >( trec, "FLOAT_DATA" ) ;
145 else if ( useData_ )
146 dataRF = new RecordFieldPtr< Array<Complex> >( trec, "DATA" ) ;
147 RecordFieldPtr< Array<Bool> > flagRF( trec, "FLAG" ) ;
148 RecordFieldPtr<Bool> flagrowRF( trec, "FLAG_ROW" ) ;
149 RecordFieldPtr<Double> timeRF( trec, "TIME" ) ;
150 RecordFieldPtr<Double> timecRF( trec, "TIME_CENTROID" ) ;
151 RecordFieldPtr<Double> intervalRF( trec, "INTERVAL" ) ;
152 RecordFieldPtr<Double> exposureRF( trec, "EXPOSURE" ) ;
153 RecordFieldPtr< Array<Float> > weightRF( trec, "WEIGHT" ) ;
154 RecordFieldPtr< Array<Float> > sigmaRF( trec, "SIGMA" ) ;
155 RecordFieldPtr<Int> ddidRF( trec, "DATA_DESC_ID" ) ;
156 RecordFieldPtr<Int> stateidRF( trec, "STATE_ID" ) ;
157 RecordFieldPtr< Array<Bool> > flagcatRF( trec, "FLAG_CATEGORY" ) ;
158
159 // OBSERVATION_ID is always 0
160 RecordFieldPtr<Int> intRF( trec, "OBSERVATION_ID" ) ;
161 *intRF = 0 ;
162
163 // ANTENNA1 and ANTENNA2 are always 0
164 intRF.attachToRecord( trec, "ANTENNA1" ) ;
165 *intRF = 0 ;
166 intRF.attachToRecord( trec, "ANTENNA2" ) ;
167 *intRF = 0 ;
168
169 // ARRAY_ID is tentatively set to 0
170 intRF.attachToRecord( trec, "ARRAY_ID" ) ;
171 *intRF = 0 ;
172
173 // PROCESSOR_ID is tentatively set to 0
174 intRF.attachToRecord( trec, "PROCESSOR_ID" ) ;
175 *intRF = 0 ;
176
177 // UVW is always [0,0,0]
178 RecordFieldPtr< Array<Double> > uvwRF( trec, "UVW" ) ;
179 *uvwRF = Vector<Double>( 3, 0.0 ) ;
180
181 //
182 // ITERATION: FIELDNAME
183 //
184 TableIterator iter0( table_->table(), "FIELDNAME" ) ;
185 while( !iter0.pastEnd() ) {
186 //Table t0( iter0.table() ) ;
187 Table t0 = iter0.table() ;
188 ROTableColumn sharedCol( t0, "FIELDNAME" ) ;
189 String fieldName = sharedCol.asString( 0 ) ;
190 sharedCol.attach( t0, "SRCNAME" ) ;
191 String srcName = sharedCol.asString( 0 ) ;
192 sharedCol.attach( t0, "TIME" ) ;
193 Double minTime = (Double)sharedCol.asdouble( 0 ) * 86400.0 ; // day->sec
194 ROArrayColumn<Double> scanRateCol( t0, "SCANRATE" ) ;
195 Vector<Double> scanRate = scanRateCol( 0 ) ;
196 String::size_type pos = fieldName.find( "__" ) ;
197 Int fieldId = -1 ;
198 if ( pos != String::npos ) {
199// os_ << "fieldName.substr( pos+2 )=" << fieldName.substr( pos+2 ) << LogIO::POST ;
200 fieldId = String::toInt( fieldName.substr( pos+2 ) ) ;
201 fieldName = fieldName.substr( 0, pos ) ;
202 }
203 else {
204// os_ << "use default field id" << LogIO::POST ;
205 fieldId = defaultFieldId ;
206 defaultFieldId++ ;
207 }
208// os_ << "fieldId" << fieldId << ": " << fieldName << LogIO::POST ;
209
210 // FIELD_ID
211 intRF.attachToRecord( trec, "FIELD_ID" ) ;
212 *intRF = fieldId ;
213
214 //
215 // ITERATION: BEAMNO
216 //
217 TableIterator iter1( t0, "BEAMNO" ) ;
218 while( !iter1.pastEnd() ) {
219 Table t1 = iter1.table() ;
220 sharedCol.attach( t1, "BEAMNO" ) ;
221 uInt beamNo = sharedCol.asuInt( 0 ) ;
222// os_ << "beamNo = " << beamNo << LogIO::POST ;
223
224 // FEED1 and FEED2
225 intRF.attachToRecord( trec, "FEED1" ) ;
226 *intRF = beamNo ;
227 intRF.attachToRecord( trec, "FEED2" ) ;
228 *intRF = beamNo ;
229
230 //
231 // ITERATION: SCANNO
232 //
233 TableIterator iter2( t1, "SCANNO" ) ;
234 while( !iter2.pastEnd() ) {
235 Table t2 = iter2.table() ;
236 sharedCol.attach( t2, "SCANNO" ) ;
237 uInt scanNo = sharedCol.asuInt( 0 ) ;
238// os_ << "scanNo = " << scanNo << LogIO::POST ;
239
240 // SCAN_NUMBER
241 // MS: 1-based
242 // Scantable: 0-based
243 intRF.attachToRecord( trec, "SCAN_NUMBER" ) ;
244 *intRF = scanNo + 1 ;
245
246 //
247 // ITERATION: IFNO
248 //
249 TableIterator iter3( t2, "IFNO" ) ;
250 while( !iter3.pastEnd() ) {
251 Table t3 = iter3.table() ;
252 sharedCol.attach( t3, "IFNO" ) ;
253 uInt ifNo = sharedCol.asuInt( 0 ) ;
254// os_ << "ifNo = " << ifNo << LogIO::POST ;
255 sharedCol.attach( t3, "FREQ_ID" ) ;
256 uInt freqId = sharedCol.asuInt( 0 ) ;
257// os_ << "freqId = " << freqId << LogIO::POST ;
258 Int subscan = 1 ; // 1-base
259 //
260 // ITERATION: SRCTYPE
261 //
262 TableIterator iter4( t3, "SRCTYPE" ) ;
263 while( !iter4.pastEnd() ) {
264 Table t4 = iter4.table() ;
265 sharedCol.attach( t4, "SRCTYPE" ) ;
266 Int srcType = sharedCol.asInt( 0 ) ;
267 Int stateId = addState( srcType, subscan ) ;
268 *stateidRF = stateId ;
269 //
270 // ITERATION: CYCLENO and TIME
271 //
272 Block<String> cols( 2 ) ;
273 cols[0] = "CYCLENO" ;
274 cols[1] = "TIME" ;
275 TableIterator iter5( t4, cols ) ;
276 while( !iter5.pastEnd() ) {
277 Table t5 = iter5.table().sort("POLNO") ;
278 //sharedCol.attach( t5, "CYCLENO" ) ;
279 //uInt cycleNo = sharedCol.asuInt( 0 ) ;
280 Int nrow = t5.nrow() ;
281// os_ << "nrow = " << nrow << LogIO::POST ;
282
283 Vector<Int> polnos( nrow ) ;
284 indgen( polnos, 0 ) ;
285 Int polid = addPolarization( polnos ) ;
286// os_ << "polid = " << polid << LogIO::POST ;
287
288 // DATA/FLOAT_DATA
289 ROArrayColumn<Float> specCol( t5, "SPECTRA" ) ;
290 ROArrayColumn<uChar> flagCol( t5, "FLAGTRA" ) ;
291 uInt nchan = specCol( 0 ).size() ;
292 IPosition cellshape( 2, nrow, nchan ) ;
293 if ( useFloatData_ ) {
294 // FLOAT_DATA
295 Matrix<Float> dataArr( cellshape ) ;
296 Matrix<Bool> flagArr( cellshape ) ;
297 Vector<Bool> tmpB ;
298 for ( Int ipol = 0 ; ipol < nrow ; ipol++ ) {
299 dataArr.row( ipol ) = specCol( ipol ) ;
300 tmpB.reference( flagArr.row( ipol ) ) ;
301 convertArray( tmpB, flagCol( ipol ) ) ;
302 }
303 ((RecordFieldPtr< Array<Float> > *)dataRF)->define( dataArr ) ;
304
305 // FLAG
306 flagRF.define( flagArr ) ;
307 }
308 else if ( useData_ ) {
309 // DATA
310 // assume nrow = 4
311 Matrix<Complex> dataArr( cellshape ) ;
312 Vector<Float> zeroIm( nchan, 0 ) ;
313 Matrix<Float> dummy( IPosition( 2, 2, nchan ) ) ;
314 dummy.row( 0 ) = specCol( 0 ) ;
315 dummy.row( 1 ) = zeroIm ;
316 dataArr.row( 0 ) = RealToComplex( dummy ) ;
317 dummy.row( 0 ) = specCol( 1 ) ;
318 dataArr.row( 3 ) = RealToComplex( dummy ) ;
319 dummy.row( 0 ) = specCol( 2 ) ;
320 dummy.row( 1 ) = specCol( 3 ) ;
321 dataArr.row( 1 ) = RealToComplex( dummy ) ;
322 dataArr.row( 2 ) = conj( dataArr.row( 1 ) ) ;
323 ((RecordFieldPtr< Array<Complex> > *)dataRF)->define( dataArr ) ;
324
325
326 // FLAG
327 Matrix<Bool> flagArr( cellshape ) ;
328 Vector<Bool> tmpB ;
329 tmpB.reference( flagArr.row( 0 ) ) ;
330 convertArray( tmpB, flagCol( 0 ) ) ;
331 tmpB.reference( flagArr.row( 3 ) ) ;
332 convertArray( tmpB, flagCol( 3 ) ) ;
333 tmpB.reference( flagArr.row( 1 ) ) ;
334 convertArray( tmpB, ( flagCol( 2 ) | flagCol( 3 ) ) ) ;
335 flagArr.row( 2 ) = flagArr.row( 1 ) ;
336 flagRF.define( flagArr ) ;
337 }
338
339 // FLAG_ROW
340 sharedCol.attach( t5, "FLAGROW" ) ;
341 Vector<uInt> flagRowArr( nrow ) ;
342 for ( Int irow = 0 ; irow < nrow ; irow++ )
343 flagRowArr[irow] = sharedCol.asuInt( irow ) ;
344 *flagrowRF = anyNE( flagRowArr, (uInt)0 ) ;
345
346 // TIME and TIME_CENTROID
347 sharedCol.attach( t5, "TIME" ) ;
348 Double mTimeV = (Double)sharedCol.asdouble( 0 ) * 86400.0 ; // day -> sec
349 *timeRF = mTimeV ;
350 *timecRF = mTimeV ;
351
352 // INTERVAL and EXPOSURE
353 sharedCol.attach( t5, "INTERVAL" ) ;
354 Double interval = (Double)sharedCol.asdouble( 0 ) ;
355 *intervalRF = interval ;
356 *exposureRF = interval ;
357
358 // WEIGHT and SIGMA
359 // always 1 at the moment
360 Vector<Float> wArr( nrow, 1.0 ) ;
361 weightRF.define( wArr ) ;
362 sigmaRF.define( wArr ) ;
363
364 // add DATA_DESCRIPTION row
365 Int ddid = addDataDescription( polid, ifNo ) ;
366// os_ << "ddid = " << ddid << LogIO::POST ;
367 *ddidRF = ddid ;
368
369 // for SYSCAL table
370 sharedCol.attach( t5, "TCAL_ID" ) ;
371 Vector<uInt> tcalIdArr( nrow ) ;
372 for ( Int irow = 0 ; irow < nrow ; irow++ )
373 tcalIdArr[irow] = sharedCol.asuInt( irow ) ;
374// os_ << "tcalIdArr = " << tcalIdArr << LogIO::POST ;
375 String key = String::toString( tcalIdArr[0] ) ;
376 if ( !tcalIdRec_.isDefined( key ) ) {
377 tcalIdRec_.define( key, tcalIdArr ) ;
378 tcalRowRec_.define( key, t5.rowNumbers() ) ;
379 }
380 else {
381 Vector<uInt> pastrows = tcalRowRec_.asArrayuInt( key ) ;
382 tcalRowRec_.define( key, concatenateArray( pastrows, t5.rowNumbers() ) ) ;
383 }
384
385 // for POINTING table
386 if ( ptTabName_ == "" ) {
387 ROArrayColumn<Double> dirCol( t5, "DIRECTION" ) ;
388 Vector<Double> dir = dirCol( 0 ) ;
389 dirCol.attach( t5, "SCANRATE" ) ;
390 Vector<Double> rate = dirCol( 0 ) ;
391 Matrix<Double> msDir( 2, 1 ) ;
392 msDir.column( 0 ) = dir ;
393 if ( anyNE( rate, 0.0 ) ) {
394 msDir.resize( 2, 2 ) ;
395 msDir.column( 1 ) = rate ;
396 }
397 addPointing( fieldName, mTimeV, interval, msDir ) ;
398 }
399
400 // FLAG_CATEGORY is tentatively set
401 flagcatRF.define( Cube<Bool>( nrow, nchan, 1, False ) ) ;
402
403 // add row
404 mstable_->addRow( 1, True ) ;
405 row.put( mstable_->nrow()-1 ) ;
406
407 iter5.next() ;
408 }
409
410 iter4.next() ;
411 }
412
413 // add SPECTRAL_WINDOW row
414 if ( allNE( processedFreqId, freqId ) ) {
415 uInt vsize = processedFreqId.size() ;
416 processedFreqId.resize( vsize+1, True ) ;
417 processedFreqId[vsize] = freqId ;
418 addSpectralWindow( ifNo, freqId ) ;
419 }
420
421 iter3.next() ;
422 }
423
424 iter2.next() ;
425 }
426
427 // add FEED row
428 addFeed( beamNo ) ;
429
430 iter1.next() ;
431 }
432
433 // add FIELD row
434 addField( fieldId, fieldName, srcName, minTime, scanRate ) ;
435
436 iter0.next() ;
437 }
438
439// delete tpoolr ;
440 delete dataRF ;
441
442 // SYSCAL
443 if ( isTcal_ )
444 fillSysCal() ;
445
446 // fill empty SPECTRAL_WINDOW rows
447 infillSpectralWindow() ;
448
449 // ASDM tables
450 const TableRecord &stKeys = table_->table().keywordSet() ;
451 TableRecord &msKeys = mstable_->rwKeywordSet() ;
452 uInt nfields = stKeys.nfields() ;
453 for ( uInt ifield = 0 ; ifield < nfields ; ifield++ ) {
454 String kname = stKeys.name( ifield ) ;
455 if ( kname.find( "ASDM" ) != String::npos ) {
456 String asdmpath = stKeys.asString( ifield ) ;
457 os_ << "found ASDM table: " << asdmpath << LogIO::POST ;
458 if ( Table::isReadable( asdmpath ) ) {
459 Table newAsdmTab( asdmpath, Table::Old ) ;
460 newAsdmTab.copy( filename_+"/"+kname, Table::New ) ;
461 os_ << "add subtable: " << kname << LogIO::POST ;
462 msKeys.defineTable( kname, Table( filename_+"/"+kname, Table::Old ) ) ;
463 }
464 }
465 }
466
467 // replace POINTING table with original one if exists
468 if ( ptTabName_ != "" ) {
469 delete mstable_ ;
470 mstable_ = 0 ;
471 Table newPtTab( ptTabName_, Table::Old ) ;
472 newPtTab.copy( filename_+"/POINTING", Table::New ) ;
473 }
474
475// double endSec = gettimeofday_sec() ;
476// os_ << "end MSWriter::write() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
477
478 return True ;
479}
480
481void MSWriter::init()
482{
483// os_.origin( LogOrigin( "MSWriter", "init()", WHERE ) ) ;
484// double startSec = gettimeofday_sec() ;
485// os_ << "start MSWriter::init() startSec=" << startSec << LogIO::POST ;
486
487 // access to scantable
488 header_ = table_->getHeader() ;
489
490 // FLOAT_DATA? or DATA?
491 if ( header_.npol > 2 ) {
492 useFloatData_ = False ;
493 useData_ = True ;
494 }
495 else {
496 useFloatData_ = True ;
497 useData_ = False ;
498 }
499
500 // polarization type
501 polType_ = header_.poltype ;
502 if ( polType_ == "" )
503 polType_ = "stokes" ;
504 else if ( polType_.find( "linear" ) != String::npos )
505 polType_ = "linear" ;
506 else if ( polType_.find( "circular" ) != String::npos )
507 polType_ = "circular" ;
508 else if ( polType_.find( "stokes" ) != String::npos )
509 polType_ = "stokes" ;
510 else if ( polType_.find( "linpol" ) != String::npos )
511 polType_ = "linpol" ;
512 else
513 polType_ = "notype" ;
514
515 // Check if some subtables are exists
516 if ( table_->tcal().table().nrow() != 0 ) {
517 ROTableColumn col( table_->tcal().table(), "TCAL" ) ;
518 if ( col.isDefined( 0 ) ) {
519 os_ << "TCAL table exists: nrow=" << table_->tcal().table().nrow() << LogIO::POST ;
520 isTcal_ = True ;
521 }
522 else {
523 os_ << "No TCAL rows" << LogIO::POST ;
524 }
525 }
526 else {
527 os_ << "No TCAL rows" << LogIO::POST ;
528 }
529 if ( table_->weather().table().nrow() != 0 ) {
530 ROTableColumn col( table_->weather().table(), "TEMPERATURE" ) ;
531 if ( col.isDefined( 0 ) ) {
532 os_ << "WEATHER table exists: nrow=" << table_->weather().table().nrow() << LogIO::POST ;
533 isWeather_ =True ;
534 }
535 else {
536 os_ << "No WEATHER rows" << LogIO::POST ;
537 }
538 }
539 else {
540 os_ << "No WEATHER rows" << LogIO::POST ;
541 }
542
543 // Are TCAL_SPECTRUM and TSYS_SPECTRUM necessary?
544 if ( isTcal_ && header_.nchan != 1 ) {
545 // examine TCAL subtable
546 Table tcaltab = table_->tcal().table() ;
547 ROArrayColumn<Float> tcalCol( tcaltab, "TCAL" ) ;
548 for ( uInt irow = 0 ; irow < tcaltab.nrow() ; irow++ ) {
549 if ( tcalCol( irow ).size() != 1 )
550 tcalSpec_ = True ;
551 }
552 // examine spectral data
553 TableIterator iter0( table_->table(), "IFNO" ) ;
554 while( !iter0.pastEnd() ) {
555 Table t0( iter0.table() ) ;
556 ROArrayColumn<Float> sharedFloatArrCol( t0, "SPECTRA" ) ;
557 uInt len = sharedFloatArrCol( 0 ).size() ;
558 if ( len != 1 ) {
559 sharedFloatArrCol.attach( t0, "TSYS" ) ;
560 if ( sharedFloatArrCol( 0 ).size() != 1 )
561 tsysSpec_ = True ;
562 }
563 iter0.next() ;
564 }
565 }
566
567 // check if reference for POINTING table exists
568 const TableRecord &rec = table_->table().keywordSet() ;
569 if ( rec.isDefined( "POINTING" ) ) {
570 ptTabName_ = rec.asString( "POINTING" ) ;
571 if ( !Table::isReadable( ptTabName_ ) ) {
572 ptTabName_ = "" ;
573 }
574 }
575
576// double endSec = gettimeofday_sec() ;
577// os_ << "end MSWriter::init() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
578}
579
580void MSWriter::setupMS()
581{
582// os_.origin( LogOrigin( "MSWriter", "setupMS()", WHERE ) ) ;
583// double startSec = gettimeofday_sec() ;
584// os_ << "start MSWriter::setupMS() startSec=" << startSec << LogIO::POST ;
585
586 TableDesc msDesc = MeasurementSet::requiredTableDesc() ;
587 if ( useFloatData_ )
588 MeasurementSet::addColumnToDesc( msDesc, MSMainEnums::FLOAT_DATA, 2 ) ;
589 else if ( useData_ )
590 MeasurementSet::addColumnToDesc( msDesc, MSMainEnums::DATA, 2 ) ;
591
592 SetupNewTable newtab( filename_, msDesc, Table::New ) ;
593
594 mstable_ = new MeasurementSet( newtab ) ;
595
596 // create subtables
597 TableDesc antennaDesc = MSAntenna::requiredTableDesc() ;
598 SetupNewTable antennaTab( mstable_->antennaTableName(), antennaDesc, Table::New ) ;
599 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::ANTENNA ), Table( antennaTab ) ) ;
600
601 TableDesc dataDescDesc = MSDataDescription::requiredTableDesc() ;
602 SetupNewTable dataDescTab( mstable_->dataDescriptionTableName(), dataDescDesc, Table::New ) ;
603 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::DATA_DESCRIPTION ), Table( dataDescTab ) ) ;
604
605 TableDesc dopplerDesc = MSDoppler::requiredTableDesc() ;
606 SetupNewTable dopplerTab( mstable_->dopplerTableName(), dopplerDesc, Table::New ) ;
607 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::DOPPLER ), Table( dopplerTab ) ) ;
608
609 TableDesc feedDesc = MSFeed::requiredTableDesc() ;
610 SetupNewTable feedTab( mstable_->feedTableName(), feedDesc, Table::New ) ;
611 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FEED ), Table( feedTab ) ) ;
612
613 TableDesc fieldDesc = MSField::requiredTableDesc() ;
614 SetupNewTable fieldTab( mstable_->fieldTableName(), fieldDesc, Table::New ) ;
615 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FIELD ), Table( fieldTab ) ) ;
616
617 TableDesc flagCmdDesc = MSFlagCmd::requiredTableDesc() ;
618 SetupNewTable flagCmdTab( mstable_->flagCmdTableName(), flagCmdDesc, Table::New ) ;
619 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FLAG_CMD ), Table( flagCmdTab ) ) ;
620
621 TableDesc freqOffsetDesc = MSFreqOffset::requiredTableDesc() ;
622 SetupNewTable freqOffsetTab( mstable_->freqOffsetTableName(), freqOffsetDesc, Table::New ) ;
623 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FREQ_OFFSET ), Table( freqOffsetTab ) ) ;
624
625 TableDesc historyDesc = MSHistory::requiredTableDesc() ;
626 SetupNewTable historyTab( mstable_->historyTableName(), historyDesc, Table::New ) ;
627 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::HISTORY ), Table( historyTab ) ) ;
628
629 TableDesc observationDesc = MSObservation::requiredTableDesc() ;
630 SetupNewTable observationTab( mstable_->observationTableName(), observationDesc, Table::New ) ;
631 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::OBSERVATION ), Table( observationTab ) ) ;
632
633 TableDesc pointingDesc = MSPointing::requiredTableDesc() ;
634 SetupNewTable pointingTab( mstable_->pointingTableName(), pointingDesc, Table::New ) ;
635 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::POINTING ), Table( pointingTab ) ) ;
636
637 TableDesc polarizationDesc = MSPolarization::requiredTableDesc() ;
638 SetupNewTable polarizationTab( mstable_->polarizationTableName(), polarizationDesc, Table::New ) ;
639 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::POLARIZATION ), Table( polarizationTab ) ) ;
640
641 TableDesc processorDesc = MSProcessor::requiredTableDesc() ;
642 SetupNewTable processorTab( mstable_->processorTableName(), processorDesc, Table::New ) ;
643 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::PROCESSOR ), Table( processorTab ) ) ;
644
645 TableDesc sourceDesc = MSSource::requiredTableDesc() ;
646 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::TRANSITION, 1 ) ;
647 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::REST_FREQUENCY, 1 ) ;
648 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::SYSVEL, 1 ) ;
649 SetupNewTable sourceTab( mstable_->sourceTableName(), sourceDesc, Table::New ) ;
650 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SOURCE ), Table( sourceTab ) ) ;
651
652 TableDesc spwDesc = MSSpectralWindow::requiredTableDesc() ;
653 SetupNewTable spwTab( mstable_->spectralWindowTableName(), spwDesc, Table::New ) ;
654 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SPECTRAL_WINDOW ), Table( spwTab ) ) ;
655
656 TableDesc stateDesc = MSState::requiredTableDesc() ;
657 SetupNewTable stateTab( mstable_->stateTableName(), stateDesc, Table::New ) ;
658 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::STATE ), Table( stateTab ) ) ;
659
660 // TODO: add TCAL_SPECTRUM and TSYS_SPECTRUM if necessary
661 TableDesc sysCalDesc = MSSysCal::requiredTableDesc() ;
662 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TCAL, 2 ) ;
663 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TSYS, 2 ) ;
664 if ( tcalSpec_ )
665 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TCAL_SPECTRUM, 2 ) ;
666 if ( tsysSpec_ )
667 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TSYS_SPECTRUM, 2 ) ;
668 SetupNewTable sysCalTab( mstable_->sysCalTableName(), sysCalDesc, Table::New ) ;
669 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SYSCAL ), Table( sysCalTab ) ) ;
670
671 TableDesc weatherDesc = MSWeather::requiredTableDesc() ;
672 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::TEMPERATURE ) ;
673 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::PRESSURE ) ;
674 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::REL_HUMIDITY ) ;
675 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::WIND_SPEED ) ;
676 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::WIND_DIRECTION ) ;
677 SetupNewTable weatherTab( mstable_->weatherTableName(), weatherDesc, Table::New ) ;
678 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::WEATHER ), Table( weatherTab ) ) ;
679
680 mstable_->initRefs() ;
681
682// double endSec = gettimeofday_sec() ;
683// os_ << "end MSWriter::setupMS() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
684}
685
686void MSWriter::fillObservation()
687{
688// double startSec = gettimeofday_sec() ;
689// os_ << "start MSWriter::fillObservation() startSec=" << startSec << LogIO::POST ;
690
691 // only 1 row
692 mstable_->observation().addRow( 1, True ) ;
693 MSObservationColumns msObsCols( mstable_->observation() ) ;
694 msObsCols.observer().put( 0, header_.observer ) ;
695 // tentatively put antennaname (from ANTENNA subtable)
696 String hAntennaName = header_.antennaname ;
697 String::size_type pos = hAntennaName.find( "//" ) ;
698 String telescopeName ;
699 if ( pos != String::npos ) {
700 telescopeName = hAntennaName.substr( 0, pos ) ;
701 }
702 else {
703 pos = hAntennaName.find( "@" ) ;
704 telescopeName = hAntennaName.substr( 0, pos ) ;
705 }
706// os_ << "telescopeName = " << telescopeName << LogIO::POST ;
707 msObsCols.telescopeName().put( 0, telescopeName ) ;
708 msObsCols.project().put( 0, header_.project ) ;
709 //ScalarMeasColumn<MEpoch> timeCol( table_->table().sort("TIME"), "TIME" ) ;
710 Table sortedtable = table_->table().sort("TIME") ;
711 ScalarMeasColumn<MEpoch> timeCol( sortedtable, "TIME" ) ;
712 Vector<MEpoch> trange( 2 ) ;
713 trange[0] = timeCol( 0 ) ;
714 trange[1] = timeCol( table_->nrow()-1 ) ;
715 msObsCols.timeRangeMeas().put( 0, trange ) ;
716
717// double endSec = gettimeofday_sec() ;
718// os_ << "end MSWriter::fillObservation() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
719}
720
721void MSWriter::fillAntenna()
722{
723// double startSec = gettimeofday_sec() ;
724// os_ << "start MSWriter::fillAntenna() startSec=" << startSec << LogIO::POST ;
725
726 // only 1 row
727 mstable_->antenna().addRow( 1, True ) ;
728 MSAntennaColumns msAntCols( mstable_->antenna() ) ;
729
730 String hAntennaName = header_.antennaname ;
731 String::size_type pos = hAntennaName.find( "//" ) ;
732 String antennaName ;
733 String stationName ;
734 if ( pos != String::npos ) {
735 hAntennaName = hAntennaName.substr( pos+2 ) ;
736 }
737 pos = hAntennaName.find( "@" ) ;
738 if ( pos != String::npos ) {
739 antennaName = hAntennaName.substr( 0, pos ) ;
740 stationName = hAntennaName.substr( pos+1 ) ;
741 }
742 else {
743 antennaName = hAntennaName ;
744 stationName = hAntennaName ;
745 }
746// os_ << "antennaName = " << antennaName << LogIO::POST ;
747// os_ << "stationName = " << stationName << LogIO::POST ;
748
749 msAntCols.name().put( 0, antennaName ) ;
750 msAntCols.station().put( 0, stationName ) ;
751
752// os_ << "antennaPosition = " << header_.antennaposition << LogIO::POST ;
753
754 msAntCols.position().put( 0, header_.antennaposition ) ;
755
756 // MOUNT is set to "ALT-AZ"
757 msAntCols.mount().put( 0, "ALT-AZ" ) ;
758
759 Double diameter = getDishDiameter( antennaName ) ;
760 msAntCols.dishDiameterQuant().put( 0, Quantity( diameter, "m" ) ) ;
761
762// double endSec = gettimeofday_sec() ;
763// os_ << "end MSWriter::fillAntenna() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
764}
765
766void MSWriter::fillProcessor()
767{
768// double startSec = gettimeofday_sec() ;
769// os_ << "start MSWriter::fillProcessor() startSec=" << startSec << LogIO::POST ;
770
771 // only add empty 1 row
772 MSProcessor msProc = mstable_->processor() ;
773 msProc.addRow( 1, True ) ;
774
775// double endSec = gettimeofday_sec() ;
776// os_ << "end MSWriter::fillProcessor() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
777}
778
779void MSWriter::fillSource()
780{
781// double startSec = gettimeofday_sec() ;
782// os_ << "start MSWriter::fillSource() startSec=" << startSec << LogIO::POST ;
783
784 // access to MS SOURCE subtable
785 MSSource msSrc = mstable_->source() ;
786
787 // access to MOLECULE subtable
788 STMolecules stm = table_->molecules() ;
789
790 Int srcId = 0 ;
791 Vector<Double> restFreq ;
792 Vector<String> molName ;
793 Vector<String> fMolName ;
794
795 // row based
796 TableRow row( msSrc ) ;
797 TableRecord &rec = row.record() ;
798 RecordFieldPtr<Int> srcidRF( rec, "SOURCE_ID" ) ;
799 RecordFieldPtr<String> nameRF( rec, "NAME" ) ;
800 RecordFieldPtr< Array<Double> > srcpmRF( rec, "PROPER_MOTION" ) ;
801 RecordFieldPtr< Array<Double> > srcdirRF( rec, "DIRECTION" ) ;
802 RecordFieldPtr<Int> numlineRF( rec, "NUM_LINES" ) ;
803 RecordFieldPtr< Array<Double> > restfreqRF( rec, "REST_FREQUENCY" ) ;
804 RecordFieldPtr< Array<Double> > sysvelRF( rec, "SYSVEL" ) ;
805 RecordFieldPtr< Array<String> > transitionRF( rec, "TRANSITION" ) ;
806 RecordFieldPtr<Double> timeRF( rec, "TIME" ) ;
807 RecordFieldPtr<Double> intervalRF( rec, "INTERVAL" ) ;
808 RecordFieldPtr<Int> spwidRF( rec, "SPECTRAL_WINDOW_ID" ) ;
809
810 //
811 // ITERATION: SRCNAME
812 //
813 TableIterator iter0( table_->table(), "SRCNAME" ) ;
814 while( !iter0.pastEnd() ) {
815 //Table t0( iter0.table() ) ;
816 Table t0 = iter0.table() ;
817
818 // get necessary information
819 ROScalarColumn<String> srcNameCol( t0, "SRCNAME" ) ;
820 String srcName = srcNameCol( 0 ) ;
821 ROArrayColumn<Double> sharedDArrRCol( t0, "SRCPROPERMOTION" ) ;
822 Vector<Double> srcPM = sharedDArrRCol( 0 ) ;
823 sharedDArrRCol.attach( t0, "SRCDIRECTION" ) ;
824 Vector<Double> srcDir = sharedDArrRCol( 0 ) ;
825 ROScalarColumn<Double> srcVelCol( t0, "SRCVELOCITY" ) ;
826 Double srcVel = srcVelCol( 0 ) ;
827
828 // NAME
829 *nameRF = srcName ;
830
831 // SOURCE_ID
832 *srcidRF = srcId ;
833
834 // PROPER_MOTION
835 *srcpmRF = srcPM ;
836
837 // DIRECTION
838 *srcdirRF = srcDir ;
839
840 //
841 // ITERATION: MOLECULE_ID
842 //
843 TableIterator iter1( t0, "MOLECULE_ID" ) ;
844 while( !iter1.pastEnd() ) {
845 //Table t1( iter1.table() ) ;
846 Table t1 = iter1.table() ;
847
848 // get necessary information
849 ROScalarColumn<uInt> molIdCol( t1, "MOLECULE_ID" ) ;
850 uInt molId = molIdCol( 0 ) ;
851 stm.getEntry( restFreq, molName, fMolName, molId ) ;
852
853 uInt numFreq = restFreq.size() ;
854
855 // NUM_LINES
856 *numlineRF = numFreq ;
857
858 // REST_FREQUENCY
859 *restfreqRF = restFreq ;
860
861 // TRANSITION
862 //*transitionRF = fMolName ;
863 Vector<String> transition ;
864 if ( fMolName.size() != 0 ) {
865 transition = fMolName ;
866 }
867 else if ( molName.size() != 0 ) {
868 transition = molName ;
869 }
870 else {
871 transition.resize( numFreq ) ;
872 transition = "" ;
873 }
874 *transitionRF = transition ;
875
876 // SYSVEL
877 Vector<Double> sysvelArr( numFreq, srcVel ) ;
878 *sysvelRF = sysvelArr ;
879
880 //
881 // ITERATION: IFNO
882 //
883 TableIterator iter2( t1, "IFNO" ) ;
884 while( !iter2.pastEnd() ) {
885 //Table t2( iter2.table() ) ;
886 Table t2 = iter2.table() ;
887 uInt nrow = msSrc.nrow() ;
888
889 // get necessary information
890 ROScalarColumn<uInt> ifNoCol( t2, "IFNO" ) ;
891 uInt ifno = ifNoCol( 0 ) ; // IFNO = SPECTRAL_WINDOW_ID
892 Double midTime ;
893 Double interval ;
894 getValidTimeRange( midTime, interval, t2 ) ;
895
896 // fill SPECTRAL_WINDOW_ID
897 *spwidRF = ifno ;
898
899 // fill TIME, INTERVAL
900 *timeRF = midTime ;
901 *intervalRF = interval ;
902
903 // add row
904 msSrc.addRow( 1, True ) ;
905 row.put( nrow ) ;
906
907 iter2.next() ;
908 }
909
910 iter1.next() ;
911 }
912
913 // increment srcId if SRCNAME changed
914 srcId++ ;
915
916 iter0.next() ;
917 }
918
919// double endSec = gettimeofday_sec() ;
920// os_ << "end MSWriter::fillSource() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
921}
922
923void MSWriter::fillWeather()
924{
925// double startSec = gettimeofday_sec() ;
926// os_ << "start MSWriter::fillWeather() startSec=" << startSec << LogIO::POST ;
927
928 // access to MS WEATHER subtable
929 MSWeather msw = mstable_->weather() ;
930
931 // access to WEATHER subtable
932 Table stw = table_->weather().table() ;
933 uInt nrow = stw.nrow() ;
934
935 if ( nrow == 0 )
936 return ;
937
938 msw.addRow( nrow, True ) ;
939 MSWeatherColumns mswCols( msw ) ;
940
941 // ANTENNA_ID is always 0
942 Vector<Int> antIdArr( nrow, 0 ) ;
943 mswCols.antennaId().putColumn( antIdArr ) ;
944
945 // fill weather status
946 ROScalarColumn<Float> sharedFloatCol( stw, "TEMPERATURE" ) ;
947 mswCols.temperature().putColumn( sharedFloatCol ) ;
948 sharedFloatCol.attach( stw, "PRESSURE" ) ;
949 mswCols.pressure().putColumn( sharedFloatCol ) ;
950 sharedFloatCol.attach( stw, "HUMIDITY" ) ;
951 mswCols.relHumidity().putColumn( sharedFloatCol ) ;
952 sharedFloatCol.attach( stw, "WINDSPEED" ) ;
953 mswCols.windSpeed().putColumn( sharedFloatCol ) ;
954 sharedFloatCol.attach( stw, "WINDAZ" ) ;
955 mswCols.windDirection().putColumn( sharedFloatCol ) ;
956
957 // fill TIME and INTERVAL
958 Double midTime ;
959 Double interval ;
960 Vector<Double> intervalArr( nrow, 0.0 ) ;
961 TableIterator iter( table_->table(), "WEATHER_ID" ) ;
962 while( !iter.pastEnd() ) {
963 //Table tab( iter.table() ) ;
964 Table tab = iter.table() ;
965
966 ROScalarColumn<uInt> widCol( tab, "WEATHER_ID" ) ;
967 uInt wid = widCol( 0 ) ;
968
969 getValidTimeRange( midTime, interval, tab ) ;
970 mswCols.time().put( wid, midTime ) ;
971 intervalArr[wid] = interval ;
972
973 iter.next() ;
974 }
975 mswCols.interval().putColumn( intervalArr ) ;
976
977// double endSec = gettimeofday_sec() ;
978// os_ << "end MSWriter::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
979}
980
981void MSWriter::fillSysCal()
982{
983// double startSec = gettimeofday_sec() ;
984// os_ << "start MSWriter::fillSysCal() startSec=" << startSec << LogIO::POST ;
985
986 //tcalIdRec_.print( cout ) ;
987
988 // access to MS SYSCAL subtable
989 MSSysCal mssc = mstable_->sysCal() ;
990
991 // access to TCAL subtable
992 Table stt = table_->tcal().table() ;
993 uInt nrow = stt.nrow() ;
994
995 // access to MAIN table
996 Block<String> cols( 6 ) ;
997 cols[0] = "TIME" ;
998 cols[1] = "TCAL_ID" ;
999 cols[2] = "TSYS" ;
1000 cols[3] = "BEAMNO" ;
1001 cols[4] = "IFNO" ;
1002 cols[5] = "INTERVAL" ;
1003 Table tab = table_->table().project( cols ) ;
1004
1005 if ( nrow == 0 )
1006 return ;
1007
1008 nrow = tcalIdRec_.nfields() ;
1009
1010 Double midTime ;
1011 Double interval ;
1012 String timeStr ;
1013
1014 // row base
1015 TableRow row( mssc ) ;
1016 TableRecord &trec = row.record() ;
1017 RecordFieldPtr<Int> antennaRF( trec, "ANTENNA_ID" ) ;
1018 RecordFieldPtr<Int> feedRF( trec, "FEED_ID" ) ;
1019 RecordFieldPtr<Int> spwRF( trec, "SPECTRAL_WINDOW_ID" ) ;
1020 RecordFieldPtr<Double> timeRF( trec, "TIME" ) ;
1021 RecordFieldPtr<Double> intervalRF( trec, "INTERVAL" ) ;
1022 RecordFieldPtr< Array<Float> > tsysRF( trec, "TSYS" ) ;
1023 RecordFieldPtr< Array<Float> > tcalRF( trec, "TCAL" ) ;
1024 RecordFieldPtr< Array<Float> > tsysspRF ;
1025 RecordFieldPtr< Array<Float> > tcalspRF ;
1026 if ( tsysSpec_ )
1027 tsysspRF.attachToRecord( trec, "TSYS_SPECTRUM" ) ;
1028 if ( tcalSpec_ )
1029 tcalspRF.attachToRecord( trec, "TCAL_SPECTRUM" ) ;
1030
1031 // ANTENNA_ID is always 0
1032 *antennaRF = 0 ;
1033
1034 Table sortedstt = stt.sort( "ID" ) ;
1035 ROArrayColumn<Float> tcalCol( sortedstt, "TCAL" ) ;
1036 ROTableColumn idCol( sortedstt, "ID" ) ;
1037 ROArrayColumn<Float> tsysCol( tab, "TSYS" ) ;
1038 ROTableColumn tcalidCol( tab, "TCAL_ID" ) ;
1039 ROTableColumn timeCol( tab, "TIME" ) ;
1040 ROTableColumn intervalCol( tab, "INTERVAL" ) ;
1041 ROTableColumn beamnoCol( tab, "BEAMNO" ) ;
1042 ROTableColumn ifnoCol( tab, "IFNO" ) ;
1043 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
1044// double t1 = gettimeofday_sec() ;
1045 Vector<uInt> ids = tcalIdRec_.asArrayuInt( irow ) ;
1046// os_ << "ids = " << ids << LogIO::POST ;
1047 uInt npol = ids.size() ;
1048 Vector<uInt> rows = tcalRowRec_.asArrayuInt( irow ) ;
1049// os_ << "rows = " << rows << LogIO::POST ;
1050 Vector<Double> atime( rows.nelements() ) ;
1051 Vector<Double> ainterval( rows.nelements() ) ;
1052 Vector<uInt> atcalid( rows.nelements() ) ;
1053 for( uInt jrow = 0 ; jrow < rows.nelements() ; jrow++ ) {
1054 atime[jrow] = (Double)timeCol.asdouble( rows[jrow] ) ;
1055 ainterval[jrow] = (Double)intervalCol.asdouble( rows[jrow] ) ;
1056 atcalid[jrow] = tcalidCol.asuInt( rows[jrow] ) ;
1057 }
1058 Vector<Float> dummy = tsysCol( rows[0] ) ;
1059 Matrix<Float> tsys( npol,dummy.nelements() ) ;
1060 tsys.row( 0 ) = dummy ;
1061 for ( uInt jrow = 1 ; jrow < npol ; jrow++ )
1062 tsys.row( jrow ) = tsysCol( rows[jrow] ) ;
1063
1064 // FEED_ID
1065 *feedRF = beamnoCol.asuInt( rows[0] ) ;
1066
1067 // SPECTRAL_WINDOW_ID
1068 *spwRF = ifnoCol.asuInt( rows[0] ) ;
1069
1070 // TIME and INTERVAL
1071 getValidTimeRange( midTime, interval, atime, ainterval ) ;
1072 *timeRF = midTime ;
1073 *intervalRF = interval ;
1074
1075 // TCAL and TSYS
1076 Matrix<Float> tcal ;
1077 Table t ;
1078 if ( idCol.asuInt( ids[0] ) == ids[0] ) {
1079// os_ << "sorted at irow=" << irow << " ids[0]=" << ids[0] << LogIO::POST ;
1080 Vector<Float> dummyC = tcalCol( ids[0] ) ;
1081 tcal.resize( npol, dummyC.size() ) ;
1082 tcal.row( 0 ) = dummyC ;
1083 }
1084 else {
1085// os_ << "NOT sorted at irow=" << irow << " ids[0]=" << ids[0] << LogIO::POST ;
1086 t = stt( stt.col("ID") == ids[0] ) ;
1087 Vector<Float> dummyC = tcalCol( 0 ) ;
1088 tcal.resize( npol, dummyC.size() ) ;
1089 tcal.row( 0 ) = dummyC ;
1090 }
1091 if ( npol == 2 ) {
1092 if ( idCol.asuInt( ids[1] ) == ids[1] ) {
1093// os_ << "sorted at irow=" << irow << " ids[1]=" << ids[1] << LogIO::POST ;
1094 tcal.row( 1 ) = tcalCol( ids[1] ) ;
1095 }
1096 else {
1097// os_ << "NOT sorted at irow=" << irow << " ids[1]=" << ids[1] << LogIO::POST ;
1098 t = stt( stt.col("ID") == ids[1] ) ;
1099 tcalCol.attach( t, "TCAL" ) ;
1100 tcal.row( 1 ) = tcalCol( 1 ) ;
1101 }
1102 }
1103 else if ( npol == 3 ) {
1104 if ( idCol.asuInt( ids[2] ) == ids[2] )
1105 tcal.row( 1 ) = tcalCol( ids[2] ) ;
1106 else {
1107 t = stt( stt.col("ID") == ids[2] ) ;
1108 tcalCol.attach( t, "TCAL" ) ;
1109 tcal.row( 1 ) = tcalCol( 0 ) ;
1110 }
1111 if ( idCol.asuInt( ids[1] ) == ids[1] )
1112 tcal.row( 2 ) = tcalCol( ids[1] ) ;
1113 else {
1114 t = stt( stt.col("ID") == ids[1] ) ;
1115 tcalCol.attach( t, "TCAL" ) ;
1116 tcal.row( 2 ) = tcalCol( 0 ) ;
1117 }
1118 }
1119 else if ( npol == 4 ) {
1120 if ( idCol.asuInt( ids[2] ) == ids[2] )
1121 tcal.row( 1 ) = tcalCol( ids[2] ) ;
1122 else {
1123 t = stt( stt.col("ID") == ids[2] ) ;
1124 tcalCol.attach( t, "TCAL" ) ;
1125 tcal.row( 1 ) = tcalCol( 0 ) ;
1126 }
1127 if ( idCol.asuInt( ids[3] ) == ids[3] )
1128 tcal.row( 2 ) = tcalCol( ids[3] ) ;
1129 else {
1130 t = stt( stt.col("ID") == ids[3] ) ;
1131 tcalCol.attach( t, "TCAL" ) ;
1132 tcal.row( 2 ) = tcalCol( 0 ) ;
1133 }
1134 if ( idCol.asuInt( ids[1] ) == ids[1] )
1135 tcal.row( 2 ) = tcalCol( ids[1] ) ;
1136 else {
1137 t = stt( stt.col("ID") == ids[1] ) ;
1138 tcalCol.attach( t, "TCAL" ) ;
1139 tcal.row( 3 ) = tcalCol( 0 ) ;
1140 }
1141 }
1142 if ( tcalSpec_ ) {
1143 // put TCAL_SPECTRUM
1144 //*tcalspRF = tcal ;
1145 tcalspRF.define( tcal ) ;
1146 // set TCAL (mean of TCAL_SPECTRUM)
1147 Matrix<Float> tcalMean( npol, 1 ) ;
1148 for ( uInt iid = 0 ; iid < npol ; iid++ ) {
1149 tcalMean( iid, 0 ) = mean( tcal.row(iid) ) ;
1150 }
1151 // put TCAL
1152 *tcalRF = tcalMean ;
1153 }
1154 else {
1155 // put TCAL
1156 *tcalRF = tcal ;
1157 }
1158
1159 if ( tsysSpec_ ) {
1160 // put TSYS_SPECTRUM
1161 //*tsysspRF = tsys ;
1162 tsysspRF.define( tsys ) ;
1163 // set TSYS (mean of TSYS_SPECTRUM)
1164 Matrix<Float> tsysMean( npol, 1 ) ;
1165 for ( uInt iid = 0 ; iid < npol ; iid++ ) {
1166 tsysMean( iid, 0 ) = mean( tsys.row(iid) ) ;
1167 }
1168 // put TSYS
1169 *tsysRF = tsysMean ;
1170 }
1171 else {
1172 // put TSYS
1173 *tsysRF = tsys ;
1174 }
1175
1176 // add row
1177 mssc.addRow( 1, True ) ;
1178 row.put( mssc.nrow()-1 ) ;
1179
1180// double t2 = gettimeofday_sec() ;
1181// os_ << irow << "th loop elapsed time = " << t2-t1 << "sec" << LogIO::POST ;
1182 }
1183
1184// double endSec = gettimeofday_sec() ;
1185// os_ << "end MSWriter::fillSysCal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1186}
1187
1188void MSWriter::addFeed( Int id )
1189{
1190// double startSec = gettimeofday_sec() ;
1191// os_ << "start MSWriter::addFeed() startSec=" << startSec << LogIO::POST ;
1192
1193 // add row
1194 MSFeed msFeed = mstable_->feed() ;
1195 msFeed.addRow( 1, True ) ;
1196 Int nrow = msFeed.nrow() ;
1197 Int numReceptors = 2 ;
1198 Vector<String> polType( numReceptors ) ;
1199 Matrix<Double> beamOffset( 2, numReceptors ) ;
1200 beamOffset = 0.0 ;
1201 Vector<Double> receptorAngle( numReceptors, 0.0 ) ;
1202 if ( polType_ == "linear" ) {
1203 polType[0] = "X" ;
1204 polType[1] = "Y" ;
1205 }
1206 else if ( polType_ == "circular" ) {
1207 polType[0] = "R" ;
1208 polType[1] = "L" ;
1209 }
1210 else {
1211 polType[0] = "X" ;
1212 polType[1] = "Y" ;
1213 }
1214 Matrix<Complex> polResponse( numReceptors, numReceptors, 0.0 ) ;
1215 for ( Int i = 0 ; i < numReceptors ; i++ )
1216 polResponse( i, i ) = 0.0 ;
1217
1218 MSFeedColumns msFeedCols( mstable_->feed() ) ;
1219
1220 msFeedCols.feedId().put( nrow-1, id ) ;
1221 msFeedCols.antennaId().put( nrow-1, 0 ) ;
1222 msFeedCols.numReceptors().put( nrow-1, numReceptors ) ;
1223 msFeedCols.polarizationType().put( nrow-1, polType ) ;
1224 msFeedCols.beamOffset().put( nrow-1, beamOffset ) ;
1225 msFeedCols.receptorAngle().put( nrow-1, receptorAngle ) ;
1226 msFeedCols.polResponse().put( nrow-1, polResponse ) ;
1227
1228// double endSec = gettimeofday_sec() ;
1229// os_ << "end MSWriter::addFeed() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1230}
1231
1232void MSWriter::addSpectralWindow( Int spwid, Int freqid )
1233{
1234// double startSec = gettimeofday_sec() ;
1235// os_ << "start MSWriter::addSpectralWindow() startSec=" << startSec << LogIO::POST ;
1236
1237 // add row
1238 MSSpectralWindow msSpw = mstable_->spectralWindow() ;
1239 while( (Int)msSpw.nrow() <= spwid ) {
1240 msSpw.addRow( 1, True ) ;
1241 }
1242
1243 MSSpWindowColumns msSpwCols( msSpw ) ;
1244
1245 STFrequencies stf = table_->frequencies() ;
1246
1247 // MEAS_FREQ_REF
1248 msSpwCols.measFreqRef().put( spwid, stf.getFrame( True ) ) ;
1249
1250 Double refpix ;
1251 Double refval ;
1252 Double inc ;
1253 stf.getEntry( refpix, refval, inc, (uInt)freqid ) ;
1254
1255 // NUM_CHAN
1256 Int nchan = (Int)(refpix * 2) + 1 ;
1257 if ( nchan == 0 )
1258 nchan = 1 ;
1259 msSpwCols.numChan().put( spwid, nchan ) ;
1260
1261 // TOTAL_BANDWIDTH
1262 Double bw = nchan * inc ;
1263 msSpwCols.totalBandwidth().put( spwid, bw ) ;
1264
1265 // REF_FREQUENCY
1266 Double refFreq = refval - refpix * inc ;
1267 msSpwCols.refFrequency().put( spwid, refFreq ) ;
1268
1269 // NET_SIDEBAND
1270 // tentative: USB->0, LSB->1
1271 Int netSideband = 0 ;
1272 if ( inc < 0 )
1273 netSideband = 1 ;
1274 msSpwCols.netSideband().put( spwid, netSideband ) ;
1275
1276 // RESOLUTION, CHAN_WIDTH, EFFECTIVE_BW
1277 Vector<Double> sharedDoubleArr( nchan, inc ) ;
1278 msSpwCols.resolution().put( spwid, sharedDoubleArr ) ;
1279 msSpwCols.chanWidth().put( spwid, sharedDoubleArr ) ;
1280 msSpwCols.effectiveBW().put( spwid, sharedDoubleArr ) ;
1281
1282 // CHAN_FREQ
1283 indgen( sharedDoubleArr, refFreq, inc ) ;
1284 msSpwCols.chanFreq().put( spwid, sharedDoubleArr ) ;
1285
1286// double endSec = gettimeofday_sec() ;
1287// os_ << "end MSWriter::addSpectralWindow() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1288}
1289
1290void MSWriter::addField( Int fid, String fieldname, String srcname, Double t, Vector<Double> rate )
1291{
1292// double startSec = gettimeofday_sec() ;
1293// os_ << "start MSWriter::addField() startSec=" << startSec << LogIO::POST ;
1294
1295 MSField msField = mstable_->field() ;
1296 while( (Int)msField.nrow() <= fid ) {
1297 msField.addRow( 1, True ) ;
1298 }
1299 MSFieldColumns msFieldCols( msField ) ;
1300
1301 // Access to SOURCE table
1302 MSSource msSrc = mstable_->source() ;
1303
1304 // fill target row
1305 msFieldCols.name().put( fid, fieldname ) ;
1306 msFieldCols.time().put( fid, t ) ;
1307 Int numPoly = 0 ;
1308 if ( anyNE( rate, 0.0 ) )
1309 numPoly = 1 ;
1310 msFieldCols.numPoly().put( fid, numPoly ) ;
1311 MSSourceIndex msSrcIdx( msSrc ) ;
1312 Int srcId = -1 ;
1313 Vector<Int> srcIdArr = msSrcIdx.matchSourceName( srcname ) ;
1314 if ( srcIdArr.size() != 0 ) {
1315 srcId = srcIdArr[0] ;
1316 MSSource msSrcSel = msSrc( msSrc.col("SOURCE_ID") == srcId ) ;
1317 ROMSSourceColumns msSrcCols( msSrcSel ) ;
1318 Vector<Double> srcDir = msSrcCols.direction()( 0 ) ;
1319 Matrix<Double> srcDirA( IPosition( 2, 2, 1+numPoly ) ) ;
1320// os_ << "srcDirA = " << srcDirA << LogIO::POST ;
1321// os_ << "sliced srcDirA = " << srcDirA.column( 0 ) << LogIO::POST ;
1322 srcDirA.column( 0 ) = srcDir ;
1323// os_ << "srcDirA = " << srcDirA << LogIO::POST ;
1324 if ( numPoly != 0 )
1325 srcDirA.column( 1 ) = rate ;
1326 msFieldCols.phaseDir().put( fid, srcDirA ) ;
1327 msFieldCols.referenceDir().put( fid, srcDirA ) ;
1328 msFieldCols.delayDir().put( fid, srcDirA ) ;
1329 }
1330 msFieldCols.sourceId().put( fid, srcId ) ;
1331
1332// double endSec = gettimeofday_sec() ;
1333// os_ << "end MSWriter::addField() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1334}
1335
1336void MSWriter::addPointing( String &name, Double &me, Double &interval, Matrix<Double> &dir )
1337{
1338// double startSec = gettimeofday_sec() ;
1339// os_ << "start MSWriter::addPointing() startSec=" << startSec << LogIO::POST ;
1340
1341 // access to POINTING subtable
1342 MSPointing msp = mstable_->pointing() ;
1343 uInt nrow = msp.nrow() ;
1344
1345 // add row
1346 msp.addRow( 1, True ) ;
1347
1348 // fill row
1349 TableRow row( msp ) ;
1350 TableRecord &rec = row.record() ;
1351 RecordFieldPtr<Int> antennaRF( rec, "ANTENNA_ID" ) ;
1352 *antennaRF = 0 ;
1353 RecordFieldPtr<Int> numpolyRF( rec, "NUM_POLY" ) ;
1354 *numpolyRF = dir.ncolumn() ;
1355 RecordFieldPtr<Double> timeRF( rec, "TIME" ) ;
1356 *timeRF = me ;
1357 RecordFieldPtr<Double> toriginRF( rec, "TIME_ORIGIN" ) ;
1358 *toriginRF = me ;
1359 RecordFieldPtr<Double> intervalRF( rec, "INTERVAL" ) ;
1360 *intervalRF = interval ;
1361 RecordFieldPtr<String> nameRF( rec, "NAME" ) ;
1362 *nameRF = name ;
1363 RecordFieldPtr<Bool> trackRF( rec, "TRACKING" ) ;
1364 *trackRF = True ;
1365 RecordFieldPtr< Array<Double> > dirRF( rec, "DIRECTION" ) ;
1366 *dirRF = dir ;
1367 RecordFieldPtr< Array<Double> > targetRF( rec, "TARGET" ) ;
1368 *dirRF = dir ;
1369 row.put( nrow ) ;
1370
1371// double endSec = gettimeofday_sec() ;
1372// os_ << "end MSWriter::addPointing() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1373}
1374
1375Int MSWriter::addPolarization( Vector<Int> polnos )
1376{
1377// double startSec = gettimeofday_sec() ;
1378// os_ << "start MSWriter::addPolarization() startSec=" << startSec << LogIO::POST ;
1379
1380// os_ << "polnos = " << polnos << LogIO::POST ;
1381 MSPolarization msPol = mstable_->polarization() ;
1382 uInt nrow = msPol.nrow() ;
1383
1384// // only 1 POLARIZATION row for 1 scantable
1385// if ( nrow > 0 )
1386// return 0 ;
1387
1388 Vector<Int> corrType = toCorrType( polnos ) ;
1389
1390 ROArrayColumn<Int> corrtCol( msPol, "CORR_TYPE" ) ;
1391 //Matrix<Int> corrTypeArr = corrtCol.getColumn() ;
1392 Int polid = -1 ;
1393 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
1394 Vector<Int> corrTypeArr = corrtCol( irow ) ;
1395 if ( corrType.nelements() == corrTypeArr.nelements()
1396 && allEQ( corrType, corrTypeArr ) ) {
1397 polid = irow ;
1398 break ;
1399 }
1400 }
1401
1402 if ( polid == -1 ) {
1403 MSPolarizationColumns msPolCols( msPol ) ;
1404
1405 // add row
1406 msPol.addRow( 1, True ) ;
1407 polid = (Int)nrow ;
1408
1409 // CORR_TYPE
1410 msPolCols.corrType().put( nrow, corrType ) ;
1411
1412 // NUM_CORR
1413 uInt npol = corrType.size() ;
1414 msPolCols.numCorr().put( nrow, npol ) ;
1415
1416 // CORR_PRODUCT
1417 Matrix<Int> corrProd( 2, npol, -1 ) ;
1418 if ( npol == 1 ) {
1419 corrProd = 0 ;
1420 }
1421 else if ( npol == 2 ) {
1422 corrProd.column( 0 ) = 0 ;
1423 corrProd.column( 1 ) = 1 ;
1424 }
1425 else {
1426 corrProd.column( 0 ) = 0 ;
1427 corrProd.column( 3 ) = 1 ;
1428 corrProd( 0,1 ) = 0 ;
1429 corrProd( 1,1 ) = 1 ;
1430 corrProd( 0,2 ) = 1 ;
1431 corrProd( 1,2 ) = 0 ;
1432 }
1433 msPolCols.corrProduct().put( nrow, corrProd ) ;
1434 }
1435
1436// double endSec = gettimeofday_sec() ;
1437// os_ << "end MSWriter::addPolarization() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1438
1439 return polid ;
1440}
1441
1442Int MSWriter::addDataDescription( Int polid, Int spwid )
1443{
1444// double startSec = gettimeofday_sec() ;
1445// os_ << "start MSWriter::addDataDescription() startSec=" << startSec << LogIO::POST ;
1446
1447 MSDataDescription msDataDesc = mstable_->dataDescription() ;
1448 uInt nrow = msDataDesc.nrow() ;
1449
1450 // only 1 POLARIZATION_ID for 1 scantable
1451 Int ddid = -1 ;
1452 ROScalarColumn<Int> spwCol( msDataDesc, "SPECTRAL_WINDOW_ID" ) ;
1453 Vector<Int> spwIds = spwCol.getColumn() ;
1454 //ROScalarColumn<Int> polCol( msDataDesc, "POLARIZATION_ID" ) ;
1455 //Vector<Int> polIds = polCol.getColumn() ;
1456 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
1457 //if ( spwid == spwIds[irow] && polid == polIds[irow] ) {
1458 if ( spwid == spwIds[irow] ) {
1459 ddid = irow ;
1460 break ;
1461 }
1462 }
1463// os_ << "ddid = " << ddid << LogIO::POST ;
1464
1465
1466 if ( ddid == -1 ) {
1467 msDataDesc.addRow( 1, True ) ;
1468 MSDataDescColumns msDataDescCols( msDataDesc ) ;
1469 msDataDescCols.polarizationId().put( nrow, polid ) ;
1470 msDataDescCols.spectralWindowId().put( nrow, spwid ) ;
1471 ddid = (Int)nrow ;
1472 }
1473
1474// double endSec = gettimeofday_sec() ;
1475// os_ << "end MSWriter::addDataDescription() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1476
1477 return ddid ;
1478}
1479
1480Int MSWriter::addState( Int st, Int &subscan )
1481{
1482// double startSec = gettimeofday_sec() ;
1483// os_ << "start MSWriter::addState() startSec=" << startSec << LogIO::POST ;
1484
1485 // access to STATE subtable
1486 MSState msState = mstable_->state() ;
1487 uInt nrow = msState.nrow() ;
1488
1489 String obsMode ;
1490 Bool isSignal ;
1491 Double tnoise ;
1492 Double tload ;
1493 queryType( st, obsMode, isSignal, tnoise, tload ) ;
1494// os_ << "obsMode = " << obsMode << " isSignal = " << isSignal << LogIO::POST ;
1495
1496 Int idx = -1 ;
1497 ROScalarColumn<String> obsModeCol( msState, "OBS_MODE" ) ;
1498 ROScalarColumn<Int> subscanCol( msState, "SUB_SCAN" ) ;
1499 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
1500 if ( obsModeCol(irow) == obsMode
1501 //&& sigCol(irow) == isSignal
1502 //&& refCol(irow) != isSignal
1503 && subscanCol(irow) == subscan ) {
1504 idx = irow ;
1505 break ;
1506 }
1507 }
1508 if ( idx == -1 ) {
1509 msState.addRow( 1, True ) ;
1510 TableRow row( msState ) ;
1511 TableRecord &rec = row.record() ;
1512 RecordFieldPtr<String> obsmodeRF( rec, "OBS_MODE" ) ;
1513 *obsmodeRF = obsMode ;
1514 RecordFieldPtr<Bool> sigRF( rec, "SIG" ) ;
1515 *sigRF = isSignal ;
1516 RecordFieldPtr<Bool> refRF( rec, "REF" ) ;
1517 *refRF = !isSignal ;
1518 RecordFieldPtr<Int> subscanRF( rec, "SUB_SCAN" ) ;
1519 *subscanRF = subscan ;
1520 RecordFieldPtr<Double> noiseRF( rec, "CAL" ) ;
1521 *noiseRF = tnoise ;
1522 RecordFieldPtr<Double> loadRF( rec, "LOAD" ) ;
1523 *loadRF = tload ;
1524 row.put( nrow ) ;
1525 idx = nrow ;
1526 }
1527 subscan++ ;
1528
1529// double endSec = gettimeofday_sec() ;
1530// os_ << "end MSWriter::addState() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1531
1532 return idx ;
1533}
1534
1535Vector<Int> MSWriter::toCorrType( Vector<Int> polnos )
1536{
1537// double startSec = gettimeofday_sec() ;
1538// os_ << "start MSWriter::toCorrType() startSec=" << startSec << LogIO::POST ;
1539
1540 uInt npol = polnos.size() ;
1541 Vector<Int> corrType( npol, Stokes::Undefined ) ;
1542
1543 if ( npol == 4 ) {
1544 if ( polType_ == "linear" ) {
1545 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
1546 if ( polnos[ipol] == 0 )
1547 corrType[ipol] = Stokes::XX ;
1548 else if ( polnos[ipol] == 1 )
1549 corrType[ipol] = Stokes::XY ;
1550 else if ( polnos[ipol] == 2 )
1551 corrType[ipol] = Stokes::YX ;
1552 else if ( polnos[ipol] == 3 )
1553 corrType[ipol] = Stokes::YY ;
1554 }
1555 }
1556 else if ( polType_ == "circular" ) {
1557 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
1558 if ( polnos[ipol] == 0 )
1559 corrType[ipol] = Stokes::RR ;
1560 else if ( polnos[ipol] == 1 )
1561 corrType[ipol] = Stokes::RL ;
1562 else if ( polnos[ipol] == 2 )
1563 corrType[ipol] = Stokes::LR ;
1564 else if ( polnos[ipol] == 3 )
1565 corrType[ipol] = Stokes::LL ;
1566 }
1567 }
1568 else if ( polType_ == "stokes" ) {
1569 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
1570 if ( polnos[ipol] == 0 )
1571 corrType[ipol] = Stokes::I ;
1572 else if ( polnos[ipol] == 1 )
1573 corrType[ipol] = Stokes::Q ;
1574 else if ( polnos[ipol] == 2 )
1575 corrType[ipol] = Stokes::U ;
1576 else if ( polnos[ipol] == 3 )
1577 corrType[ipol] = Stokes::V ;
1578 }
1579 }
1580 }
1581 else if ( npol == 2 ) {
1582 if ( polType_ == "linear" ) {
1583 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
1584 if ( polnos[ipol] == 0 )
1585 corrType[ipol] = Stokes::XX ;
1586 else if ( polnos[ipol] == 1 )
1587 corrType[ipol] = Stokes::YY ;
1588 }
1589 }
1590 else if ( polType_ == "circular" ) {
1591 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
1592 if ( polnos[ipol] == 0 )
1593 corrType[ipol] = Stokes::RR ;
1594 else if ( polnos[ipol] == 1 )
1595 corrType[ipol] = Stokes::LL ;
1596 }
1597 }
1598 else if ( polType_ == "stokes" ) {
1599 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
1600 if ( polnos[ipol] == 0 )
1601 corrType[ipol] = Stokes::I ;
1602 else if ( polnos[ipol] == 1 )
1603 corrType[ipol] = Stokes::V ;
1604 }
1605 }
1606 else if ( polType_ == "linpol" ) {
1607 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
1608 if ( polnos[ipol] == 1 )
1609 corrType[ipol] = Stokes::Plinear ;
1610 else if ( polnos[ipol] == 2 )
1611 corrType[ipol] = Stokes::Pangle ;
1612 }
1613 }
1614 }
1615 else if ( npol == 1 ) {
1616 if ( polType_ == "linear" )
1617 corrType[0] = Stokes::XX ;
1618 else if ( polType_ == "circular" )
1619 corrType[0] = Stokes::RR ;
1620 else if ( polType_ == "stokes" )
1621 corrType[0] = Stokes::I ;
1622 }
1623
1624// double endSec = gettimeofday_sec() ;
1625// os_ << "end MSWriter::toCorrType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1626
1627 return corrType ;
1628}
1629
1630void MSWriter::getValidTimeRange( Double &me, Double &interval, Table &tab )
1631{
1632// double startSec = gettimeofday_sec() ;
1633// os_ << "start MSWriter::getVaridTimeRange() startSec=" << startSec << LogIO::POST ;
1634
1635 // sort table
1636 //Table stab = tab.sort( "TIME" ) ;
1637
1638 ROScalarColumn<Double> timeCol( tab, "TIME" ) ;
1639 Vector<Double> timeArr = timeCol.getColumn() ;
1640 Double minTime ;
1641 Double maxTime ;
1642 minMax( minTime, maxTime, timeArr ) ;
1643 Double midTime = 0.5 * ( minTime + maxTime ) * 86400.0 ;
1644 // unit for TIME
1645 // Scantable: "d"
1646 // MS: "s"
1647 me = midTime ;
1648 interval = ( maxTime - minTime ) * 86400.0 ;
1649
1650// double endSec = gettimeofday_sec() ;
1651// os_ << "end MSWriter::getValidTimeRange() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1652}
1653
1654void MSWriter::getValidTimeRange( Double &me, Double &interval, Vector<Double> &atime, Vector<Double> &ainterval )
1655{
1656// double startSec = gettimeofday_sec() ;
1657// os_ << "start MSWriter::getVaridTimeRange() startSec=" << startSec << LogIO::POST ;
1658
1659 // sort table
1660 //Table stab = tab.sort( "TIME" ) ;
1661
1662 Double minTime ;
1663 Double maxTime ;
1664 minMax( minTime, maxTime, atime ) ;
1665 Double midTime = 0.5 * ( minTime + maxTime ) * 86400.0 ;
1666 // unit for TIME
1667 // Scantable: "d"
1668 // MS: "s"
1669 me = midTime ;
1670 interval = ( maxTime - minTime ) * 86400.0 + mean( ainterval ) ;
1671
1672// double endSec = gettimeofday_sec() ;
1673// os_ << "end MSWriter::getValidTimeRange() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1674}
1675
1676//void MSWriter::queryType( Int type, String &stype, Bool &b )
1677void MSWriter::queryType( Int type, String &stype, Bool &b, Double &t, Double &l )
1678{
1679// double startSec = gettimeofday_sec() ;
1680// os_ << "start MSWriter::queryType() startSec=" << startSec << LogIO::POST ;
1681
1682 // 2011/03/14 TN
1683 // OBS_MODE string of MS created by importasdm task is slightly
1684 // (but critically) changed.
1685 switch ( type ) {
1686 case SrcType::PSON:
1687 stype = "OBSERVE_TARGET_ON_SOURCE,POSITION_SWITCH" ;
1688 b = True ;
1689 t = 0.0 ;
1690 l = 0.0 ;
1691 break ;
1692 case SrcType::PSOFF:
1693 stype = "OBSERVE_TARGET_OFF_SOURCE,POSITION_SWITCH" ;
1694 b = False ;
1695 t = 0.0 ;
1696 l = 0.0 ;
1697 break ;
1698 case SrcType::NOD:
1699 stype = "OBSERVE_TARGET_ON_SOURCE,NOD" ;
1700 b = True ;
1701 t = 0.0 ;
1702 l = 0.0 ;
1703 break ;
1704 case SrcType::FSON:
1705 stype = "OBSERVE_TARGET_ON_SOURCE,FREQUENCY_SWITCH_SIG" ;
1706 b = True ;
1707 t = 0.0 ;
1708 l = 0.0 ;
1709 break ;
1710 case SrcType::FSOFF:
1711 stype = "OBSERVE_TARGET_ON_SOURCE,FREQUENCY_SWITCH_REF" ;
1712 b = False ;
1713 t = 0.0 ;
1714 l = 0.0 ;
1715 break ;
1716 case SrcType::SKY:
1717 stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,UNSPECIFIED" ;
1718 b = False ;
1719 t = 0.0 ;
1720 l = 1.0 ;
1721 break ;
1722 case SrcType::HOT:
1723 stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,UNSPECIFIED" ;
1724 b = False ;
1725 t = 0.0 ;
1726 l = 1.0 ;
1727 break ;
1728 case SrcType::WARM:
1729 stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,UNSPECIFIED" ;
1730 t = 0.0 ;
1731 b = False ;
1732 l = 1.0 ;
1733 break ;
1734 case SrcType::COLD:
1735 stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,UNSPECIFIED" ;
1736 b = False ;
1737 t = 0.0 ;
1738 l = 1.0 ;
1739 break ;
1740 case SrcType::PONCAL:
1741 stype = "CALIBRATE_TEMPERATURE_ON_SOURCE,POSITION_SWITCH" ;
1742 b = True ;
1743 t = 1.0 ;
1744 l = 0.0 ;
1745 break ;
1746 case SrcType::POFFCAL:
1747 stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,POSITION_SWITCH" ;
1748 b = False ;
1749 t = 1.0 ;
1750 l = 0.0 ;
1751 break ;
1752 case SrcType::NODCAL:
1753 stype = "CALIBRATE_TEMPERATURE_ON_SOURCE,NOD" ;
1754 b = True ;
1755 t = 1.0 ;
1756 l = 0.0 ;
1757 break ;
1758 case SrcType::FONCAL:
1759 stype = "CALIBRATE_TEMPERATURE_ON_SOURCE,FREQUENCY_SWITCH_SIG" ;
1760 b = True ;
1761 t = 1.0 ;
1762 l = 0.0 ;
1763 break ;
1764 case SrcType::FOFFCAL:
1765 stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,FREQUENCY_SWITCH_REF" ;
1766 b = False ;
1767 t = 1.0 ;
1768 l = 0.0 ;
1769 break ;
1770 case SrcType::FSLO:
1771 stype = "OBSERVE_TARGET_ON_SOURCE,FREQUENCY_SWITCH_LOWER" ;
1772 b = True ;
1773 t = 0.0 ;
1774 l = 0.0 ;
1775 break ;
1776 case SrcType::FLOOFF:
1777 stype = "OBSERVE_TARGET_OFF_SOURCE,FREQUENCY_SWITCH_LOWER" ;
1778 b = False ;
1779 t = 0.0 ;
1780 l = 0.0 ;
1781 break ;
1782 case SrcType::FLOSKY:
1783 stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,FREQUENCY_SWITCH_LOWER" ;
1784 b = False ;
1785 t = 0.0 ;
1786 l = 1.0 ;
1787 break ;
1788 case SrcType::FLOHOT:
1789 stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,FREQUENCY_SWITCH_LOWER" ;
1790 b = False ;
1791 t = 0.0 ;
1792 l = 1.0 ;
1793 break ;
1794 case SrcType::FLOWARM:
1795 stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,FREQUENCY_SWITCH_LOWER" ;
1796 b = False ;
1797 t = 0.0 ;
1798 l = 1.0 ;
1799 break ;
1800 case SrcType::FLOCOLD:
1801 stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,FREQUENCY_SWITCH_LOWER" ;
1802 b = False ;
1803 t = 0.0 ;
1804 l = 1.0 ;
1805 break ;
1806 case SrcType::FSHI:
1807 stype = "OBSERVE_TARGET_ON_SOURCE,FREQUENCY_SWITCH_HIGHER" ;
1808 b = True ;
1809 t = 0.0 ;
1810 l = 0.0 ;
1811 break ;
1812 case SrcType::FHIOFF:
1813 stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,FREQUENCY_SWITCH_HIGHER" ;
1814 b = False ;
1815 t = 0.0 ;
1816 l = 0.0 ;
1817 break ;
1818 case SrcType::FHISKY:
1819 stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,FREQUENCY_SWITCH_HIGHER" ;
1820 b = False ;
1821 t = 0.0 ;
1822 l = 1.0 ;
1823 break ;
1824 case SrcType::FHIHOT:
1825 stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,FREQUENCY_SWITCH_HIGHER" ;
1826 b = False ;
1827 t = 0.0 ;
1828 l = 1.0 ;
1829 break ;
1830 case SrcType::FHIWARM:
1831 stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,FREQUENCY_SWITCH_HIGHER" ;
1832 b = False ;
1833 t = 0.0 ;
1834 l = 1.0 ;
1835 break ;
1836 case SrcType::FHICOLD:
1837 stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,FREQUENCY_SWITCH_HIGHER" ;
1838 b = False ;
1839 t = 0.0 ;
1840 l = 1.0 ;
1841 break ;
1842 case SrcType::SIG:
1843 stype = "OBSERVE_TARGET_ON_SOURCE,UNSPECIFIED" ;
1844 b = True ;
1845 t = 0.0 ;
1846 l = 0.0 ;
1847 break ;
1848 case SrcType::REF:
1849 stype = "OBSERVE_TARGET_ON_SOURCE,UNSPECIFIED" ;
1850 b = False ;
1851 t = 0.0 ;
1852 l = 0.0 ;
1853 break ;
1854 default:
1855 stype = "UNSPECIFIED" ;
1856 b = True ;
1857 t = 0.0 ;
1858 l = 0.0 ;
1859 break ;
1860 }
1861
1862// double endSec = gettimeofday_sec() ;
1863// os_ << "end MSWriter::queryType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1864}
1865
1866Double MSWriter::getDishDiameter( String antname )
1867{
1868 Double diameter = 0.0 ;
1869
1870 antname.upcase() ;
1871
1872 if ( antname.matches( Regex( "DV[0-9]+$" ) )
1873 || antname.matches( Regex( "DA[0-9]+$" ) )
1874 || antname.matches( Regex( "PM[0-9]+$" ) ) )
1875 diameter = 12.0 ;
1876 else if ( antname.matches( Regex( "CM[0-9]+$" ) ) )
1877 diameter = 7.0 ;
1878 else if ( antname.contains( "GBT" ) )
1879 diameter = 104.9 ;
1880 else if ( antname.contains( "MOPRA" ) )
1881 diameter = 22.0 ;
1882 else if ( antname.contains( "PKS" ) || antname.contains( "PARKS" ) )
1883 diameter = 64.0 ;
1884 else if ( antname.contains( "TIDBINBILLA" ) )
1885 diameter = 70.0 ;
1886 else if ( antname.contains( "CEDUNA" ) )
1887 diameter = 30.0 ;
1888 else if ( antname.contains( "HOBART" ) )
1889 diameter = 26.0 ;
1890 else if ( antname.contains( "APEX" ) )
1891 diameter = 12.0 ;
1892 else if ( antname.contains( "ASTE" ) )
1893 diameter = 10.0 ;
1894 else if ( antname.contains( "NRO" ) )
1895 diameter = 45.0 ;
1896 else
1897 diameter = 1.0 ;
1898
1899 return diameter ;
1900}
1901
1902void MSWriter::infillSpectralWindow()
1903{
1904 MSSpectralWindow msSpw = mstable_->spectralWindow() ;
1905 MSSpWindowColumns msSpwCols( msSpw ) ;
1906 uInt nrow = msSpw.nrow() ;
1907
1908 ScalarColumn<Int> measFreqRefCol = msSpwCols.measFreqRef() ;
1909 ArrayColumn<Double> chanFreqCol = msSpwCols.chanFreq() ;
1910 ArrayColumn<Double> chanWidthCol = msSpwCols.chanWidth() ;
1911 ArrayColumn<Double> effectiveBWCol = msSpwCols.effectiveBW() ;
1912 ArrayColumn<Double> resolutionCol = msSpwCols.resolution() ;
1913 Vector<Double> dummy( 1, 0.0 ) ;
1914 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
1915 if ( !(chanFreqCol.isDefined( irow )) ) {
1916 measFreqRefCol.put( irow, 1 ) ;
1917 chanFreqCol.put( irow, dummy ) ;
1918 chanWidthCol.put( irow, dummy ) ;
1919 effectiveBWCol.put( irow, dummy ) ;
1920 resolutionCol.put( irow, dummy ) ;
1921 }
1922 }
1923
1924}
1925
1926}
Note: See TracBrowser for help on using the repository browser.