source: trunk/src/MSWriter.cpp@ 1975

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

New Development: Yes

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...

Writer improvement.
File size: 31.5 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
20#include <tables/Tables/ExprNode.h>
21#include <tables/Tables/TableDesc.h>
22#include <tables/Tables/SetupNewTab.h>
23#include <tables/Tables/TableIter.h>
24#include <tables/Tables/RefRows.h>
25
26#include <ms/MeasurementSets/MeasurementSet.h>
27#include <ms/MeasurementSets/MSColumns.h>
28#include <ms/MeasurementSets/MSPolIndex.h>
29#include <ms/MeasurementSets/MSDataDescIndex.h>
30#include <ms/MeasurementSets/MSSourceIndex.h>
31
32#include "MSWriter.h"
33#include "STHeader.h"
34#include "STFrequencies.h"
35#include "STMolecules.h"
36
37using namespace casa ;
38
39namespace asap
40{
41
42MSWriter::MSWriter(CountedPtr<Scantable> stable)
43 : table_(stable)
44{
45 os_ = LogIO() ;
46 os_.origin( LogOrigin( "MSWriter", "MSWriter()", WHERE ) ) ;
47 os_ << "MSWriter::MSWriter()" << LogIO::POST ;
48
49 // initialize writer
50 init() ;
51}
52
53MSWriter::~MSWriter()
54{
55 os_.origin( LogOrigin( "MSWriter", "~MSWriter()", WHERE ) ) ;
56 os_ << "MSWriter::~MSWriter()" << LogIO::POST ;
57}
58
59bool MSWriter::write(const string& filename, const Record& rec)
60{
61 os_.origin( LogOrigin( "MSWriter", "write()", WHERE ) ) ;
62 os_ << "MSWriter::write()" << LogIO::POST ;
63
64 filename_ = filename ;
65
66 // parsing MS options
67 Bool overwrite = False ;
68 if ( rec.isDefined( "ms" ) ) {
69 Record msrec = rec.asRecord( "ms" ) ;
70 if ( msrec.isDefined( "overwrite" ) ) {
71 overwrite = msrec.asBool( "overwrite" ) ;
72 }
73 }
74
75 os_ << "Parsing MS options" << endl ;
76 os_ << " overwrite = " << overwrite << LogIO::POST ;
77
78 File file( filename_ ) ;
79 if ( file.exists() ) {
80 if ( overwrite ) {
81 os_ << filename_ << " exists. Overwrite existing data... " << LogIO::POST ;
82 if ( file.isRegular() ) RegularFile(file).remove() ;
83 else if ( file.isDirectory() ) Directory(file).removeRecursive() ;
84 else SymLink(file).remove() ;
85 }
86 else {
87 os_ << LogIO::SEVERE << "ERROR: " << filename_ << " exists..." << LogIO::POST ;
88 return False ;
89 }
90 }
91
92 // set up MS
93 setupMS() ;
94
95 // subtables
96 // OBSERVATION
97 fillObservation() ;
98
99 // ANTENNA
100 fillAntenna() ;
101
102 // PROCESSOR
103 fillProcessor() ;
104
105 // SOURCE
106 fillSource() ;
107
108 // MAIN
109 // Iterate over several ids
110 Vector<uInt> processedFreqId( 0 ) ;
111 Int defaultFieldId = 0 ;
112 //
113 // ITERATION: FIELDNAME
114 //
115 Int added0 = 0 ;
116 Int current0 = mstable_->nrow() ;
117 TableIterator iter0( table_->table(), "FIELDNAME" ) ;
118 while( !iter0.pastEnd() ) {
119 Table t0( iter0.table() ) ;
120 ROScalarColumn<String> sharedStrCol( t0, "FIELDNAME" ) ;
121 String fieldName = sharedStrCol( 0 ) ;
122 sharedStrCol.attach( t0, "SRCNAME" ) ;
123 String srcName = sharedStrCol( 0 ) ;
124 ROScalarColumn<Double> timeCol( t0, "TIME" ) ;
125 Double minTime = min( timeCol.getColumn() ) ;
126 ROArrayColumn<Double> scanRateCol( t0, "SCANRATE" ) ;
127 Vector<Double> scanRate = scanRateCol.getColumn()[0] ;
128 String::size_type pos = fieldName.find( "__" ) ;
129 Int fieldId = -1 ;
130 if ( pos != String::npos ) {
131 os_ << "fieldName.substr( pos+2 )=" << fieldName.substr( pos+2 ) << LogIO::POST ;
132 fieldId = String::toInt( fieldName.substr( pos+2 ) ) ;
133 fieldName = fieldName.substr( 0, pos ) ;
134 }
135 else {
136 os_ << "use default field id" << LogIO::POST ;
137 fieldId = defaultFieldId ;
138 defaultFieldId++ ;
139 }
140 os_ << "fieldId" << fieldId << ": " << fieldName << LogIO::POST ;
141 //
142 // ITERATION: BEAMNO
143 //
144 Int added1 = 0 ;
145 Int current1 = mstable_->nrow() ;
146 TableIterator iter1( t0, "BEAMNO" ) ;
147 while( !iter1.pastEnd() ) {
148 Table t1( iter1.table() ) ;
149 ROScalarColumn<uInt> beamNoCol( t1, "BEAMNO" ) ;
150 uInt beamNo = beamNoCol( 0 ) ;
151 os_ << "beamNo = " << beamNo << LogIO::POST ;
152 //
153 // ITERATION: SCANNO
154 //
155 Int added2 = 0 ;
156 Int current2 = mstable_->nrow() ;
157 TableIterator iter2( t1, "SCANNO" ) ;
158 while( !iter2.pastEnd() ) {
159 Table t2( iter2.table() ) ;
160 ROScalarColumn<uInt> scanNoCol( t2, "SCANNO" ) ;
161 uInt scanNo = scanNoCol( 0 ) ;
162 os_ << "scanNo = " << scanNo << LogIO::POST ;
163 //
164 // ITERATION: CYCLENO
165 //
166 Int added3 = 0 ;
167 Int current3 = mstable_->nrow() ;
168 TableIterator iter3( t2, "CYCLENO" ) ;
169 while( !iter3.pastEnd() ) {
170 Table t3( iter3.table() ) ;
171 //
172 // ITERATION: IFNO
173 //
174 Int added4 = 0 ;
175 Int current4 = mstable_->nrow() ;
176 TableIterator iter4( t3, "IFNO" ) ;
177 while( !iter4.pastEnd() ) {
178 Table t4( iter4.table() ) ;
179 ROScalarColumn<uInt> ifNoCol( t4, "IFNO" ) ;
180 uInt ifNo = ifNoCol( 0 ) ;
181 os_ << "ifNo = " << ifNo << LogIO::POST ;
182 ROScalarColumn<uInt> freqIdCol( t4, "FREQ_ID" ) ;
183 uInt freqId = freqIdCol( 0 ) ;
184 os_ << "freqId = " << freqId << LogIO::POST ;
185 //
186 // ITERATION: TIME
187 //
188 Int added5 = 0 ;
189 Int current5 = mstable_->nrow() ;
190 TableIterator iter5( t4, "TIME" ) ;
191 while( !iter5.pastEnd() ) {
192 Table t5( iter5.table().sort("POLNO") ) ;
193 Int prevnr = mstable_->nrow() ;
194 Int nrow = t5.nrow() ;
195 os_ << "nrow = " << nrow << LogIO::POST ;
196
197 // add row
198 mstable_->addRow( 1, True ) ;
199
200 Vector<Int> polnos( nrow ) ;
201 indgen( polnos, 0 ) ;
202 Int polid = addPolarization( polnos ) ;
203 os_ << "polid = " << polid << LogIO::POST ;
204 //
205 // LOOP: POLNO
206 //
207 for ( Int ipol = 0 ; ipol < nrow ; ipol++ ) {
208 }
209
210 // TIME and TIME_CENTROID
211 ROScalarMeasColumn<MEpoch> timeCol( t5, "TIME" ) ;
212 ScalarMeasColumn<MEpoch> msTimeCol( *mstable_, "TIME" ) ;
213 msTimeCol.put( prevnr, timeCol( 0 ) ) ;
214 msTimeCol.attach( *mstable_, "TIME_CENTROID" ) ;
215 msTimeCol.put( prevnr, timeCol( 0 ) ) ;
216
217 // INTERVAL and EXPOSURE
218 ROScalarColumn<Double> intervalCol( t5, "INTERVAL" ) ;
219 ScalarColumn<Double> msIntervalCol( *mstable_, "INTERVAL" ) ;
220 msIntervalCol.put( prevnr, intervalCol( 0 ) ) ;
221 msIntervalCol.attach( *mstable_, "EXPOSURE" ) ;
222 msIntervalCol.put( prevnr, intervalCol( 0 ) ) ;
223
224 // add DATA_DESCRIPTION row
225 Int ddid = addDataDescription( polid, ifNo ) ;
226 os_ << "ddid = " << ddid << LogIO::POST ;
227 ScalarColumn<Int> ddIdCol( *mstable_, "DATA_DESC_ID" ) ;
228 ddIdCol.put( prevnr, ddid ) ;
229
230 added5 += 1 ;
231 os_ << "added5 = " << added5 << " current5 = " << current5 << LogIO::POST ;
232 iter5.next() ;
233 }
234
235 // add SPECTRAL_WINDOW row
236 if ( allNE( processedFreqId, freqId ) ) {
237 uInt vsize = processedFreqId.size() ;
238 processedFreqId.resize( vsize+1, True ) ;
239 processedFreqId[vsize] = freqId ;
240 addSpectralWindow( ifNo, freqId ) ;
241 }
242
243 added4 += added5 ;
244 os_ << "added4 = " << added4 << " current4 = " << current4 << LogIO::POST ;
245 iter4.next() ;
246 }
247 added3 += added4 ;
248 os_ << "added3 = " << added3 << " current3 = " << current3 << LogIO::POST ;
249 iter3.next() ;
250 }
251
252 // SCAN_NUMBER
253 RefRows rows3( current3, current3+added3-1 ) ;
254 Vector<Int> scanNum( added3, scanNo ) ;
255 ScalarColumn<Int> scanNumCol( *mstable_, "SCAN_NUMBER" ) ;
256 scanNumCol.putColumnCells( rows3, scanNum ) ;
257
258 added2 += added3 ;
259 os_ << "added2 = " << added2 << " current2 = " << current2 << LogIO::POST ;
260 iter2.next() ;
261 }
262
263 // FEED1 and FEED2
264 RefRows rows2( current2, current2+added2-1 ) ;
265 Vector<Int> feedId( added2, beamNo ) ;
266 ScalarColumn<Int> feedCol( *mstable_, "FEED1" ) ;
267 feedCol.putColumnCells( rows2, feedId ) ;
268 feedCol.attach( *mstable_, "FEED2" ) ;
269 feedCol.putColumnCells( rows2, feedId ) ;
270
271 // add FEED row
272 addFeed( beamNo ) ;
273
274 added1 += added2 ;
275 os_ << "added1 = " << added1 << " current1 = " << current1 << LogIO::POST ;
276 iter1.next() ;
277 }
278
279 // FIELD_ID
280 RefRows rows1( current1, current1+added1-1 ) ;
281 Vector<Int> fieldIds( added1, fieldId ) ;
282 ScalarColumn<Int> fieldIdCol( *mstable_, "FIELD_ID" ) ;
283 fieldIdCol.putColumnCells( rows1, fieldIds ) ;
284
285 // add FIELD row
286 addField( fieldId, fieldName, srcName, minTime, scanRate ) ;
287
288 added0 += added1 ;
289 os_ << "added0 = " << added0 << " current0 = " << current0 << LogIO::POST ;
290 iter0.next() ;
291 }
292
293 // OBSERVATION_ID is always 0
294 ScalarColumn<Int> sharedIntCol( *mstable_, "OBSERVATION_ID" ) ;
295 Vector<Int> sharedIntArr( added0, 0 ) ;
296 sharedIntCol.putColumn( sharedIntArr ) ;
297
298 // ANTENNA1 and ANTENNA2 are always 0
299 sharedIntArr = 0 ;
300 sharedIntCol.attach( *mstable_, "ANTENNA1" ) ;
301 sharedIntCol.putColumn( sharedIntArr ) ;
302 sharedIntCol.attach( *mstable_, "ANTENNA2" ) ;
303 sharedIntCol.putColumn( sharedIntArr ) ;
304
305 // ARRAY_ID is tentatively set to 0
306 sharedIntArr = 0 ;
307 sharedIntCol.attach( *mstable_, "ARRAY_ID" ) ;
308 sharedIntCol.putColumn( sharedIntArr ) ;
309
310 // PROCESSOR_ID is tentatively set to 0
311 sharedIntArr = 0 ;
312 sharedIntCol.attach( *mstable_, "PROCESSOR_ID" ) ;
313 sharedIntCol.putColumn( sharedIntArr ) ;
314
315 return True ;
316}
317
318void MSWriter::init()
319{
320// os_.origin( LogOrigin( "MSWriter", "init()", WHERE ) ) ;
321 os_ << "MSWriter::init()" << LogIO::POST ;
322
323 // access to scantable
324 header_ = table_->getHeader() ;
325
326 // FLOAT_DATA? or DATA?
327 if ( header_.npol > 2 ) {
328 isFloatData_ = False ;
329 isData_ = True ;
330 }
331 else {
332 isFloatData_ = True ;
333 isData_ = False ;
334 }
335
336 // polarization type
337 polType_ = header_.poltype ;
338 if ( polType_ == "" )
339 polType_ = "stokes" ;
340 else if ( polType_.find( "linear" ) != String::npos )
341 polType_ = "linear" ;
342 else if ( polType_.find( "circular" ) != String::npos )
343 polType_ = "circular" ;
344 else if ( polType_.find( "stokes" ) != String::npos )
345 polType_ = "stokes" ;
346 else if ( polType_.find( "linpol" ) != String::npos )
347 polType_ = "linpol" ;
348 else
349 polType_ = "notype" ;
350
351}
352
353void MSWriter::setupMS()
354{
355// os_.origin( LogOrigin( "MSWriter", "setupMS()", WHERE ) ) ;
356 os_ << "MSWriter::setupMS()" << LogIO::POST ;
357
358 TableDesc msDesc = MeasurementSet::requiredTableDesc() ;
359
360 // any additional columns?
361 // TODO: add FLOAT_DATA column if npol == 2 otherwise add DATA column
362
363 SetupNewTable newtab( filename_, msDesc, Table::New ) ;
364
365 mstable_ = new MeasurementSet( newtab ) ;
366
367 // create subtables
368 TableDesc antennaDesc = MSAntenna::requiredTableDesc() ;
369 SetupNewTable antennaTab( mstable_->antennaTableName(), antennaDesc, Table::New ) ;
370 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::ANTENNA ), Table( antennaTab ) ) ;
371
372 TableDesc dataDescDesc = MSDataDescription::requiredTableDesc() ;
373 SetupNewTable dataDescTab( mstable_->dataDescriptionTableName(), dataDescDesc, Table::New ) ;
374 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::DATA_DESCRIPTION ), Table( dataDescTab ) ) ;
375
376 TableDesc dopplerDesc = MSDoppler::requiredTableDesc() ;
377 SetupNewTable dopplerTab( mstable_->dopplerTableName(), dopplerDesc, Table::New ) ;
378 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::DOPPLER ), Table( dopplerTab ) ) ;
379
380 TableDesc feedDesc = MSFeed::requiredTableDesc() ;
381 SetupNewTable feedTab( mstable_->feedTableName(), feedDesc, Table::New ) ;
382 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FEED ), Table( feedTab ) ) ;
383
384 TableDesc fieldDesc = MSField::requiredTableDesc() ;
385 SetupNewTable fieldTab( mstable_->fieldTableName(), fieldDesc, Table::New ) ;
386 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FIELD ), Table( fieldTab ) ) ;
387
388 TableDesc flagCmdDesc = MSFlagCmd::requiredTableDesc() ;
389 SetupNewTable flagCmdTab( mstable_->flagCmdTableName(), flagCmdDesc, Table::New ) ;
390 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FLAG_CMD ), Table( flagCmdTab ) ) ;
391
392 TableDesc freqOffsetDesc = MSFreqOffset::requiredTableDesc() ;
393 SetupNewTable freqOffsetTab( mstable_->freqOffsetTableName(), freqOffsetDesc, Table::New ) ;
394 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FREQ_OFFSET ), Table( freqOffsetTab ) ) ;
395
396 TableDesc historyDesc = MSHistory::requiredTableDesc() ;
397 SetupNewTable historyTab( mstable_->historyTableName(), historyDesc, Table::New ) ;
398 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::HISTORY ), Table( historyTab ) ) ;
399
400 TableDesc observationDesc = MSObservation::requiredTableDesc() ;
401 SetupNewTable observationTab( mstable_->observationTableName(), observationDesc, Table::New ) ;
402 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::OBSERVATION ), Table( observationTab ) ) ;
403
404 TableDesc pointingDesc = MSPointing::requiredTableDesc() ;
405 SetupNewTable pointingTab( mstable_->pointingTableName(), pointingDesc, Table::New ) ;
406 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::POINTING ), Table( pointingTab ) ) ;
407
408 TableDesc polarizationDesc = MSPolarization::requiredTableDesc() ;
409 SetupNewTable polarizationTab( mstable_->polarizationTableName(), polarizationDesc, Table::New ) ;
410 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::POLARIZATION ), Table( polarizationTab ) ) ;
411
412 TableDesc processorDesc = MSProcessor::requiredTableDesc() ;
413 SetupNewTable processorTab( mstable_->processorTableName(), processorDesc, Table::New ) ;
414 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::PROCESSOR ), Table( processorTab ) ) ;
415
416 TableDesc sourceDesc = MSSource::requiredTableDesc() ;
417 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::TRANSITION, 1 ) ;
418 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::REST_FREQUENCY, 1 ) ;
419 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::SYSVEL, 1 ) ;
420 SetupNewTable sourceTab( mstable_->sourceTableName(), sourceDesc, Table::New ) ;
421 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SOURCE ), Table( sourceTab ) ) ;
422
423 TableDesc spwDesc = MSSpectralWindow::requiredTableDesc() ;
424 SetupNewTable spwTab( mstable_->spectralWindowTableName(), spwDesc, Table::New ) ;
425 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SPECTRAL_WINDOW ), Table( spwTab ) ) ;
426
427 TableDesc stateDesc = MSState::requiredTableDesc() ;
428 SetupNewTable stateTab( mstable_->stateTableName(), stateDesc, Table::New ) ;
429 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::STATE ), Table( stateTab ) ) ;
430
431 // TODO: add TCAL_SPECTRUM and TSYS_SPECTRUM if necessary
432 TableDesc sysCalDesc = MSSysCal::requiredTableDesc() ;
433 SetupNewTable sysCalTab( mstable_->sysCalTableName(), sysCalDesc, Table::New ) ;
434 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SYSCAL ), Table( sysCalTab ) ) ;
435
436 TableDesc weatherDesc = MSWeather::requiredTableDesc() ;
437 SetupNewTable weatherTab( mstable_->weatherTableName(), weatherDesc, Table::New ) ;
438 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::WEATHER ), Table( weatherTab ) ) ;
439
440 mstable_->initRefs() ;
441
442}
443
444void MSWriter::fillObservation()
445{
446 os_ << "set up OBSERVATION subtable" << LogIO::POST ;
447 // only 1 row
448 mstable_->observation().addRow( 1, True ) ;
449 MSObservationColumns msObsCols( mstable_->observation() ) ;
450 msObsCols.observer().put( 0, header_.observer ) ;
451 // tentatively put antennaname (from ANTENNA subtable)
452 String hAntennaName = header_.antennaname ;
453 String::size_type pos = hAntennaName.find( "//" ) ;
454 String telescopeName ;
455 if ( pos != String::npos ) {
456 telescopeName = hAntennaName.substr( 0, pos ) ;
457 }
458 else {
459 pos = hAntennaName.find( "@" ) ;
460 telescopeName = hAntennaName.substr( 0, pos ) ;
461 }
462 os_ << "telescopeName = " << telescopeName << LogIO::POST ;
463 msObsCols.telescopeName().put( 0, telescopeName ) ;
464 msObsCols.project().put( 0, header_.project ) ;
465 //ScalarMeasColumn<MEpoch> timeCol( table_->table().sort("TIME"), "TIME" ) ;
466 Table sortedtable = table_->table().sort("TIME") ;
467 ScalarMeasColumn<MEpoch> timeCol( sortedtable, "TIME" ) ;
468 Vector<MEpoch> trange( 2 ) ;
469 trange[0] = timeCol( 0 ) ;
470 trange[1] = timeCol( table_->nrow()-1 ) ;
471 msObsCols.timeRangeMeas().put( 0, trange ) ;
472}
473
474void MSWriter::fillAntenna()
475{
476 os_ << "set up ANTENNA subtable" << LogIO::POST ;
477 // only 1 row
478 mstable_->antenna().addRow( 1, True ) ;
479 MSAntennaColumns msAntCols( mstable_->antenna() ) ;
480
481 String hAntennaName = header_.antennaname ;
482 String::size_type pos = hAntennaName.find( "//" ) ;
483 String antennaName ;
484 String stationName ;
485 if ( pos != String::npos ) {
486 hAntennaName = hAntennaName.substr( pos+2 ) ;
487 }
488 pos = hAntennaName.find( "@" ) ;
489 if ( pos != String::npos ) {
490 antennaName = hAntennaName.substr( 0, pos ) ;
491 stationName = hAntennaName.substr( pos+1 ) ;
492 }
493 else {
494 antennaName = hAntennaName ;
495 stationName = hAntennaName ;
496 }
497 os_ << "antennaName = " << antennaName << LogIO::POST ;
498 os_ << "stationName = " << stationName << LogIO::POST ;
499
500 msAntCols.name().put( 0, antennaName ) ;
501 msAntCols.station().put( 0, stationName ) ;
502
503 os_ << "antennaPosition = " << header_.antennaposition << LogIO::POST ;
504
505 msAntCols.position().put( 0, header_.antennaposition ) ;
506}
507
508void MSWriter::fillProcessor()
509{
510 os_ << "set up PROCESSOR subtable" << LogIO::POST ;
511
512 // only add empty 1 row
513 MSProcessor msProc = mstable_->processor() ;
514 msProc.addRow( 1, True ) ;
515}
516
517void MSWriter::fillSource()
518{
519 os_ << "set up SOURCE subtable" << LogIO::POST ;
520
521
522 // access to MS SOURCE subtable
523 MSSource msSrc = mstable_->source() ;
524
525 // access to MOLECULE subtable
526 STMolecules stm = table_->molecules() ;
527
528 Int srcId = 0 ;
529 Vector<Double> restFreq ;
530 Vector<String> molName ;
531 Vector<String> fMolName ;
532 //
533 // ITERATION: SRCNAME
534 //
535 Int added0 = 0 ;
536 Int current0 = msSrc.nrow() ;
537 TableIterator iter0( table_->table(), "SRCNAME" ) ;
538 while( !iter0.pastEnd() ) {
539 Table t0( iter0.table() ) ;
540
541 // get necessary information
542 ROScalarColumn<String> srcNameCol( t0, "SRCNAME" ) ;
543 String srcName = srcNameCol( 0 ) ;
544 ROArrayColumn<Double> sharedDArrRCol( t0, "SRCPROPERMOTION" ) ;
545 Vector<Double> srcPM = sharedDArrRCol( 0 ) ;
546 sharedDArrRCol.attach( t0, "SRCDIRECTION" ) ;
547 Vector<Double> srcDir = sharedDArrRCol( 0 ) ;
548 ROScalarColumn<Double> srcVelCol( t0, "SRCVELOCITY" ) ;
549 Double srcVel = srcVelCol( 0 ) ;
550
551 //
552 // ITERATION: MOLECULE_ID
553 //
554 Int added1 = 0 ;
555 Int current1 = msSrc.nrow() ;
556 TableIterator iter1( t0, "MOLECULE_ID" ) ;
557 while( !iter1.pastEnd() ) {
558 Table t1( iter1.table() ) ;
559
560 // get necessary information
561 ROScalarColumn<uInt> molIdCol( t1, "MOLECULE_ID" ) ;
562 uInt molId = molIdCol( 0 ) ;
563 stm.getEntry( restFreq, molName, fMolName, molId ) ;
564
565 //
566 // ITERATION: IFNO
567 //
568 Int added2 = 0 ;
569 Int current2 = msSrc.nrow() ;
570 TableIterator iter2( t1, "IFNO" ) ;
571 while( !iter2.pastEnd() ) {
572 Table t2( iter2.table() ) ;
573 uInt prevnr = msSrc.nrow() ;
574
575 // get necessary information
576 ROScalarColumn<uInt> ifNoCol( t2, "IFNO" ) ;
577 uInt ifno = ifNoCol( 0 ) ; // IFNO = SPECTRAL_WINDOW_ID
578 MEpoch midTimeMeas ;
579 Double interval ;
580 getValidTimeRange( midTimeMeas, interval, t2 ) ;
581
582 // add row
583 msSrc.addRow( 1, True ) ;
584
585 // fill SPECTRAL_WINDOW_ID
586 ScalarColumn<Int> sSpwIdCol( msSrc, "SPECTRAL_WINDOW_ID" ) ;
587 sSpwIdCol.put( prevnr, ifno ) ;
588
589 // fill TIME, INTERVAL
590 ScalarMeasColumn<MEpoch> sTimeMeasCol( msSrc, "TIME" ) ;
591 sTimeMeasCol.put( prevnr, midTimeMeas ) ;
592 ScalarColumn<Double> sIntervalCol( msSrc, "INTERVAL" ) ;
593 sIntervalCol.put( prevnr, interval ) ;
594
595 added2++ ;
596 iter2.next() ;
597 }
598
599 // fill NUM_LINES, REST_FREQUENCY, TRANSITION, SYSVEL
600 RefRows rows2( current2, current2+added2-1 ) ;
601 uInt numFreq = restFreq.size() ;
602 Vector<Int> numLines( added2, numFreq ) ;
603 ScalarColumn<Int> numLinesCol( msSrc, "NUM_LINES" ) ;
604 numLinesCol.putColumnCells( rows2, numLines ) ;
605 if ( numFreq != 0 ) {
606 Array<Double> rf( IPosition( 2, numFreq, added2 ) ) ;
607 Array<String> trans( IPosition( 2, numFreq, added2 ) ) ;
608 Array<Double> srcVelArr( IPosition( 2, numFreq, added2 ), srcVel ) ;
609 for ( uInt ifreq = 0 ; ifreq < numFreq ; ifreq++ ) {
610 Slicer slice( Slice(ifreq),Slice(0,added2,1) ) ;
611 rf( slice ) = restFreq[ifreq] ;
612 String transStr = fMolName[ifreq] ;
613 if ( transStr.size() == 0 ) {
614 transStr = molName[ifreq] ;
615 }
616 trans( slice ) = transStr ;
617 }
618 ArrayColumn<Double> sharedDArrCol( msSrc, "REST_FREQUENCY" ) ;
619 sharedDArrCol.putColumnCells( rows2, rf ) ;
620 sharedDArrCol.attach( msSrc, "SYSVEL" ) ;
621 sharedDArrCol.putColumnCells( rows2, srcVelArr ) ;
622 ArrayColumn<String> transCol( msSrc, "TRANSITION" ) ;
623 transCol.putColumnCells( rows2, trans ) ;
624 }
625
626 added1 += added2 ;
627 iter1.next() ;
628 }
629
630 // fill NAME, SOURCE_ID
631 RefRows rows1( current1, current1+added1-1 ) ;
632 Vector<String> nameArr( added1, srcName ) ;
633 Vector<Int> srcIdArr( added1, srcId ) ;
634 ScalarColumn<String> sNameCol( msSrc, "NAME" ) ;
635 ScalarColumn<Int> srcIdCol( msSrc, "SOURCE_ID" ) ;
636 sNameCol.putColumnCells( rows1, nameArr ) ;
637 srcIdCol.putColumnCells( rows1, srcIdArr ) ;
638
639 // fill DIRECTION, PROPER_MOTION
640 Array<Double> srcPMArr( IPosition( 2, 2, added1 ) ) ;
641 Array<Double> srcDirArr( IPosition( 2, 2, added1 ) ) ;
642 for ( uInt i = 0 ; i < 2 ; i++ ) {
643 Slicer slice( Slice(i),Slice(0,added1,1) ) ;
644 srcPMArr( slice ) = srcPM[i] ;
645 srcDirArr( slice ) = srcDir[i] ;
646 }
647 ArrayColumn<Double> sharedDArrCol( msSrc, "DIRECTION" ) ;
648 sharedDArrCol.putColumnCells( rows1, srcDirArr ) ;
649 sharedDArrCol.attach( msSrc, "PROPER_MOTION" ) ;
650 sharedDArrCol.putColumnCells( rows1, srcPMArr ) ;
651
652 // fill TIME, INTERVAL
653
654 // increment srcId if SRCNAME changed
655 srcId++ ;
656
657 added0 += added1 ;
658 iter0.next() ;
659 }
660}
661
662void MSWriter::addFeed( Int id )
663{
664 os_ << "set up FEED subtable" << LogIO::POST ;
665
666 // add row
667 MSFeed msFeed = mstable_->feed() ;
668 msFeed.addRow( 1, True ) ;
669 Int nrow = msFeed.nrow() ;
670
671 MSFeedColumns msFeedCols( mstable_->feed() ) ;
672
673 msFeedCols.feedId().put( nrow-1, id ) ;
674 msFeedCols.antennaId().put( nrow-1, 0 ) ;
675}
676
677void MSWriter::addSpectralWindow( Int spwid, Int freqid )
678{
679 os_ << "set up SPECTRAL_WINDOW subtable" << LogIO::POST ;
680
681 // add row
682 MSSpectralWindow msSpw = mstable_->spectralWindow() ;
683 while( (Int)msSpw.nrow() <= spwid ) {
684 msSpw.addRow( 1, True ) ;
685 }
686
687 MSSpWindowColumns msSpwCols( msSpw ) ;
688
689 STFrequencies stf = table_->frequencies() ;
690
691 // MEAS_FREQ_REF
692 msSpwCols.measFreqRef().put( spwid, stf.getFrame( True ) ) ;
693
694 Double refpix ;
695 Double refval ;
696 Double inc ;
697 stf.getEntry( refpix, refval, inc, (uInt)freqid ) ;
698
699 // NUM_CHAN
700 Int nchan = refpix * 2 ;
701 msSpwCols.numChan().put( spwid, nchan ) ;
702
703 // TOTAL_BANDWIDTH
704 Double bw = nchan * inc ;
705 msSpwCols.totalBandwidth().put( spwid, bw ) ;
706
707 // REF_FREQUENCY
708 Double refFreq = refval - refpix * inc ;
709 msSpwCols.refFrequency().put( spwid, refFreq ) ;
710
711 // NET_SIDEBAND
712 // tentative: USB->0, LSB->1
713 Int netSideband = 0 ;
714 if ( inc < 0 )
715 netSideband = 1 ;
716 msSpwCols.netSideband().put( spwid, netSideband ) ;
717
718 // RESOLUTION, CHAN_WIDTH, EFFECTIVE_BW
719 Vector<Double> sharedDoubleArr( nchan, inc ) ;
720 msSpwCols.resolution().put( spwid, sharedDoubleArr ) ;
721 msSpwCols.chanWidth().put( spwid, sharedDoubleArr ) ;
722 msSpwCols.effectiveBW().put( spwid, sharedDoubleArr ) ;
723
724 // CHAN_FREQ
725 indgen( sharedDoubleArr, refFreq, inc ) ;
726 msSpwCols.chanFreq().put( spwid, sharedDoubleArr ) ;
727}
728
729void MSWriter::addField( Int fid, String fieldname, String srcname, Double t, Vector<Double> rate )
730{
731 os_ << "set up FIELD subtable" << LogIO::POST ;
732
733 MSField msField = mstable_->field() ;
734 while( (Int)msField.nrow() <= fid ) {
735 msField.addRow( 1, True ) ;
736 }
737 MSFieldColumns msFieldCols( msField ) ;
738
739 // Access to SOURCE table
740 MSSource msSrc = mstable_->source() ;
741
742 // fill target row
743 msFieldCols.name().put( fid, fieldname ) ;
744 msFieldCols.time().put( fid, t ) ;
745 Int numPoly = 0 ;
746 if ( anyNE( rate, 0.0 ) )
747 numPoly = 1 ;
748 msFieldCols.numPoly().put( fid, numPoly ) ;
749 MSSourceIndex msSrcIdx( msSrc ) ;
750 Int srcId = -1 ;
751 Vector<Int> srcIdArr = msSrcIdx.matchSourceName( srcname ) ;
752 if ( srcIdArr.size() != 0 ) {
753 srcId = srcIdArr[0] ;
754 MSSource msSrcSel = msSrc( msSrc.col("SOURCE_ID") == srcId ) ;
755 ROMSSourceColumns msSrcCols( msSrcSel ) ;
756 Vector<Double> srcDir = msSrcCols.direction()( 0 ) ;
757 Array<Double> srcDirA( IPosition( 2, 2, 1+numPoly ) ) ;
758 os_ << "srcDirA = " << srcDirA << LogIO::POST ;
759 os_ << "sliced srcDirA = " << srcDirA( Slicer( Slice(0,2,1), Slice(0) ) ) << LogIO::POST ;
760 srcDirA( Slicer( Slice(0,2,1), Slice(0) ) ) = srcDir.reform( IPosition(2,2,1) ) ;
761 os_ << "srcDirA = " << srcDirA << LogIO::POST ;
762 if ( numPoly != 0 )
763 srcDirA( Slicer( Slice(0,2,1), Slice(1) ) ) = rate ;
764 msFieldCols.phaseDir().put( fid, srcDirA ) ;
765 msFieldCols.referenceDir().put( fid, srcDirA ) ;
766 msFieldCols.delayDir().put( fid, srcDirA ) ;
767 }
768 msFieldCols.sourceId().put( fid, srcId ) ;
769}
770
771Int MSWriter::addPolarization( Vector<Int> polnos )
772{
773 os_ << "set up POLARIZATION subtable" << LogIO::POST ;
774
775 os_ << "polnos = " << polnos << LogIO::POST ;
776 MSPolarization msPol = mstable_->polarization() ;
777 uInt nrow = msPol.nrow() ;
778
779 Vector<Int> corrType = toCorrType( polnos ) ;
780 os_ << "corrType = " << corrType << LogIO::POST ;
781
782 MSPolarizationIndex msPolIdx( msPol ) ;
783 Vector<Int> polids = msPolIdx.matchCorrType( corrType ) ;
784 os_ << "polids = " << polids << LogIO::POST ;
785
786 Int polid = -1 ;
787
788 if ( polids.size() == 0 ) {
789 // add row
790 msPol.addRow( 1, True ) ;
791 polid = (Int)nrow ;
792
793 MSPolarizationColumns msPolCols( msPol ) ;
794
795 // CORR_TYPE
796 msPolCols.corrType().put( nrow, corrType ) ;
797
798 // NUM_CORR
799 uInt npol = corrType.size() ;
800 msPolCols.numCorr().put( nrow, npol ) ;
801
802 // CORR_PRODUCT
803 Matrix<Int> corrProd( 2, npol, -1 ) ;
804 if ( npol == 1 ) {
805 corrProd = 0 ;
806 }
807 else if ( npol == 2 ) {
808 corrProd.column( 0 ) = 0 ;
809 corrProd.column( 1 ) = 1 ;
810 }
811 else {
812 corrProd.column( 0 ) = 0 ;
813 corrProd.column( 3 ) = 1 ;
814 corrProd( 0,1 ) = 0 ;
815 corrProd( 1,1 ) = 1 ;
816 corrProd( 0,2 ) = 1 ;
817 corrProd( 1,2 ) = 0 ;
818 }
819 msPolCols.corrProduct().put( nrow, corrProd ) ;
820
821 }
822 else {
823 polid = polids[0] ;
824 }
825
826 return polid ;
827}
828
829Int MSWriter::addDataDescription( Int polid, Int spwid )
830{
831 os_ << "set up SPECTRAL_WINDOW subtable" << LogIO::POST ;
832
833 MSDataDescription msDataDesc = mstable_->dataDescription() ;
834 uInt nrow = msDataDesc.nrow() ;
835
836 MSDataDescIndex msDataDescIdx( msDataDesc ) ;
837
838 Vector<Int> ddids = msDataDescIdx.matchSpwIdAndPolznId( spwid, polid ) ;
839 os_ << "ddids = " << ddids << LogIO::POST ;
840
841 Int ddid = -1 ;
842 if ( ddids.size() == 0 ) {
843 msDataDesc.addRow( 1, True ) ;
844 MSDataDescColumns msDataDescCols( msDataDesc ) ;
845 msDataDescCols.polarizationId().put( nrow, polid ) ;
846 msDataDescCols.spectralWindowId().put( nrow, spwid ) ;
847 ddid = (Int)nrow ;
848 }
849 else {
850 ddid = ddids[0] ;
851 }
852
853 return ddid ;
854}
855
856Vector<Int> MSWriter::toCorrType( Vector<Int> polnos )
857{
858 uInt npol = polnos.size() ;
859 Vector<Int> corrType( npol, Stokes::Undefined ) ;
860
861 if ( npol == 4 ) {
862 if ( polType_ == "linear" ) {
863 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
864 if ( polnos[ipol] == 0 )
865 corrType[ipol] = Stokes::XX ;
866 else if ( polnos[ipol] == 1 )
867 corrType[ipol] = Stokes::XY ;
868 else if ( polnos[ipol] == 2 )
869 corrType[ipol] = Stokes::YX ;
870 else if ( polnos[ipol] == 3 )
871 corrType[ipol] = Stokes::YY ;
872 }
873 }
874 else if ( polType_ == "circular" ) {
875 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
876 if ( polnos[ipol] == 0 )
877 corrType[ipol] = Stokes::RR ;
878 else if ( polnos[ipol] == 1 )
879 corrType[ipol] = Stokes::RL ;
880 else if ( polnos[ipol] == 2 )
881 corrType[ipol] = Stokes::LR ;
882 else if ( polnos[ipol] == 3 )
883 corrType[ipol] = Stokes::LL ;
884 }
885 }
886 else if ( polType_ == "stokes" ) {
887 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
888 if ( polnos[ipol] == 0 )
889 corrType[ipol] = Stokes::I ;
890 else if ( polnos[ipol] == 1 )
891 corrType[ipol] = Stokes::Q ;
892 else if ( polnos[ipol] == 2 )
893 corrType[ipol] = Stokes::U ;
894 else if ( polnos[ipol] == 3 )
895 corrType[ipol] = Stokes::V ;
896 }
897 }
898 }
899 else if ( npol == 2 ) {
900 if ( polType_ == "linear" ) {
901 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
902 if ( polnos[ipol] == 0 )
903 corrType[ipol] = Stokes::XX ;
904 else if ( polnos[ipol] == 1 )
905 corrType[ipol] = Stokes::YY ;
906 }
907 }
908 else if ( polType_ == "circular" ) {
909 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
910 if ( polnos[ipol] == 0 )
911 corrType[ipol] = Stokes::RR ;
912 else if ( polnos[ipol] == 1 )
913 corrType[ipol] = Stokes::LL ;
914 }
915 }
916 else if ( polType_ == "stokes" ) {
917 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
918 if ( polnos[ipol] == 0 )
919 corrType[ipol] = Stokes::I ;
920 else if ( polnos[ipol] == 1 )
921 corrType[ipol] = Stokes::V ;
922 }
923 }
924 else if ( polType_ == "linpol" ) {
925 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
926 if ( polnos[ipol] == 1 )
927 corrType[ipol] = Stokes::Plinear ;
928 else if ( polnos[ipol] == 2 )
929 corrType[ipol] = Stokes::Pangle ;
930 }
931 }
932 }
933 else if ( npol == 1 ) {
934 if ( polType_ == "linear" )
935 corrType[0] = Stokes::XX ;
936 else if ( polType_ == "circular" )
937 corrType[0] = Stokes::RR ;
938 else if ( polType_ == "stokes" )
939 corrType[0] = Stokes::I ;
940 }
941
942 return corrType ;
943}
944
945void MSWriter::getValidTimeRange( MEpoch &me, Double &interval, Table &tab )
946{
947 ROScalarMeasColumn<MEpoch> timeMeasCol( tab, "TIME" ) ;
948 ROScalarColumn<Double> timeCol( tab, "TIME" ) ;
949 String refStr = timeMeasCol( 0 ).getRefString() ;
950 Vector<Double> timeArr = timeCol.getColumn() ;
951 MEpoch::Types meType ;
952 MEpoch::getType( meType, refStr ) ;
953 Unit tUnit = timeMeasCol.measDesc().getUnits()( 0 ) ;
954 Double minTime ;
955 Double maxTime ;
956 minMax( minTime, maxTime, timeArr ) ;
957 Double midTime = 0.5 * ( minTime + maxTime ) ;
958 me = MEpoch( Quantity( midTime, tUnit ), meType ) ;
959 interval = Quantity( maxTime-minTime, tUnit ).getValue( "s" ) ;
960}
961
962}
Note: See TracBrowser for help on using the repository browser.