source: trunk/src/MSWriter.cpp@ 1974

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

New Development: Yes

JIRA Issue: Yes CAS-2718

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: New class msfiller and mswriter added

Test Programs: List test programs

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

New filler/writer specific for Scantable-MS conversion defined.
This is not called from scantable constructor right now.
However, you can call it by explicitly importing msfiller/mswriter.
File size: 23.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/TableDesc.h>
21#include <tables/Tables/SetupNewTab.h>
22#include <tables/Tables/TableIter.h>
23#include <tables/Tables/RefRows.h>
24
25#include <ms/MeasurementSets/MeasurementSet.h>
26#include <ms/MeasurementSets/MSColumns.h>
27#include <ms/MeasurementSets/MSPolIndex.h>
28#include <ms/MeasurementSets/MSDataDescIndex.h>
29
30#include "MSWriter.h"
31#include "STHeader.h"
32#include "STFrequencies.h"
33
34using namespace casa ;
35
36namespace asap
37{
38
39MSWriter::MSWriter(CountedPtr<Scantable> stable)
40 : table_(stable)
41{
42 os_ = LogIO() ;
43 os_.origin( LogOrigin( "MSWriter", "MSWriter()", WHERE ) ) ;
44 os_ << "MSWriter::MSWriter()" << LogIO::POST ;
45
46 // initialize writer
47 init() ;
48}
49
50MSWriter::~MSWriter()
51{
52 os_.origin( LogOrigin( "MSWriter", "~MSWriter()", WHERE ) ) ;
53 os_ << "MSWriter::~MSWriter()" << LogIO::POST ;
54}
55
56bool MSWriter::write(const string& filename, const Record& rec)
57{
58 os_.origin( LogOrigin( "MSWriter", "write()", WHERE ) ) ;
59 os_ << "MSWriter::write()" << LogIO::POST ;
60
61 filename_ = filename ;
62
63 // parsing MS options
64 Bool overwrite = False ;
65 if ( rec.isDefined( "ms" ) ) {
66 Record msrec = rec.asRecord( "ms" ) ;
67 if ( msrec.isDefined( "overwrite" ) ) {
68 overwrite = msrec.asBool( "overwrite" ) ;
69 }
70 }
71
72 os_ << "Parsing MS options" << endl ;
73 os_ << " overwrite = " << overwrite << LogIO::POST ;
74
75 File file( filename_ ) ;
76 if ( file.exists() ) {
77 if ( overwrite ) {
78 os_ << filename_ << " exists. Overwrite existing data... " << LogIO::POST ;
79 if ( file.isRegular() ) RegularFile(file).remove() ;
80 else if ( file.isDirectory() ) Directory(file).removeRecursive() ;
81 else SymLink(file).remove() ;
82 }
83 else {
84 os_ << LogIO::SEVERE << "ERROR: " << filename_ << " exists..." << LogIO::POST ;
85 return False ;
86 }
87 }
88
89 // set up MS
90 setupMS() ;
91
92 // subtables
93 // OBSERVATION
94 fillObservation() ;
95
96 // ANTENNA
97 fillAntenna() ;
98
99 // MAIN
100 // Iterate over several ids
101 Vector<uInt> processedFreqId( 0 ) ;
102 Int defaultFieldId = 0 ;
103 //
104 // ITERATION: FIELDNAME
105 //
106 Int added0 = 0 ;
107 Int current0 = mstable_->nrow() ;
108 TableIterator iter0( table_->table(), "FIELDNAME" ) ;
109 while( !iter0.pastEnd() ) {
110 Table t0( iter0.table() ) ;
111 ROScalarColumn<String> fieldNameCol( t0, "FIELDNAME" ) ;
112 String fieldName = fieldNameCol(0) ;
113 String::size_type pos = fieldName.find( "__" ) ;
114 Int fieldId = -1 ;
115 if ( pos != String::npos ) {
116 os_ << "fieldName.substr( pos+2 )=" << fieldName.substr( pos+2 ) << LogIO::POST ;
117 fieldId = String::toInt( fieldName.substr( pos+2 ) ) ;
118 fieldName = fieldName.substr( 0, pos ) ;
119 }
120 else {
121 os_ << "use default field id" << LogIO::POST ;
122 fieldId = defaultFieldId ;
123 defaultFieldId++ ;
124 }
125 os_ << "fieldId" << fieldId << ": " << fieldName << LogIO::POST ;
126 //
127 // ITERATION: BEAMNO
128 //
129 Int added1 = 0 ;
130 Int current1 = mstable_->nrow() ;
131 TableIterator iter1( t0, "BEAMNO" ) ;
132 while( !iter1.pastEnd() ) {
133 Table t1( iter1.table() ) ;
134 ROScalarColumn<uInt> beamNoCol( t1, "BEAMNO" ) ;
135 uInt beamNo = beamNoCol(0) ;
136 os_ << "beamNo = " << beamNo << LogIO::POST ;
137 //
138 // ITERATION: SCANNO
139 //
140 Int added2 = 0 ;
141 Int current2 = mstable_->nrow() ;
142 TableIterator iter2( t1, "SCANNO" ) ;
143 while( !iter2.pastEnd() ) {
144 Table t2( iter2.table() ) ;
145 ROScalarColumn<uInt> scanNoCol( t2, "SCANNO" ) ;
146 uInt scanNo = scanNoCol(0) ;
147 os_ << "scanNo = " << scanNo << LogIO::POST ;
148 //
149 // ITERATION: CYCLENO
150 //
151 Int added3 = 0 ;
152 Int current3 = mstable_->nrow() ;
153 TableIterator iter3( t2, "CYCLENO" ) ;
154 while( !iter3.pastEnd() ) {
155 Table t3( iter3.table() ) ;
156 //
157 // ITERATION: IFNO
158 //
159 Int added4 = 0 ;
160 Int current4 = mstable_->nrow() ;
161 TableIterator iter4( t3, "IFNO" ) ;
162 while( !iter4.pastEnd() ) {
163 Table t4( iter4.table() ) ;
164 ROScalarColumn<uInt> ifNoCol( t4, "IFNO" ) ;
165 uInt ifNo = ifNoCol(0) ;
166 os_ << "ifNo = " << ifNo << LogIO::POST ;
167 ROScalarColumn<uInt> freqIdCol( t4, "FREQ_ID" ) ;
168 uInt freqId = freqIdCol(0) ;
169 os_ << "freqId = " << freqId << LogIO::POST ;
170 //
171 // ITERATION: TIME
172 //
173 Int added5 = 0 ;
174 Int current5 = mstable_->nrow() ;
175 TableIterator iter5( t4, "TIME" ) ;
176 while( !iter5.pastEnd() ) {
177 Table t5( iter5.table().sort("POLNO") ) ;
178 Int prevnr = mstable_->nrow() ;
179 Int nrow = t5.nrow() ;
180 os_ << "nrow = " << nrow << LogIO::POST ;
181
182 // add row
183 mstable_->addRow( 1, True ) ;
184
185 Vector<Int> polnos( nrow ) ;
186 indgen( polnos, 0 ) ;
187 Int polid = addPolarization( polnos ) ;
188 os_ << "polid = " << polid << LogIO::POST ;
189 //
190 // LOOP: POLNO
191 //
192 for ( Int ipol = 0 ; ipol < nrow ; ipol++ ) {
193 }
194
195 // TIME and TIME_CENTROID
196 ROScalarMeasColumn<MEpoch> timeCol( t5, "TIME" ) ;
197 ScalarMeasColumn<MEpoch> msTimeCol( *mstable_, "TIME" ) ;
198 msTimeCol.put( prevnr, timeCol(0) ) ;
199 msTimeCol.attach( *mstable_, "TIME_CENTROID" ) ;
200 msTimeCol.put( prevnr, timeCol(0) ) ;
201
202 // INTERVAL and EXPOSURE
203 ROScalarColumn<Double> intervalCol( t5, "INTERVAL" ) ;
204 ScalarColumn<Double> msIntervalCol( *mstable_, "INTERVAL" ) ;
205 msIntervalCol.put( prevnr, intervalCol(0) ) ;
206 msIntervalCol.attach( *mstable_, "EXPOSURE" ) ;
207 msIntervalCol.put( prevnr, intervalCol(0) ) ;
208
209 // add DATA_DESCRIPTION row
210 Int ddid = addDataDescription( polid, ifNo ) ;
211 os_ << "ddid = " << ddid << LogIO::POST ;
212 ScalarColumn<Int> ddIdCol( *mstable_, "DATA_DESC_ID" ) ;
213 ddIdCol.put( prevnr, ddid ) ;
214
215 added5 += 1 ;
216 os_ << "added5 = " << added5 << " current5 = " << current5 << LogIO::POST ;
217 iter5.next() ;
218 }
219
220 // add SPECTRAL_WINDOW row
221 if ( allNE( processedFreqId, freqId ) ) {
222 uInt vsize = processedFreqId.size() ;
223 processedFreqId.resize( vsize+1, True ) ;
224 processedFreqId[vsize] = freqId ;
225 addSpectralWindow( ifNo, freqId ) ;
226 }
227
228 added4 += added5 ;
229 os_ << "added4 = " << added4 << " current4 = " << current4 << LogIO::POST ;
230 iter4.next() ;
231 }
232 added3 += added4 ;
233 os_ << "added3 = " << added3 << " current3 = " << current3 << LogIO::POST ;
234 iter3.next() ;
235 }
236
237 // SCAN_NUMBER
238 RefRows rows3( current3, current3+added3-1 ) ;
239 Vector<Int> scanNum( added3, scanNo ) ;
240 ScalarColumn<Int> scanNumCol( *mstable_, "SCAN_NUMBER" ) ;
241 scanNumCol.putColumnCells( rows3, scanNum ) ;
242
243 added2 += added3 ;
244 os_ << "added2 = " << added2 << " current2 = " << current2 << LogIO::POST ;
245 iter2.next() ;
246 }
247
248 // FEED1 and FEED2
249 RefRows rows2( current2, current2+added2-1 ) ;
250 Vector<Int> feedId( added2, beamNo ) ;
251 ScalarColumn<Int> feedCol( *mstable_, "FEED1" ) ;
252 feedCol.putColumnCells( rows2, feedId ) ;
253 feedCol.attach( *mstable_, "FEED2" ) ;
254 feedCol.putColumnCells( rows2, feedId ) ;
255
256 // add FEED row
257 addFeed( beamNo ) ;
258
259 added1 += added2 ;
260 os_ << "added1 = " << added1 << " current1 = " << current1 << LogIO::POST ;
261 iter1.next() ;
262 }
263
264 // FIELD_ID
265 RefRows rows1( current1, current1+added1-1 ) ;
266 Vector<Int> fieldIds( added1, fieldId ) ;
267 ScalarColumn<Int> fieldIdCol( *mstable_, "FIELD_ID" ) ;
268 fieldIdCol.putColumnCells( rows1, fieldIds ) ;
269
270 added0 += added1 ;
271 os_ << "added0 = " << added0 << " current0 = " << current0 << LogIO::POST ;
272 iter0.next() ;
273 }
274
275 // OBSERVATION_ID is always 0
276 ScalarColumn<Int> sharedIntCol( *mstable_, "OBSERVATION_ID" ) ;
277 Vector<Int> sharedIntArr( added0, 0 ) ;
278 sharedIntCol.putColumn( sharedIntArr ) ;
279
280 // ANTENNA1 and ANTENNA2 are always 0
281 sharedIntArr = 0 ;
282 sharedIntCol.attach( *mstable_, "ANTENNA1" ) ;
283 sharedIntCol.putColumn( sharedIntArr ) ;
284 sharedIntCol.attach( *mstable_, "ANTENNA2" ) ;
285 sharedIntCol.putColumn( sharedIntArr ) ;
286
287 // ARRAY_ID is tentatively set to 0
288 sharedIntArr = 0 ;
289 sharedIntCol.attach( *mstable_, "ARRAY_ID" ) ;
290 sharedIntCol.putColumn( sharedIntArr ) ;
291
292 return True ;
293}
294
295void MSWriter::init()
296{
297// os_.origin( LogOrigin( "MSWriter", "init()", WHERE ) ) ;
298 os_ << "MSWriter::init()" << LogIO::POST ;
299
300 // access to scantable
301 header_ = table_->getHeader() ;
302
303 // FLOAT_DATA? or DATA?
304 if ( header_.npol > 2 ) {
305 isFloatData_ = False ;
306 isData_ = True ;
307 }
308 else {
309 isFloatData_ = True ;
310 isData_ = False ;
311 }
312
313 // polarization type
314 polType_ = header_.poltype ;
315 if ( polType_ == "" )
316 polType_ = "stokes" ;
317 else if ( polType_.find( "linear" ) != String::npos )
318 polType_ = "linear" ;
319 else if ( polType_.find( "circular" ) != String::npos )
320 polType_ = "circular" ;
321 else if ( polType_.find( "stokes" ) != String::npos )
322 polType_ = "stokes" ;
323 else if ( polType_.find( "linpol" ) != String::npos )
324 polType_ = "linpol" ;
325 else
326 polType_ = "notype" ;
327
328}
329
330void MSWriter::setupMS()
331{
332// os_.origin( LogOrigin( "MSWriter", "setupMS()", WHERE ) ) ;
333 os_ << "MSWriter::setupMS()" << LogIO::POST ;
334
335 TableDesc msDesc = MeasurementSet::requiredTableDesc() ;
336
337 // any additional columns?
338 // TODO: add FLOAT_DATA column if npol == 2 otherwise add DATA column
339
340 SetupNewTable newtab( filename_, msDesc, Table::New ) ;
341
342 mstable_ = new MeasurementSet( newtab ) ;
343
344 // create subtables
345 TableDesc antennaDesc = MSAntenna::requiredTableDesc() ;
346 SetupNewTable antennaTab( mstable_->antennaTableName(), antennaDesc, Table::New ) ;
347 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::ANTENNA ), Table( antennaTab ) ) ;
348
349 TableDesc dataDescDesc = MSDataDescription::requiredTableDesc() ;
350 SetupNewTable dataDescTab( mstable_->dataDescriptionTableName(), dataDescDesc, Table::New ) ;
351 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::DATA_DESCRIPTION ), Table( dataDescTab ) ) ;
352
353 TableDesc dopplerDesc = MSDoppler::requiredTableDesc() ;
354 SetupNewTable dopplerTab( mstable_->dopplerTableName(), dopplerDesc, Table::New ) ;
355 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::DOPPLER ), Table( dopplerTab ) ) ;
356
357 TableDesc feedDesc = MSFeed::requiredTableDesc() ;
358 SetupNewTable feedTab( mstable_->feedTableName(), feedDesc, Table::New ) ;
359 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FEED ), Table( feedTab ) ) ;
360
361 TableDesc fieldDesc = MSField::requiredTableDesc() ;
362 SetupNewTable fieldTab( mstable_->fieldTableName(), fieldDesc, Table::New ) ;
363 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FIELD ), Table( fieldTab ) ) ;
364
365 TableDesc flagCmdDesc = MSFlagCmd::requiredTableDesc() ;
366 SetupNewTable flagCmdTab( mstable_->flagCmdTableName(), flagCmdDesc, Table::New ) ;
367 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FLAG_CMD ), Table( flagCmdTab ) ) ;
368
369 TableDesc freqOffsetDesc = MSFreqOffset::requiredTableDesc() ;
370 SetupNewTable freqOffsetTab( mstable_->freqOffsetTableName(), freqOffsetDesc, Table::New ) ;
371 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FREQ_OFFSET ), Table( freqOffsetTab ) ) ;
372
373 TableDesc historyDesc = MSHistory::requiredTableDesc() ;
374 SetupNewTable historyTab( mstable_->historyTableName(), historyDesc, Table::New ) ;
375 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::HISTORY ), Table( historyTab ) ) ;
376
377 TableDesc observationDesc = MSObservation::requiredTableDesc() ;
378 SetupNewTable observationTab( mstable_->observationTableName(), observationDesc, Table::New ) ;
379 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::OBSERVATION ), Table( observationTab ) ) ;
380
381 TableDesc pointingDesc = MSPointing::requiredTableDesc() ;
382 SetupNewTable pointingTab( mstable_->pointingTableName(), pointingDesc, Table::New ) ;
383 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::POINTING ), Table( pointingTab ) ) ;
384
385 TableDesc polarizationDesc = MSPolarization::requiredTableDesc() ;
386 SetupNewTable polarizationTab( mstable_->polarizationTableName(), polarizationDesc, Table::New ) ;
387 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::POLARIZATION ), Table( polarizationTab ) ) ;
388
389 TableDesc processorDesc = MSProcessor::requiredTableDesc() ;
390 SetupNewTable processorTab( mstable_->processorTableName(), processorDesc, Table::New ) ;
391 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::PROCESSOR ), Table( processorTab ) ) ;
392
393 TableDesc sourceDesc = MSSource::requiredTableDesc() ;
394 SetupNewTable sourceTab( mstable_->sourceTableName(), sourceDesc, Table::New ) ;
395 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SOURCE ), Table( sourceTab ) ) ;
396
397 TableDesc spwDesc = MSSpectralWindow::requiredTableDesc() ;
398 SetupNewTable spwTab( mstable_->spectralWindowTableName(), spwDesc, Table::New ) ;
399 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SPECTRAL_WINDOW ), Table( spwTab ) ) ;
400
401 TableDesc stateDesc = MSState::requiredTableDesc() ;
402 SetupNewTable stateTab( mstable_->stateTableName(), stateDesc, Table::New ) ;
403 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::STATE ), Table( stateTab ) ) ;
404
405 // TODO: add TCAL_SPECTRUM and TSYS_SPECTRUM if necessary
406 TableDesc sysCalDesc = MSSysCal::requiredTableDesc() ;
407 SetupNewTable sysCalTab( mstable_->sysCalTableName(), sysCalDesc, Table::New ) ;
408 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SYSCAL ), Table( sysCalTab ) ) ;
409
410 TableDesc weatherDesc = MSWeather::requiredTableDesc() ;
411 SetupNewTable weatherTab( mstable_->weatherTableName(), weatherDesc, Table::New ) ;
412 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::WEATHER ), Table( weatherTab ) ) ;
413
414 mstable_->initRefs() ;
415
416}
417
418void MSWriter::fillObservation()
419{
420 os_ << "set up OBSERVATION subtable" << LogIO::POST ;
421 // only 1 row
422 mstable_->observation().addRow( 1, True ) ;
423 MSObservationColumns msObsCols( mstable_->observation() ) ;
424 msObsCols.observer().put( 0, header_.observer ) ;
425 // tentatively put antennaname (from ANTENNA subtable)
426 String hAntennaName = header_.antennaname ;
427 String::size_type pos = hAntennaName.find( "//" ) ;
428 String telescopeName ;
429 if ( pos != String::npos ) {
430 telescopeName = hAntennaName.substr( 0, pos ) ;
431 }
432 else {
433 pos = hAntennaName.find( "@" ) ;
434 telescopeName = hAntennaName.substr( 0, pos ) ;
435 }
436 os_ << "telescopeName = " << telescopeName << LogIO::POST ;
437 msObsCols.telescopeName().put( 0, telescopeName ) ;
438 msObsCols.project().put( 0, header_.project ) ;
439 //ScalarMeasColumn<MEpoch> timeCol( table_->table().sort("TIME"), "TIME" ) ;
440 Table sortedtable = table_->table().sort("TIME") ;
441 ScalarMeasColumn<MEpoch> timeCol( sortedtable, "TIME" ) ;
442 Vector<MEpoch> trange( 2 ) ;
443 trange[0] = timeCol( 0 ) ;
444 trange[1] = timeCol( table_->nrow()-1 ) ;
445 msObsCols.timeRangeMeas().put( 0, trange ) ;
446}
447
448void MSWriter::fillAntenna()
449{
450 os_ << "set up ANTENNA subtable" << LogIO::POST ;
451 // only 1 row
452 mstable_->antenna().addRow( 1, True ) ;
453 MSAntennaColumns msAntCols( mstable_->antenna() ) ;
454
455 String hAntennaName = header_.antennaname ;
456 String::size_type pos = hAntennaName.find( "//" ) ;
457 String antennaName ;
458 String stationName ;
459 if ( pos != String::npos ) {
460 hAntennaName = hAntennaName.substr( pos+2 ) ;
461 }
462 pos = hAntennaName.find( "@" ) ;
463 if ( pos != String::npos ) {
464 antennaName = hAntennaName.substr( 0, pos ) ;
465 stationName = hAntennaName.substr( pos+1 ) ;
466 }
467 else {
468 antennaName = hAntennaName ;
469 stationName = hAntennaName ;
470 }
471 os_ << "antennaName = " << antennaName << LogIO::POST ;
472 os_ << "stationName = " << stationName << LogIO::POST ;
473
474 msAntCols.name().put( 0, antennaName ) ;
475 msAntCols.station().put( 0, stationName ) ;
476
477 os_ << "antennaPosition = " << header_.antennaposition << LogIO::POST ;
478
479 msAntCols.position().put( 0, header_.antennaposition ) ;
480}
481
482void MSWriter::addFeed( Int id )
483{
484 os_ << "set up FEED subtable" << LogIO::POST ;
485
486 // add row
487 MSFeed msFeed = mstable_->feed() ;
488 msFeed.addRow( 1, True ) ;
489 Int nrow = msFeed.nrow() ;
490
491 MSFeedColumns msFeedCols( mstable_->feed() ) ;
492
493 msFeedCols.feedId().put( nrow-1, id ) ;
494 msFeedCols.antennaId().put( nrow-1, 0 ) ;
495}
496
497void MSWriter::addSpectralWindow( Int spwid, Int freqid )
498{
499 os_ << "set up SPECTRAL_WINDOW subtable" << LogIO::POST ;
500
501 // add row
502 MSSpectralWindow msSpw = mstable_->spectralWindow() ;
503 while( (Int)msSpw.nrow() <= spwid ) {
504 msSpw.addRow( 1, True ) ;
505 }
506
507 MSSpWindowColumns msSpwCols( msSpw ) ;
508
509 STFrequencies stf = table_->frequencies() ;
510
511 // MEAS_FREQ_REF
512 msSpwCols.measFreqRef().put( spwid, stf.getFrame( True ) ) ;
513
514 Double refpix ;
515 Double refval ;
516 Double inc ;
517 stf.getEntry( refpix, refval, inc, (uInt)freqid ) ;
518
519 // NUM_CHAN
520 Int nchan = refpix * 2 ;
521 msSpwCols.numChan().put( spwid, nchan ) ;
522
523 // TOTAL_BANDWIDTH
524 Double bw = nchan * inc ;
525 msSpwCols.totalBandwidth().put( spwid, bw ) ;
526
527 // REF_FREQUENCY
528 Double refFreq = refval - refpix * inc ;
529 msSpwCols.refFrequency().put( spwid, refFreq ) ;
530
531 // NET_SIDEBAND
532 // tentative: USB->0, LSB->1
533 Int netSideband = 0 ;
534 if ( inc < 0 )
535 netSideband = 1 ;
536 msSpwCols.netSideband().put( spwid, netSideband ) ;
537
538 // RESOLUTION, CHAN_WIDTH, EFFECTIVE_BW
539 Vector<Double> sharedDoubleArr( nchan, inc ) ;
540 msSpwCols.resolution().put( spwid, sharedDoubleArr ) ;
541 msSpwCols.chanWidth().put( spwid, sharedDoubleArr ) ;
542 msSpwCols.effectiveBW().put( spwid, sharedDoubleArr ) ;
543
544 // CHAN_FREQ
545 indgen( sharedDoubleArr, refFreq, inc ) ;
546 msSpwCols.chanFreq().put( spwid, sharedDoubleArr ) ;
547}
548
549Int MSWriter::addPolarization( Vector<Int> polnos )
550{
551 os_ << "set up SPECTRAL_WINDOW subtable" << LogIO::POST ;
552
553 os_ << "polnos = " << polnos << LogIO::POST ;
554 MSPolarization msPol = mstable_->polarization() ;
555 uInt nrow = msPol.nrow() ;
556
557 Vector<Int> corrType = toCorrType( polnos ) ;
558 os_ << "corrType = " << corrType << LogIO::POST ;
559
560 MSPolarizationIndex msPolIdx( msPol ) ;
561 Vector<Int> polids = msPolIdx.matchCorrType( corrType ) ;
562 os_ << "polids = " << polids << LogIO::POST ;
563
564 Int polid = -1 ;
565
566 if ( polids.size() == 0 ) {
567 // add row
568 msPol.addRow( 1, True ) ;
569 polid = (Int)nrow ;
570
571 MSPolarizationColumns msPolCols( msPol ) ;
572
573 // CORR_TYPE
574 msPolCols.corrType().put( nrow, corrType ) ;
575
576 // NUM_CORR
577 uInt npol = corrType.size() ;
578 msPolCols.numCorr().put( nrow, npol ) ;
579
580 // CORR_PRODUCT
581 Matrix<Int> corrProd( 2, npol, -1 ) ;
582 if ( npol == 1 ) {
583 corrProd = 0 ;
584 }
585 else if ( npol == 2 ) {
586 corrProd.column(0) = 0 ;
587 corrProd.column(1) = 1 ;
588 }
589 else {
590 corrProd.column(0) = 0 ;
591 corrProd.column(3) = 1 ;
592 corrProd(0,1) = 0 ;
593 corrProd(1,1) = 1 ;
594 corrProd(0,2) = 1 ;
595 corrProd(1,2) = 0 ;
596 }
597 msPolCols.corrProduct().put( nrow, corrProd ) ;
598
599 }
600 else {
601 polid = polids[0] ;
602 }
603
604 return polid ;
605}
606
607Vector<Int> MSWriter::toCorrType( Vector<Int> polnos )
608{
609 uInt npol = polnos.size() ;
610 Vector<Int> corrType( npol, Stokes::Undefined ) ;
611
612 if ( npol == 4 ) {
613 if ( polType_ == "linear" ) {
614 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
615 if ( polnos[ipol] == 0 )
616 corrType[ipol] = Stokes::XX ;
617 else if ( polnos[ipol] == 1 )
618 corrType[ipol] = Stokes::XY ;
619 else if ( polnos[ipol] == 2 )
620 corrType[ipol] = Stokes::YX ;
621 else if ( polnos[ipol] == 3 )
622 corrType[ipol] = Stokes::YY ;
623 }
624 }
625 else if ( polType_ == "circular" ) {
626 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
627 if ( polnos[ipol] == 0 )
628 corrType[ipol] = Stokes::RR ;
629 else if ( polnos[ipol] == 1 )
630 corrType[ipol] = Stokes::RL ;
631 else if ( polnos[ipol] == 2 )
632 corrType[ipol] = Stokes::LR ;
633 else if ( polnos[ipol] == 3 )
634 corrType[ipol] = Stokes::LL ;
635 }
636 }
637 else if ( polType_ == "stokes" ) {
638 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
639 if ( polnos[ipol] == 0 )
640 corrType[ipol] = Stokes::I ;
641 else if ( polnos[ipol] == 1 )
642 corrType[ipol] = Stokes::Q ;
643 else if ( polnos[ipol] == 2 )
644 corrType[ipol] = Stokes::U ;
645 else if ( polnos[ipol] == 3 )
646 corrType[ipol] = Stokes::V ;
647 }
648 }
649 }
650 else if ( npol == 2 ) {
651 if ( polType_ == "linear" ) {
652 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
653 if ( polnos[ipol] == 0 )
654 corrType[ipol] = Stokes::XX ;
655 else if ( polnos[ipol] == 1 )
656 corrType[ipol] = Stokes::YY ;
657 }
658 }
659 else if ( polType_ == "circular" ) {
660 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
661 if ( polnos[ipol] == 0 )
662 corrType[ipol] = Stokes::RR ;
663 else if ( polnos[ipol] == 1 )
664 corrType[ipol] = Stokes::LL ;
665 }
666 }
667 else if ( polType_ == "stokes" ) {
668 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
669 if ( polnos[ipol] == 0 )
670 corrType[ipol] = Stokes::I ;
671 else if ( polnos[ipol] == 1 )
672 corrType[ipol] = Stokes::V ;
673 }
674 }
675 else if ( polType_ == "linpol" ) {
676 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
677 if ( polnos[ipol] == 1 )
678 corrType[ipol] = Stokes::Plinear ;
679 else if ( polnos[ipol] == 2 )
680 corrType[ipol] = Stokes::Pangle ;
681 }
682 }
683 }
684 else if ( npol == 1 ) {
685 if ( polType_ == "linear" )
686 corrType[0] = Stokes::XX ;
687 else if ( polType_ == "circular" )
688 corrType[0] = Stokes::RR ;
689 else if ( polType_ == "stokes" )
690 corrType[0] = Stokes::I ;
691 }
692
693 return corrType ;
694}
695
696Int MSWriter::addDataDescription( Int polid, Int spwid )
697{
698 os_ << "set up SPECTRAL_WINDOW subtable" << LogIO::POST ;
699
700 MSDataDescription msDataDesc = mstable_->dataDescription() ;
701 uInt nrow = msDataDesc.nrow() ;
702
703 MSDataDescIndex msDataDescIdx( msDataDesc ) ;
704
705 Vector<Int> ddids = msDataDescIdx.matchSpwIdAndPolznId( spwid, polid ) ;
706 os_ << "ddids = " << ddids << LogIO::POST ;
707
708 Int ddid = -1 ;
709 if ( ddids.size() == 0 ) {
710 msDataDesc.addRow( 1, True ) ;
711 MSDataDescColumns msDataDescCols( msDataDesc ) ;
712 msDataDescCols.polarizationId().put( nrow, polid ) ;
713 msDataDescCols.spectralWindowId().put( nrow, spwid ) ;
714 ddid = (Int)nrow ;
715 }
716 else {
717 ddid = ddids[0] ;
718 }
719
720 return ddid ;
721}
722}
Note: See TracBrowser for help on using the repository browser.