source: branches/casa-prerelease/pre-asap/src/MSWriter.cpp@ 2325

Last change on this file since 2325 was 2319, checked in by Takeshi Nakazato, 13 years ago

New Development: No

JIRA Issue: No

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

Merge bug fix on trunk (r2318).


File size: 73.3 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#include <assert.h>
14
15#include <casa/OS/File.h>
16#include <casa/OS/RegularFile.h>
17#include <casa/OS/Directory.h>
18#include <casa/OS/SymLink.h>
19#include <casa/BasicSL/String.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#include "MathUtils.h"
41#include "TableTraverse.h"
42
43using namespace casa ;
44using namespace std ;
45
46namespace asap {
47
48class PolarizedComponentHolder {
49public:
50 PolarizedComponentHolder()
51 : nchan(0),
52 maxnpol(4)
53 {
54 reset() ;
55 }
56 PolarizedComponentHolder( uInt n )
57 : nchan(n),
58 maxnpol(4)
59 {
60 reset() ;
61 }
62
63 void reset()
64 {
65 npol = 0 ;
66 data.clear() ;
67 flag.clear() ;
68 flagrow = False ;
69 polnos.resize() ;
70 }
71
72 void accumulate( uInt id, Vector<Float> sp, Vector<Bool> fl, Bool flr )
73 {
74 map< uInt,Vector<Float> >::iterator itr = data.find( id ) ;
75 if ( id < maxnpol && itr == data.end() ) {
76 addPol( id ) ;
77 accumulateData( id, sp ) ;
78 accumulateFlag( id, fl ) ;
79 accumulateFlagRow( flr ) ;
80 npol++ ;
81 }
82 }
83
84 void setNchan( uInt n ) { nchan = n ; }
85 uInt nPol() { return npol ; }
86 uInt nChan() { return nchan ; }
87 Vector<uInt> polNos() { return polnos ; }
88 Vector<Float> getWeight() { return Vector<Float>( npol, 1.0 ) ; }
89 Vector<Float> getSigma() { return Vector<Float>( npol, 1.0 ) ; }
90 Bool getFlagRow() { return flagrow ; }
91 Cube<Bool> getFlagCategory() { return Cube<Bool>( npol, nchan, 1, False ) ; }
92 Matrix<Float> getData()
93 {
94 Matrix<Float> v( npol, nchan ) ;
95 for ( map< uInt,Vector<Float> >::iterator i = data.begin() ; i != data.end() ; i++ ) {
96 v.row( i->first ) = i->second ;
97 }
98 return v ;
99 }
100 Matrix<Bool> getFlag()
101 {
102 Matrix<Bool> v( npol, nchan ) ;
103 for ( map< uInt,Vector<Bool> >::iterator i = flag.begin() ; i != flag.end() ; i++ ) {
104 v.row( i->first ) = i->second ;
105 }
106 return v ;
107 }
108 Matrix<Complex> getComplexData()
109 {
110 Matrix<Complex> v( npol, nchan ) ;
111 Matrix<Float> dummy( 2, nchan, 0.0 ) ;
112 map< uInt,Vector<Float> >::iterator itr0 = data.find( 0 ) ;
113 map< uInt,Vector<Float> >::iterator itr1 = data.find( 1 ) ;
114 if ( itr0 != data.end() ) {
115 dummy.row( 0 ) = itr0->second ;
116 v.row( 0 ) = RealToComplex( dummy ) ;
117 }
118 if ( itr1 != data.end() ) {
119 dummy.row( 0 ) = itr1->second ;
120 v.row( npol-1 ) = RealToComplex( dummy ) ;
121 }
122 itr0 = data.find( 2 ) ;
123 itr1 = data.find( 3 ) ;
124 if ( itr0 != data.end() && itr1 != data.end() ) {
125 dummy.row( 0 ) = itr0->second ;
126 dummy.row( 1 ) = itr1->second ;
127 v.row( 1 ) = RealToComplex( dummy ) ;
128 v.row( 2 ) = conj( v.row( 1 ) ) ;
129 }
130 return v ;
131 }
132
133 Matrix<Bool> getComplexFlag()
134 {
135 Matrix<Bool> tmp = getFlag() ;
136 Matrix<Bool> v( npol, nchan ) ;
137 v.row( 0 ) = tmp.row( 0 ) ;
138 if ( npol == 2 ) {
139 v.row( npol-1 ) = tmp.row( 1 ) ;
140 }
141 else if ( npol > 2 ) {
142 v.row( npol-1 ) = tmp.row( 1 ) ;
143 v.row( 1 ) = tmp.row( 2 ) || tmp.row( 3 ) ;
144 v.row( 2 ) = v.row( 1 ) ;
145 }
146 return v ;
147 }
148
149private:
150 void accumulateData( uInt &id, Vector<Float> &v )
151 {
152 data.insert( pair< uInt,Vector<Float> >( id, v ) ) ;
153 }
154 void accumulateFlag( uInt &id, Vector<Bool> &v )
155 {
156 flag.insert( pair< uInt,Vector<Bool> >( id, v ) ) ;
157 }
158 void accumulateFlagRow( Bool &v )
159 {
160 flagrow |= v ;
161 }
162 void addPol( uInt id )
163 {
164 uInt i = polnos.nelements() ;
165 polnos.resize( i+1, True ) ;
166 polnos[i] = id ;
167 }
168
169 uInt nchan;
170 const uInt maxnpol;
171 uInt npol;
172 Vector<uInt> polnos;
173
174 map< uInt,Vector<Float> > data;
175 map< uInt,Vector<Bool> > flag;
176 Bool flagrow;
177};
178
179class BaseMSWriterVisitor: public TableVisitor {
180 const String *lastFieldName;
181 uInt lastRecordNo;
182 uInt lastBeamNo, lastScanNo, lastIfNo, lastPolNo;
183 Int lastSrcType;
184 uInt lastCycleNo;
185 Double lastTime;
186protected:
187 const Table &table;
188 uInt count;
189public:
190 BaseMSWriterVisitor(const Table &table)
191 : table(table)
192 {
193 static const String dummy;
194 lastFieldName = &dummy;
195 count = 0;
196 }
197
198 virtual void enterFieldName(const uInt recordNo, const String &columnValue) {
199 }
200 virtual void leaveFieldName(const uInt recordNo, const String &columnValue) {
201 }
202 virtual void enterBeamNo(const uInt recordNo, uInt columnValue) { }
203 virtual void leaveBeamNo(const uInt recordNo, uInt columnValue) { }
204 virtual void enterScanNo(const uInt recordNo, uInt columnValue) { }
205 virtual void leaveScanNo(const uInt recordNo, uInt columnValue) { }
206 virtual void enterIfNo(const uInt recordNo, uInt columnValue) { }
207 virtual void leaveIfNo(const uInt recordNo, uInt columnValue) { }
208 virtual void enterSrcType(const uInt recordNo, Int columnValue) { }
209 virtual void leaveSrcType(const uInt recordNo, Int columnValue) { }
210 virtual void enterCycleNo(const uInt recordNo, uInt columnValue) { }
211 virtual void leaveCycleNo(const uInt recordNo, uInt columnValue) { }
212 virtual void enterTime(const uInt recordNo, Double columnValue) { }
213 virtual void leaveTime(const uInt recordNo, Double columnValue) { }
214 virtual void enterPolNo(const uInt recordNo, uInt columnValue) { }
215 virtual void leavePolNo(const uInt recordNo, uInt columnValue) { }
216
217 virtual Bool visitRecord(const uInt recordNo,
218 const String &fieldName,
219 const uInt beamNo,
220 const uInt scanNo,
221 const uInt ifNo,
222 const Int srcType,
223 const uInt cycleNo,
224 const Double time,
225 const uInt polNo) { return True ;}
226
227 virtual Bool visit(Bool isFirst, const uInt recordNo,
228 const uInt nCols, void const *const colValues[]) {
229 const String *fieldName = NULL;
230 uInt beamNo, scanNo, ifNo;
231 Int srcType;
232 uInt cycleNo;
233 Double time;
234 uInt polNo;
235 { // prologue
236 uInt i = 0;
237 {
238 const String *col = (const String*)colValues[i++];
239 fieldName = &col[recordNo];
240 }
241 {
242 const uInt *col = (const uInt *)colValues[i++];
243 beamNo = col[recordNo];
244 }
245 {
246 const uInt *col = (const uInt *)colValues[i++];
247 scanNo = col[recordNo];
248 }
249 {
250 const uInt *col = (const uInt *)colValues[i++];
251 ifNo = col[recordNo];
252 }
253 {
254 const Int *col = (const Int *)colValues[i++];
255 srcType = col[recordNo];
256 }
257 {
258 const uInt *col = (const uInt *)colValues[i++];
259 cycleNo = col[recordNo];
260 }
261 {
262 const Double *col = (const Double *)colValues[i++];
263 time = col[recordNo];
264 }
265 {
266 const Int *col = (const Int *)colValues[i++];
267 polNo = col[recordNo];
268 }
269 assert(nCols == i);
270 }
271
272 if (isFirst) {
273 enterFieldName(recordNo, *fieldName);
274 enterBeamNo(recordNo, beamNo);
275 enterScanNo(recordNo, scanNo);
276 enterIfNo(recordNo, ifNo);
277 enterSrcType(recordNo, srcType);
278 enterCycleNo(recordNo, cycleNo);
279 enterTime(recordNo, time);
280 enterPolNo(recordNo, polNo);
281 } else {
282 if (lastFieldName->compare(*fieldName) != 0) {
283 leavePolNo(lastRecordNo, lastPolNo);
284 leaveTime(lastRecordNo, lastTime);
285 leaveCycleNo(lastRecordNo, lastCycleNo);
286 leaveSrcType(lastRecordNo, lastSrcType);
287 leaveIfNo(lastRecordNo, lastIfNo);
288 leaveScanNo(lastRecordNo, lastScanNo);
289 leaveBeamNo(lastRecordNo, lastBeamNo);
290 leaveFieldName(lastRecordNo, *lastFieldName);
291
292 enterFieldName(recordNo, *fieldName);
293 enterBeamNo(recordNo, beamNo);
294 enterScanNo(recordNo, scanNo);
295 enterIfNo(recordNo, ifNo);
296 enterSrcType(recordNo, srcType);
297 enterCycleNo(recordNo, cycleNo);
298 enterTime(recordNo, time);
299 enterPolNo(recordNo, polNo);
300 } else if (lastBeamNo != beamNo) {
301 leavePolNo(lastRecordNo, lastPolNo);
302 leaveTime(lastRecordNo, lastTime);
303 leaveCycleNo(lastRecordNo, lastCycleNo);
304 leaveSrcType(lastRecordNo, lastSrcType);
305 leaveIfNo(lastRecordNo, lastIfNo);
306 leaveScanNo(lastRecordNo, lastScanNo);
307 leaveBeamNo(lastRecordNo, lastBeamNo);
308
309 enterBeamNo(recordNo, beamNo);
310 enterScanNo(recordNo, scanNo);
311 enterIfNo(recordNo, ifNo);
312 enterSrcType(recordNo, srcType);
313 enterCycleNo(recordNo, cycleNo);
314 enterTime(recordNo, time);
315 enterPolNo(recordNo, polNo);
316 } else if (lastScanNo != scanNo) {
317 leavePolNo(lastRecordNo, lastPolNo);
318 leaveTime(lastRecordNo, lastTime);
319 leaveCycleNo(lastRecordNo, lastCycleNo);
320 leaveSrcType(lastRecordNo, lastSrcType);
321 leaveIfNo(lastRecordNo, lastIfNo);
322 leaveScanNo(lastRecordNo, lastScanNo);
323
324 enterScanNo(recordNo, scanNo);
325 enterIfNo(recordNo, ifNo);
326 enterSrcType(recordNo, srcType);
327 enterCycleNo(recordNo, cycleNo);
328 enterTime(recordNo, time);
329 enterPolNo(recordNo, polNo);
330 } else if (lastIfNo != ifNo) {
331 leavePolNo(lastRecordNo, lastPolNo);
332 leaveTime(lastRecordNo, lastTime);
333 leaveCycleNo(lastRecordNo, lastCycleNo);
334 leaveSrcType(lastRecordNo, lastSrcType);
335 leaveIfNo(lastRecordNo, lastIfNo);
336
337 enterIfNo(recordNo, ifNo);
338 enterSrcType(recordNo, srcType);
339 enterCycleNo(recordNo, cycleNo);
340 enterTime(recordNo, time);
341 enterPolNo(recordNo, polNo);
342 } else if (lastSrcType != srcType) {
343 leavePolNo(lastRecordNo, lastPolNo);
344 leaveTime(lastRecordNo, lastTime);
345 leaveCycleNo(lastRecordNo, lastCycleNo);
346 leaveSrcType(lastRecordNo, lastSrcType);
347
348 enterSrcType(recordNo, srcType);
349 enterCycleNo(recordNo, cycleNo);
350 enterTime(recordNo, time);
351 enterPolNo(recordNo, polNo);
352 } else if (lastCycleNo != cycleNo) {
353 leavePolNo(lastRecordNo, lastPolNo);
354 leaveTime(lastRecordNo, lastTime);
355 leaveCycleNo(lastRecordNo, lastCycleNo);
356
357 enterCycleNo(recordNo, cycleNo);
358 enterTime(recordNo, time);
359 enterPolNo(recordNo, polNo);
360 } else if (lastTime != time) {
361 leavePolNo(lastRecordNo, lastPolNo);
362 leaveTime(lastRecordNo, lastTime);
363
364 enterTime(recordNo, time);
365 enterPolNo(recordNo, polNo);
366 } else if (lastPolNo != polNo) {
367 leavePolNo(lastRecordNo, lastPolNo);
368 enterPolNo(recordNo, polNo);
369 }
370 }
371 count++;
372 Bool result = visitRecord(recordNo, *fieldName, beamNo, scanNo, ifNo, srcType,
373 cycleNo, time, polNo);
374
375 { // epilogue
376 lastRecordNo = recordNo;
377
378 lastFieldName = fieldName;
379 lastBeamNo = beamNo;
380 lastScanNo = scanNo;
381 lastIfNo = ifNo;
382 lastSrcType = srcType;
383 lastCycleNo = cycleNo;
384 lastTime = time;
385 lastPolNo = polNo;
386 }
387 return result ;
388 }
389
390 virtual void finish() {
391 if (count > 0) {
392 leavePolNo(lastRecordNo, lastPolNo);
393 leaveTime(lastRecordNo, lastTime);
394 leaveCycleNo(lastRecordNo, lastCycleNo);
395 leaveSrcType(lastRecordNo, lastSrcType);
396 leaveIfNo(lastRecordNo, lastIfNo);
397 leaveScanNo(lastRecordNo, lastScanNo);
398 leaveBeamNo(lastRecordNo, lastBeamNo);
399 leaveFieldName(lastRecordNo, *lastFieldName);
400 }
401 }
402};
403
404class MSWriterVisitor: public BaseMSWriterVisitor, public MSWriterUtils {
405public:
406 MSWriterVisitor(const Table &table, Table &mstable)
407 : BaseMSWriterVisitor(table),
408 ms(mstable)
409 {
410 rowidx = 0 ;
411 fieldName = "" ;
412 defaultFieldId = 0 ;
413 spwId = -1 ;
414 subscan = 1 ;
415 ptName = "" ;
416 srcId = 0 ;
417
418 row = TableRow( ms ) ;
419
420 holder.reset() ;
421
422 makePolMap() ;
423 initFrequencies() ;
424 initCorrProductTemplate() ;
425
426 //
427 // add rows to MS
428 //
429 uInt addrow = table.nrow() ;
430 ms.addRow( addrow ) ;
431
432 // attach to Scantable columns
433 spectraCol.attach( table, "SPECTRA" ) ;
434 flagtraCol.attach( table, "FLAGTRA" ) ;
435 flagRowCol.attach( table, "FLAGROW" ) ;
436 tcalIdCol.attach( table, "TCAL_ID" ) ;
437 intervalCol.attach( table, "INTERVAL" ) ;
438 directionCol.attach( table, "DIRECTION" ) ;
439 scanRateCol.attach( table, "SCANRATE" ) ;
440 timeCol.attach( table, "TIME" ) ;
441 freqIdCol.attach( table, "FREQ_ID" ) ;
442 sourceNameCol.attach( table, "SRCNAME" ) ;
443 sourceDirectionCol.attach( table, "SRCDIRECTION" ) ;
444 fieldNameCol.attach( table, "FIELDNAME" ) ;
445
446 // MS subtables
447 attachSubtables() ;
448
449 // attach to MS columns
450 attachMain() ;
451 attachPointing() ;
452 }
453
454 virtual void enterFieldName(const uInt recordNo, const String &columnValue) {
455 //printf("%u: FieldName: %s\n", recordNo, columnValue.c_str());
456 fieldName = fieldNameCol.asString( recordNo ) ;
457 String::size_type pos = fieldName.find( "__" ) ;
458 if ( pos != String::npos ) {
459 fieldId = String::toInt( fieldName.substr( pos+2 ) ) ;
460 fieldName = fieldName.substr( 0, pos ) ;
461 }
462 else {
463 fieldId = defaultFieldId ;
464 defaultFieldId++ ;
465 }
466 Double tSec = timeCol.asdouble( recordNo ) * 86400.0 ;
467 Vector<Double> srcDir = sourceDirectionCol( recordNo ) ;
468 Vector<Double> srate = scanRateCol( recordNo ) ;
469 String srcName = sourceNameCol.asString( recordNo ) ;
470
471 addField( fieldId, fieldName, srcName, srcDir, srate, tSec ) ;
472
473 // put value
474 *fieldIdRF = fieldId ;
475 }
476 virtual void leaveFieldName(const uInt recordNo, const String &columnValue) {
477 }
478 virtual void enterBeamNo(const uInt recordNo, uInt columnValue) {
479 //printf("%u: BeamNo: %u\n", recordNo, columnValue);
480
481 feedId = (Int)columnValue ;
482
483 // put value
484 *feed1RF = feedId ;
485 *feed2RF = feedId ;
486 }
487 virtual void leaveBeamNo(const uInt recordNo, uInt columnValue) {
488 }
489 virtual void enterScanNo(const uInt recordNo, uInt columnValue) {
490 //printf("%u: ScanNo: %u\n", recordNo, columnValue);
491
492 // put value
493 // SCAN_NUMBER is 0-based in Scantable while 1-based in MS
494 *scanNumberRF = (Int)columnValue + 1 ;
495 }
496 virtual void leaveScanNo(const uInt recordNo, uInt columnValue) {
497 subscan = 1 ;
498 }
499 virtual void enterIfNo(const uInt recordNo, uInt columnValue) {
500 //printf("%u: IfNo: %u\n", recordNo, columnValue);
501
502 spwId = (Int)columnValue ;
503 uInt freqId = freqIdCol.asuInt( recordNo ) ;
504
505 Vector<Float> sp = spectraCol( recordNo ) ;
506 uInt nchan = sp.nelements() ;
507 holder.setNchan( nchan ) ;
508
509 addSpectralWindow( spwId, freqId ) ;
510
511 addFeed( feedId, spwId ) ;
512 }
513 virtual void leaveIfNo(const uInt recordNo, uInt columnValue) {
514 }
515 virtual void enterSrcType(const uInt recordNo, Int columnValue) {
516 //printf("%u: SrcType: %d\n", recordNo, columnValue);
517
518 Int stateId = addState( columnValue ) ;
519
520 // put value
521 *stateIdRF = stateId ;
522 }
523 virtual void leaveSrcType(const uInt recordNo, Int columnValue) {
524 }
525 virtual void enterCycleNo(const uInt recordNo, uInt columnValue) {
526 //printf("%u: CycleNo: %u\n", recordNo, columnValue);
527 }
528 virtual void leaveCycleNo(const uInt recordNo, uInt columnValue) {
529 }
530 virtual void enterTime(const uInt recordNo, Double columnValue) {
531 //printf("%u: Time: %f\n", recordNo, columnValue);
532
533 Double timeSec = columnValue * 86400.0 ;
534 Double interval = intervalCol.asdouble( recordNo ) ;
535
536 if ( ptName.empty() ) {
537 Vector<Double> dir = directionCol( recordNo ) ;
538 Vector<Double> rate = scanRateCol( recordNo ) ;
539 if ( anyNE( rate, 0.0 ) ) {
540 Matrix<Double> msdir( 2, 2 ) ;
541 msdir.column( 0 ) = dir ;
542 msdir.column( 1 ) = rate ;
543 addPointing( timeSec, interval, msdir ) ;
544 }
545 else {
546 Matrix<Double> msdir( 2, 1 ) ;
547 msdir.column( 0 ) = dir ;
548 addPointing( timeSec, interval, msdir ) ;
549 }
550 }
551
552 // put value
553 *timeRF = timeSec ;
554 *timeCentroidRF = timeSec ;
555 *intervalRF = interval ;
556 *exposureRF = interval ;
557 }
558 virtual void leaveTime(const uInt recordNo, Double columnValue) {
559 if ( holder.nPol() > 0 ) {
560 Vector<Float> w = holder.getWeight() ;
561 Cube<Bool> c = holder.getFlagCategory() ;
562 Bool flr = holder.getFlagRow() ;
563 Matrix<Bool> fl = holder.getFlag() ;
564 Vector<uInt> polnos = holder.polNos() ;
565 Int polId = addPolarization( polnos ) ;
566 Int ddId = addDataDescription( polId, spwId ) ;
567
568 // put field
569 *dataDescIdRF = ddId ;
570 *flagRowRF = flr ;
571 weightRF.define( w ) ;
572 sigmaRF.define( w ) ;
573 flagCategoryRF.define( c ) ;
574 flagRF.define( fl ) ;
575 if ( useFloat ) {
576 Matrix<Float> sp = holder.getData() ;
577 floatDataRF.define( sp ) ;
578 }
579 else {
580 Matrix<Complex> sp = holder.getComplexData() ;
581 dataRF.define( sp ) ;
582 }
583
584 // commit row
585 row.put( rowidx ) ;
586 rowidx++ ;
587
588 // reset holder
589 holder.reset() ;
590 }
591 }
592 virtual void enterPolNo(const uInt recordNo, uInt columnValue) {
593 //printf("%u: PolNo: %d\n", recordNo, columnValue);
594 }
595 virtual void leavePolNo(const uInt recordNo, uInt columnValue) {
596 }
597
598 virtual Bool visitRecord(const uInt recordNo,
599 const String &fieldName,
600 const uInt beamNo,
601 const uInt scanNo,
602 const uInt ifNo,
603 const Int srcType,
604 const uInt cycleNo,
605 const Double time,
606 const uInt polNo) {
607 //printf("%u: %s, %u, %u, %u, %d, %u, %f, %d\n", recordNo,
608 // fieldName.c_str(), beamNo, scanNo, ifNo, srcType, cycleNo, time, polNo);
609
610 Vector<Float> sp = spectraCol( recordNo ) ;
611 Vector<uChar> tmp = flagtraCol( recordNo ) ;
612 Vector<Bool> fl( tmp.shape() ) ;
613 convertArray( fl, tmp ) ;
614 Bool flr = (Bool)flagRowCol.asuInt( recordNo ) ;
615 holder.accumulate( polNo, sp, fl, flr ) ;
616
617 return True ;
618 }
619
620 virtual void finish() {
621 BaseMSWriterVisitor::finish();
622 //printf("Total: %u\n", count);
623
624 // remove rows
625 if ( ms.nrow() > rowidx ) {
626 uInt numRemove = ms.nrow() - rowidx ;
627 //cout << "numRemove = " << numRemove << endl ;
628 Vector<uInt> rows( numRemove ) ;
629 indgen( rows, rowidx ) ;
630 ms.removeRow( rows ) ;
631 }
632
633 // fill empty SPECTRAL_WINDOW rows
634 infillSpectralWindow() ;
635 }
636
637 void dataColumnName( String name )
638 {
639 TableRecord &r = row.record() ;
640 if ( name == "DATA" ) {
641 useFloat = False ;
642 dataRF.attachToRecord( r, name ) ;
643 }
644 else if ( name == "FLOAT_DATA" ) {
645 useFloat = True ;
646 floatDataRF.attachToRecord( r, name ) ;
647 }
648 }
649 void pointingTableName( String name ) {
650 ptName = name ;
651 }
652 void setSourceRecord( Record &r ) {
653 srcRec = r ;
654 }
655private:
656 void addField( Int &fid, String &fname, String &srcName,
657 Vector<Double> &sdir, Vector<Double> &srate,
658 Double &tSec )
659 {
660 uInt nrow = fieldtab.nrow() ;
661 while( (Int)nrow <= fid ) {
662 fieldtab.addRow( 1, True ) ;
663 nrow++ ;
664 }
665
666 Matrix<Double> dir ;
667 Int numPoly = 0 ;
668 if ( anyNE( srate, 0.0 ) ) {
669 dir.resize( 2, 2 ) ;
670 dir.column( 0 ) = sdir ;
671 dir.column( 1 ) = srate ;
672 numPoly = 1 ;
673 }
674 else {
675 dir.resize( 2, 1 ) ;
676 dir.column( 0 ) = sdir ;
677 }
678 srcId = srcRec.asInt( srcName ) ;
679
680 TableRow tr( fieldtab ) ;
681 TableRecord &r = tr.record() ;
682 putField( "NAME", r, fname ) ;
683 putField( "NUM_POLY", r, numPoly ) ;
684 putField( "TIME", r, tSec ) ;
685 putField( "SOURCE_ID", r, srcId ) ;
686 defineField( "DELAY_DIR", r, dir ) ;
687 defineField( "REFERENCE_DIR", r, dir ) ;
688 defineField( "PHASE_DIR", r, dir ) ;
689 tr.put( fid ) ;
690
691 // for POINTING table
692 *poNameRF = fname ;
693 }
694 Int addState( Int &id )
695 {
696 String obsMode ;
697 Bool isSignal ;
698 Double tnoise ;
699 Double tload ;
700 queryType( id, obsMode, isSignal, tnoise, tload ) ;
701
702 String key = obsMode+"_"+String::toString( subscan ) ;
703 Int idx = -1 ;
704 uInt nEntry = stateEntry.nelements() ;
705 for ( uInt i = 0 ; i < nEntry ; i++ ) {
706 if ( stateEntry[i] == key ) {
707 idx = i ;
708 break ;
709 }
710 }
711 if ( idx == -1 ) {
712 uInt nrow = statetab.nrow() ;
713 statetab.addRow( 1, True ) ;
714 TableRow tr( statetab ) ;
715 TableRecord &r = tr.record() ;
716 putField( "OBS_MODE", r, obsMode ) ;
717 putField( "SIG", r, isSignal ) ;
718 isSignal = !isSignal ;
719 putField( "REF", r, isSignal ) ;
720 putField( "CAL", r, tnoise ) ;
721 putField( "LOAD", r, tload ) ;
722 tr.put( nrow ) ;
723 idx = nrow ;
724
725 stateEntry.resize( nEntry+1, True ) ;
726 stateEntry[nEntry] = key ;
727 }
728 subscan++ ;
729
730 return idx ;
731 }
732 void addPointing( Double &tSec, Double &interval, Matrix<Double> &dir )
733 {
734 uInt nrow = potab.nrow() ;
735 potab.addRow( 1, True ) ;
736
737 *poNumPolyRF = dir.ncolumn() - 1 ;
738 *poTimeRF = tSec ;
739 *poTimeOriginRF = tSec ;
740 *poIntervalRF = interval ;
741 poDirectionRF.define( dir ) ;
742 poTargetRF.define( dir ) ;
743 porow.put( nrow ) ;
744 }
745 Int addPolarization( Vector<uInt> &nos )
746 {
747 Int idx = -1 ;
748 uInt nEntry = polEntry.size() ;
749 for ( uInt i = 0 ; i < nEntry ; i++ ) {
750 if ( polEntry[i].conform( nos ) && allEQ( polEntry[i], nos ) ) {
751 idx = i ;
752 break ;
753 }
754 }
755
756 Int numCorr ;
757 Vector<Int> corrType ;
758 Matrix<Int> corrProduct ;
759 polProperty( nos, numCorr, corrType, corrProduct ) ;
760
761 if ( idx == -1 ) {
762 uInt nrow = poltab.nrow() ;
763 poltab.addRow( 1, True ) ;
764 TableRow tr( poltab ) ;
765 TableRecord &r = tr.record() ;
766 putField( "NUM_CORR", r, numCorr ) ;
767 defineField( "CORR_TYPE", r, corrType ) ;
768 defineField( "CORR_PRODUCT", r, corrProduct ) ;
769 tr.put( nrow ) ;
770 idx = nrow ;
771
772 polEntry.resize( nEntry+1 ) ;
773 polEntry[nEntry] = nos ;
774 }
775
776 return idx ;
777 }
778 Int addDataDescription( Int pid, Int sid )
779 {
780 Int idx = -1 ;
781 uInt nEntry = ddEntry.nrow() ;
782 Vector<Int> key( 2 ) ;
783 key[0] = pid ;
784 key[1] = sid ;
785 for ( uInt i = 0 ; i < nEntry ; i++ ) {
786 if ( allEQ( ddEntry.row(i), key ) ) {
787 idx = i ;
788 break ;
789 }
790 }
791
792 if ( idx == -1 ) {
793 uInt nrow = ddtab.nrow() ;
794 ddtab.addRow( 1, True ) ;
795 TableRow tr( ddtab ) ;
796 TableRecord &r = tr.record() ;
797 putField( "POLARIZATION_ID", r, pid ) ;
798 putField( "SPECTRAL_WINDOW_ID", r, sid ) ;
799 tr.put( nrow ) ;
800 idx = nrow ;
801
802 ddEntry.resize( nEntry+1, 2, True ) ;
803 ddEntry.row(nEntry) = key ;
804 }
805
806 return idx ;
807 }
808 void infillSpectralWindow()
809 {
810 ROScalarColumn<Int> nchanCol( spwtab, "NUM_CHAN" ) ;
811 Vector<Int> nchan = nchanCol.getColumn() ;
812 TableRow tr( spwtab ) ;
813 TableRecord &r = tr.record() ;
814 Int mfr = 1 ;
815 Vector<Double> dummy( 1, 0.0 ) ;
816 putField( "MEAS_FREQ_REF", r, mfr ) ;
817 defineField( "CHAN_FREQ", r, dummy ) ;
818 defineField( "CHAN_WIDTH", r, dummy ) ;
819 defineField( "EFFECTIVE_BW", r, dummy ) ;
820 defineField( "RESOLUTION", r, dummy ) ;
821
822 for ( uInt i = 0 ; i < spwtab.nrow() ; i++ ) {
823 if ( nchan[i] == 0 )
824 tr.put( i ) ;
825 }
826 }
827 void addSpectralWindow( Int sid, uInt fid )
828 {
829 if ( !processedFreqId[fid] ) {
830 uInt nrow = spwtab.nrow() ;
831 while( (Int)nrow <= sid ) {
832 spwtab.addRow( 1, True ) ;
833 nrow++ ;
834 }
835 processedFreqId[fid] = True ;
836 }
837
838 Double rp = refpix[fid] ;
839 Double rv = refval[fid] ;
840 Double ic = increment[fid] ;
841
842 Int mfrInt = (Int)freqframe ;
843 Int nchan = holder.nChan() ;
844 Double bw = nchan * abs( ic ) ;
845 Double reffreq = rv - rp * ic ;
846 Int netsb = 0 ; // USB->0, LSB->1
847 if ( ic < 0 )
848 netsb = 1 ;
849 Vector<Double> res( nchan, abs(ic) ) ;
850 Vector<Double> chanf( nchan ) ;
851 indgen( chanf, reffreq, ic ) ;
852
853 TableRow tr( spwtab ) ;
854 TableRecord &r = tr.record() ;
855 putField( "MEAS_FREQ_REF", r, mfrInt ) ;
856 putField( "NUM_CHAN", r, nchan ) ;
857 putField( "TOTAL_BANDWIDTH", r, bw ) ;
858 putField( "REF_FREQUENCY", r, reffreq ) ;
859 putField( "NET_SIDEBAND", r, netsb ) ;
860 defineField( "RESOLUTION", r, res ) ;
861 defineField( "CHAN_WIDTH", r, res ) ;
862 defineField( "EFFECTIVE_BW", r, res ) ;
863 defineField( "CHAN_FREQ", r, chanf ) ;
864 tr.put( sid ) ;
865 }
866 void addFeed( Int fid, Int sid )
867 {
868 Int idx = -1 ;
869 uInt nEntry = feedEntry.nrow() ;
870 Vector<Int> key( 2 ) ;
871 key[0] = fid ;
872 key[1] = sid ;
873 for ( uInt i = 0 ; i < nEntry ; i++ ) {
874 if ( allEQ( feedEntry.row(i), key ) ) {
875 idx = i ;
876 break ;
877 }
878 }
879
880 if ( idx == -1 ) {
881 uInt nrow = feedtab.nrow() ;
882 feedtab.addRow( 1, True ) ;
883 Int numReceptors = 2 ;
884 Vector<String> polType( numReceptors ) ;
885 Matrix<Double> beamOffset( 2, numReceptors, 0.0 ) ;
886 Vector<Double> receptorAngle( numReceptors, 0.0 ) ;
887 if ( poltype == "linear" ) {
888 polType[0] = "X" ;
889 polType[1] = "Y" ;
890 }
891 else if ( poltype == "circular" ) {
892 polType[0] = "R" ;
893 polType[1] = "L" ;
894 }
895 else {
896 polType[0] = "X" ;
897 polType[1] = "Y" ;
898 }
899 Matrix<Complex> polResponse( numReceptors, numReceptors, 0.0 ) ;
900
901 TableRow tr( feedtab ) ;
902 TableRecord &r = tr.record() ;
903 putField( "FEED_ID", r, fid ) ;
904 putField( "BEAM_ID", r, fid ) ;
905 Int tmp = 0 ;
906 putField( "ANTENNA_ID", r, tmp ) ;
907 putField( "SPECTRAL_WINDOW_ID", r, sid ) ;
908 putField( "NUM_RECEPTORS", r, numReceptors ) ;
909 defineField( "POLARIZATION_TYPE", r, polType ) ;
910 defineField( "BEAM_OFFSET", r, beamOffset ) ;
911 defineField( "RECEPTOR_ANGLE", r, receptorAngle ) ;
912 defineField( "POL_RESPONSE", r, polResponse ) ;
913 tr.put( nrow ) ;
914
915 feedEntry.resize( nEntry+1, 2, True ) ;
916 feedEntry.row( nEntry ) = key ;
917 }
918 }
919 void makePolMap()
920 {
921 const TableRecord &keys = table.keywordSet() ;
922 poltype = keys.asString( "POLTYPE" ) ;
923
924 if ( poltype == "stokes" ) {
925 polmap.resize( 4 ) ;
926 polmap[0] = Stokes::I ;
927 polmap[1] = Stokes::Q ;
928 polmap[2] = Stokes::U ;
929 polmap[3] = Stokes::V ;
930 }
931 else if ( poltype == "linear" ) {
932 polmap.resize( 4 ) ;
933 polmap[0] = Stokes::XX ;
934 polmap[1] = Stokes::YY ;
935 polmap[2] = Stokes::XY ;
936 polmap[3] = Stokes::YX ;
937 }
938 else if ( poltype == "circular" ) {
939 polmap.resize( 4 ) ;
940 polmap[0] = Stokes::RR ;
941 polmap[1] = Stokes::LL ;
942 polmap[2] = Stokes::RL ;
943 polmap[3] = Stokes::LR ;
944 }
945 else if ( poltype == "linpol" ) {
946 polmap.resize( 2 ) ;
947 polmap[0] = Stokes::Plinear ;
948 polmap[1] = Stokes::Pangle ;
949 }
950 else {
951 polmap.resize( 0 ) ;
952 }
953 }
954 void initFrequencies()
955 {
956 const TableRecord &keys = table.keywordSet() ;
957 Table tab = keys.asTable( "FREQUENCIES" ) ;
958 ROScalarColumn<uInt> idcol( tab, "ID" ) ;
959 ROScalarColumn<Double> rpcol( tab, "REFPIX" ) ;
960 ROScalarColumn<Double> rvcol( tab, "REFVAL" ) ;
961 ROScalarColumn<Double> iccol( tab, "INCREMENT" ) ;
962 Vector<uInt> id = idcol.getColumn() ;
963 Vector<Double> rp = rpcol.getColumn() ;
964 Vector<Double> rv = rvcol.getColumn() ;
965 Vector<Double> ic = iccol.getColumn() ;
966 for ( uInt i = 0 ; i < id.nelements() ; i++ ) {
967 processedFreqId.insert( pair<uInt,Bool>( id[i], False ) ) ;
968 refpix.insert( pair<uInt,Double>( id[i], rp[i] ) ) ;
969 refval.insert( pair<uInt,Double>( id[i], rv[i] ) ) ;
970 increment.insert( pair<uInt,Double>( id[i], ic[i] ) ) ;
971 }
972 String frameStr = tab.keywordSet().asString( "BASEFRAME" ) ;
973 MFrequency::getType( freqframe, frameStr ) ;
974 }
975 void attachSubtables()
976 {
977 const TableRecord &keys = table.keywordSet() ;
978 TableRecord &mskeys = ms.rwKeywordSet() ;
979
980 // FIELD table
981 fieldtab = mskeys.asTable( "FIELD" ) ;
982
983 // SPECTRAL_WINDOW table
984 spwtab = mskeys.asTable( "SPECTRAL_WINDOW" ) ;
985
986 // POINTING table
987 potab = mskeys.asTable( "POINTING" ) ;
988
989 // POLARIZATION table
990 poltab = mskeys.asTable( "POLARIZATION" ) ;
991
992 // DATA_DESCRIPTION table
993 ddtab = mskeys.asTable( "DATA_DESCRIPTION" ) ;
994
995 // STATE table
996 statetab = mskeys.asTable( "STATE" ) ;
997
998 // FEED table
999 feedtab = mskeys.asTable( "FEED" ) ;
1000 }
1001 void attachMain()
1002 {
1003 TableRecord &r = row.record() ;
1004 dataDescIdRF.attachToRecord( r, "DATA_DESC_ID" ) ;
1005 flagRowRF.attachToRecord( r, "FLAG_ROW" ) ;
1006 weightRF.attachToRecord( r, "WEIGHT" ) ;
1007 sigmaRF.attachToRecord( r, "SIGMA" ) ;
1008 flagCategoryRF.attachToRecord( r, "FLAG_CATEGORY" ) ;
1009 flagRF.attachToRecord( r, "FLAG" ) ;
1010 timeRF.attachToRecord( r, "TIME" ) ;
1011 timeCentroidRF.attachToRecord( r, "TIME_CENTROID" ) ;
1012 intervalRF.attachToRecord( r, "INTERVAL" ) ;
1013 exposureRF.attachToRecord( r, "EXPOSURE" ) ;
1014 fieldIdRF.attachToRecord( r, "FIELD_ID" ) ;
1015 feed1RF.attachToRecord( r, "FEED1" ) ;
1016 feed2RF.attachToRecord( r, "FEED2" ) ;
1017 scanNumberRF.attachToRecord( r, "SCAN_NUMBER" ) ;
1018 stateIdRF.attachToRecord( r, "STATE_ID" ) ;
1019
1020 // constant values
1021 Int id = 0 ;
1022 RecordFieldPtr<Int> intRF( r, "OBSERVATION_ID" ) ;
1023 *intRF = 0 ;
1024 intRF.attachToRecord( r, "ANTENNA1" ) ;
1025 *intRF = 0 ;
1026 intRF.attachToRecord( r, "ANTENNA2" ) ;
1027 *intRF = 0 ;
1028 intRF.attachToRecord( r, "ARRAY_ID" ) ;
1029 *intRF = 0 ;
1030 intRF.attachToRecord( r, "PROCESSOR_ID" ) ;
1031 *intRF = 0 ;
1032 RecordFieldPtr< Vector<Double> > arrayRF( r, "UVW" ) ;
1033 arrayRF.define( Vector<Double>( 3, 0.0 ) ) ;
1034 }
1035 void attachPointing()
1036 {
1037 porow = TableRow( potab ) ;
1038 TableRecord &r = porow.record() ;
1039 poNumPolyRF.attachToRecord( r, "NUM_POLY" ) ;
1040 poTimeRF.attachToRecord( r, "TIME" ) ;
1041 poTimeOriginRF.attachToRecord( r, "TIME_ORIGIN" ) ;
1042 poIntervalRF.attachToRecord( r, "INTERVAL" ) ;
1043 poNameRF.attachToRecord( r, "NAME" ) ;
1044 poDirectionRF.attachToRecord( r, "DIRECTION" ) ;
1045 poTargetRF.attachToRecord( r, "TARGET" ) ;
1046
1047 // constant values
1048 RecordFieldPtr<Int> antIdRF( r, "ANTENNA_ID" ) ;
1049 *antIdRF = 0 ;
1050 RecordFieldPtr<Bool> trackingRF( r, "TRACKING" ) ;
1051 *trackingRF = True ;
1052 }
1053 void queryType( Int type, String &stype, Bool &b, Double &t, Double &l )
1054 {
1055 t = 0.0 ;
1056 l = 0.0 ;
1057
1058 String sep1="#" ;
1059 String sep2="," ;
1060 String target="OBSERVE_TARGET" ;
1061 String atmcal="CALIBRATE_TEMPERATURE" ;
1062 String onstr="ON_SOURCE" ;
1063 String offstr="OFF_SOURCE" ;
1064 String pswitch="POSITION_SWITCH" ;
1065 String nod="NOD" ;
1066 String fswitch="FREQUENCY_SWITCH" ;
1067 String sigstr="SIG" ;
1068 String refstr="REF" ;
1069 String unspecified="UNSPECIFIED" ;
1070 String ftlow="LOWER" ;
1071 String fthigh="HIGHER" ;
1072 switch ( type ) {
1073 case SrcType::PSON:
1074 stype = target+sep1+onstr+sep2+pswitch ;
1075 b = True ;
1076 break ;
1077 case SrcType::PSOFF:
1078 stype = target+sep1+offstr+sep2+pswitch ;
1079 b = False ;
1080 break ;
1081 case SrcType::NOD:
1082 stype = target+sep1+onstr+sep2+nod ;
1083 b = True ;
1084 break ;
1085 case SrcType::FSON:
1086 stype = target+sep1+onstr+sep2+fswitch+sep1+sigstr ;
1087 b = True ;
1088 break ;
1089 case SrcType::FSOFF:
1090 stype = target+sep1+onstr+sep2+fswitch+sep1+refstr ;
1091 b = False ;
1092 break ;
1093 case SrcType::SKY:
1094 stype = atmcal+sep1+offstr+sep2+unspecified ;
1095 b = False ;
1096 break ;
1097 case SrcType::HOT:
1098 stype = atmcal+sep1+offstr+sep2+unspecified ;
1099 b = False ;
1100 break ;
1101 case SrcType::WARM:
1102 stype = atmcal+sep1+offstr+sep2+unspecified ;
1103 b = False ;
1104 break ;
1105 case SrcType::COLD:
1106 stype = atmcal+sep1+offstr+sep2+unspecified ;
1107 b = False ;
1108 break ;
1109 case SrcType::PONCAL:
1110 stype = atmcal+sep1+onstr+sep2+pswitch ;
1111 b = True ;
1112 break ;
1113 case SrcType::POFFCAL:
1114 stype = atmcal+sep1+offstr+sep2+pswitch ;
1115 b = False ;
1116 break ;
1117 case SrcType::NODCAL:
1118 stype = atmcal+sep1+onstr+sep2+nod ;
1119 b = True ;
1120 break ;
1121 case SrcType::FONCAL:
1122 stype = atmcal+sep1+onstr+sep2+fswitch+sep1+sigstr ;
1123 b = True ;
1124 break ;
1125 case SrcType::FOFFCAL:
1126 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+refstr ;
1127 b = False ;
1128 break ;
1129 case SrcType::FSLO:
1130 stype = target+sep1+onstr+sep2+fswitch+sep1+ftlow ;
1131 b = True ;
1132 break ;
1133 case SrcType::FLOOFF:
1134 stype = target+sep1+offstr+sep2+fswitch+sep1+ftlow ;
1135 b = False ;
1136 break ;
1137 case SrcType::FLOSKY:
1138 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
1139 b = False ;
1140 break ;
1141 case SrcType::FLOHOT:
1142 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
1143 b = False ;
1144 break ;
1145 case SrcType::FLOWARM:
1146 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
1147 b = False ;
1148 break ;
1149 case SrcType::FLOCOLD:
1150 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
1151 b = False ;
1152 break ;
1153 case SrcType::FSHI:
1154 stype = target+sep1+onstr+sep2+fswitch+sep1+fthigh ;
1155 b = True ;
1156 break ;
1157 case SrcType::FHIOFF:
1158 stype = target+sep1+offstr+sep2+fswitch+sep1+fthigh ;
1159 b = False ;
1160 break ;
1161 case SrcType::FHISKY:
1162 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
1163 b = False ;
1164 break ;
1165 case SrcType::FHIHOT:
1166 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
1167 b = False ;
1168 break ;
1169 case SrcType::FHIWARM:
1170 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
1171 b = False ;
1172 break ;
1173 case SrcType::FHICOLD:
1174 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
1175 b = False ;
1176 break ;
1177 case SrcType::SIG:
1178 stype = target+sep1+onstr+sep2+unspecified ;
1179 b = True ;
1180 break ;
1181 case SrcType::REF:
1182 stype = target+sep1+offstr+sep2+unspecified ;
1183 b = False ;
1184 break ;
1185 default:
1186 stype = unspecified ;
1187 b = True ;
1188 break ;
1189 }
1190 }
1191 void polProperty( Vector<uInt> &nos, Int &n, Vector<Int> &c, Matrix<Int> &cp )
1192 {
1193 n = nos.nelements() ;
1194 c.resize( n ) ;
1195 for ( Int i = 0 ; i < n ; i++ )
1196 c[i] = (Int)polmap[nos[i]] ;
1197 if ( n == 4 &&
1198 ( poltype == "linear" || poltype == "circular" ) ) {
1199 Int tmp = c[1] ;
1200 c[1] = c[2] ;
1201 c[2] = c[3] ;
1202 c[3] = tmp ;
1203 }
1204 cp = corrProductTemplate[n] ;
1205 }
1206 void initCorrProductTemplate()
1207 {
1208 Int n = 1 ;
1209 {
1210 Matrix<Int> c( 2, n, 0 ) ;
1211 corrProductTemplate[n] = c ;
1212 }
1213 n = 2 ;
1214 {
1215 Matrix<Int> c( 2, n, 0 ) ;
1216 c.column( 1 ) = 1 ;
1217 corrProductTemplate[n] = c ;
1218 }
1219 n = 4 ;
1220 {
1221 Matrix<Int> c( 2, n, 0 ) ;
1222 c( 0, 2 ) = 1 ;
1223 c( 0, 3 ) = 1 ;
1224 c( 1, 1 ) = 1 ;
1225 c( 1, 3 ) = 1 ;
1226 corrProductTemplate[n] = c ;
1227 }
1228 }
1229
1230 Table &ms;
1231 TableRow row;
1232 uInt rowidx;
1233 String fieldName;
1234 Int fieldId;
1235 Int srcId;
1236 Int defaultFieldId;
1237 Int spwId;
1238 Int feedId;
1239 Int subscan;
1240 PolarizedComponentHolder holder;
1241 String ptName;
1242 Bool useFloat;
1243 String poltype;
1244 Vector<Stokes::StokesTypes> polmap;
1245
1246 // MS subtables
1247 Table spwtab;
1248 Table statetab;
1249 Table ddtab;
1250 Table poltab;
1251 Table fieldtab;
1252 Table feedtab;
1253 Table potab;
1254
1255 // Scantable MAIN columns
1256 ROArrayColumn<Float> spectraCol;
1257 ROArrayColumn<Double> directionCol,scanRateCol,sourceDirectionCol;
1258 ROArrayColumn<uChar> flagtraCol;
1259 ROTableColumn tcalIdCol,intervalCol,flagRowCol,timeCol,freqIdCol,
1260 sourceNameCol,fieldNameCol;
1261
1262 // MS MAIN columns
1263 RecordFieldPtr<Int> dataDescIdRF,fieldIdRF,feed1RF,feed2RF,
1264 scanNumberRF,stateIdRF;
1265 RecordFieldPtr<Bool> flagRowRF;
1266 RecordFieldPtr<Double> timeRF,timeCentroidRF,intervalRF,exposureRF;
1267 RecordFieldPtr< Vector<Float> > weightRF,sigmaRF;
1268 RecordFieldPtr< Cube<Bool> > flagCategoryRF;
1269 RecordFieldPtr< Matrix<Bool> > flagRF;
1270 RecordFieldPtr< Matrix<Float> > floatDataRF;
1271 RecordFieldPtr< Matrix<Complex> > dataRF;
1272
1273 // MS POINTING columns
1274 TableRow porow;
1275 RecordFieldPtr<Int> poNumPolyRF ;
1276 RecordFieldPtr<Double> poTimeRF,
1277 poTimeOriginRF,
1278 poIntervalRF ;
1279 RecordFieldPtr<String> poNameRF ;
1280 RecordFieldPtr< Matrix<Double> > poDirectionRF,
1281 poTargetRF ;
1282
1283 Vector<String> stateEntry;
1284 Matrix<Int> ddEntry;
1285 Matrix<Int> feedEntry;
1286 vector< Vector<uInt> > polEntry;
1287 map<uInt,Bool> processedFreqId;
1288 map<uInt,Double> refpix;
1289 map<uInt,Double> refval;
1290 map<uInt,Double> increment;
1291 MFrequency::Types freqframe;
1292 Record srcRec;
1293 map< Int, Matrix<Int> > corrProductTemplate;
1294};
1295
1296class BaseMSSysCalVisitor: public TableVisitor {
1297 uInt lastRecordNo;
1298 uInt lastBeamNo, lastIfNo, lastPolNo;
1299 Double lastTime;
1300protected:
1301 const Table &table;
1302 uInt count;
1303public:
1304 BaseMSSysCalVisitor(const Table &table)
1305 : table(table)
1306 {
1307 count = 0;
1308 }
1309
1310 virtual void enterBeamNo(const uInt recordNo, uInt columnValue) { }
1311 virtual void leaveBeamNo(const uInt recordNo, uInt columnValue) { }
1312 virtual void enterIfNo(const uInt recordNo, uInt columnValue) { }
1313 virtual void leaveIfNo(const uInt recordNo, uInt columnValue) { }
1314 virtual void enterPolNo(const uInt recordNo, uInt columnValue) { }
1315 virtual void leavePolNo(const uInt recordNo, uInt columnValue) { }
1316 virtual void enterTime(const uInt recordNo, Double columnValue) { }
1317 virtual void leaveTime(const uInt recordNo, Double columnValue) { }
1318
1319 virtual Bool visitRecord(const uInt recordNo,
1320 const uInt beamNo,
1321 const uInt ifNo,
1322 const uInt polNo,
1323 const Double time) { return True ;}
1324
1325 virtual Bool visit(Bool isFirst, const uInt recordNo,
1326 const uInt nCols, void const *const colValues[]) {
1327 uInt beamNo, ifNo, polNo;
1328 Double time;
1329 { // prologue
1330 uInt i = 0;
1331 {
1332 const uInt *col = (const uInt *)colValues[i++];
1333 beamNo = col[recordNo];
1334 }
1335 {
1336 const uInt *col = (const uInt *)colValues[i++];
1337 ifNo = col[recordNo];
1338 }
1339 {
1340 const Double *col = (const Double *)colValues[i++];
1341 time = col[recordNo];
1342 }
1343 {
1344 const uInt *col = (const uInt *)colValues[i++];
1345 polNo = col[recordNo];
1346 }
1347 assert(nCols == i);
1348 }
1349
1350 if (isFirst) {
1351 enterBeamNo(recordNo, beamNo);
1352 enterIfNo(recordNo, ifNo);
1353 enterTime(recordNo, time);
1354 enterPolNo(recordNo, polNo);
1355 } else {
1356 if (lastBeamNo != beamNo) {
1357 leavePolNo(lastRecordNo, lastPolNo);
1358 leaveTime(lastRecordNo, lastTime);
1359 leaveIfNo(lastRecordNo, lastIfNo);
1360 leaveBeamNo(lastRecordNo, lastBeamNo);
1361
1362 enterBeamNo(recordNo, beamNo);
1363 enterIfNo(recordNo, ifNo);
1364 enterTime(recordNo, time);
1365 enterPolNo(recordNo, polNo);
1366 } else if (lastIfNo != ifNo) {
1367 leavePolNo(lastRecordNo, lastPolNo);
1368 leaveTime(lastRecordNo, lastTime);
1369 leaveIfNo(lastRecordNo, lastIfNo);
1370
1371 enterIfNo(recordNo, ifNo);
1372 enterTime(recordNo, time);
1373 enterPolNo(recordNo, polNo);
1374 } else if (lastTime != time) {
1375 leavePolNo(lastRecordNo, lastPolNo);
1376 leaveTime(lastRecordNo, lastTime);
1377
1378 enterTime(recordNo, time);
1379 enterPolNo(recordNo, polNo);
1380 } else if (lastPolNo != polNo) {
1381 leavePolNo(lastRecordNo, lastPolNo);
1382 enterPolNo(recordNo, polNo);
1383 }
1384 }
1385 count++;
1386 Bool result = visitRecord(recordNo, beamNo, ifNo, polNo, time);
1387
1388 { // epilogue
1389 lastRecordNo = recordNo;
1390
1391 lastBeamNo = beamNo;
1392 lastIfNo = ifNo;
1393 lastPolNo = polNo;
1394 lastTime = time;
1395 }
1396 return result ;
1397 }
1398
1399 virtual void finish() {
1400 if (count > 0) {
1401 leavePolNo(lastRecordNo, lastPolNo);
1402 leaveTime(lastRecordNo, lastTime);
1403 leaveIfNo(lastRecordNo, lastIfNo);
1404 leaveBeamNo(lastRecordNo, lastBeamNo);
1405 }
1406 }
1407};
1408
1409class BaseTsysHolder
1410{
1411public:
1412 BaseTsysHolder( ROArrayColumn<Float> &tsysCol )
1413 : col( tsysCol ),
1414 nchan(0)
1415 {
1416 reset() ;
1417 }
1418 virtual Array<Float> getTsys() = 0 ;
1419 void setNchan( uInt n ) { nchan = n ; }
1420 void appendTsys( uInt row )
1421 {
1422 Vector<Float> v = col( row ) ;
1423 uInt len = tsys.nrow() ;
1424 tsys.resize( len+1, nchan, True ) ;
1425 if ( v.nelements() == nchan )
1426 tsys.row( len ) = v ;
1427 else
1428 tsys.row( len ) = v[0] ;
1429 }
1430 void setTsys( uInt row, uInt idx )
1431 {
1432 if ( idx >= nrow() )
1433 appendTsys( row ) ;
1434 else {
1435 Vector<Float> v = col( row ) ;
1436 if ( v.nelements() == nchan )
1437 tsys.row( idx ) = v ;
1438 else
1439 tsys.row( idx ) = v[0] ;
1440 }
1441 }
1442 void reset()
1443 {
1444 tsys.resize() ;
1445 }
1446 uInt nrow() { return tsys.nrow() ; }
1447 Bool isEffective()
1448 {
1449 return ( !(tsys.empty()) && anyNE( tsys, (Float)1.0 ) ) ;
1450 }
1451 BaseTsysHolder &operator= ( const BaseTsysHolder &v )
1452 {
1453 if ( this != &v )
1454 tsys.assign( v.tsys ) ;
1455 return *this ;
1456 }
1457protected:
1458 ROArrayColumn<Float> col ;
1459 Matrix<Float> tsys ;
1460 uInt nchan ;
1461};
1462
1463class TsysHolder : public BaseTsysHolder
1464{
1465public:
1466 TsysHolder( ROArrayColumn<Float> &tsysCol )
1467 : BaseTsysHolder( tsysCol )
1468 {}
1469 virtual Array<Float> getTsys()
1470 {
1471 return tsys.column( 0 ) ;
1472 }
1473};
1474
1475class TsysSpectrumHolder : public BaseTsysHolder
1476{
1477public:
1478 TsysSpectrumHolder( ROArrayColumn<Float> &tsysCol )
1479 : BaseTsysHolder( tsysCol )
1480 {}
1481 virtual Array<Float> getTsys()
1482 {
1483 return tsys ;
1484 }
1485};
1486
1487class BaseTcalProcessor
1488{
1489public:
1490 BaseTcalProcessor( ROArrayColumn<Float> &tcalCol )
1491 : col( tcalCol )
1492 {}
1493 void setTcalId( Vector<uInt> &tcalId ) { id.assign( tcalId ) ; }
1494 virtual Array<Float> getTcal() = 0 ;
1495protected:
1496 ROArrayColumn<Float> col ;
1497 Vector<uInt> id ;
1498};
1499
1500class TcalProcessor : public BaseTcalProcessor
1501{
1502public:
1503 TcalProcessor( ROArrayColumn<Float> &tcalCol )
1504 : BaseTcalProcessor( tcalCol )
1505 {}
1506 virtual Array<Float> getTcal()
1507 {
1508 uInt npol = id.nelements() ;
1509 Vector<Float> tcal( npol ) ;
1510 for ( uInt ipol = 0 ; ipol < npol ; ipol++ )
1511 tcal[ipol] = col( id[ipol] ).data()[0] ;
1512 //cout << "TcalProcessor: tcal = " << tcal << endl ;
1513 return tcal ;
1514 }
1515};
1516
1517class TcalSpectrumProcessor : public BaseTcalProcessor
1518{
1519public:
1520 TcalSpectrumProcessor( ROArrayColumn<Float> &tcalCol )
1521 : BaseTcalProcessor( tcalCol )
1522 {}
1523 virtual Array<Float> getTcal()
1524 {
1525 uInt npol = id.nelements() ;
1526 Vector<Float> tcal0 = col( 0 ) ;
1527 uInt nchan = tcal0.nelements() ;
1528 Matrix<Float> tcal( npol, nchan ) ;
1529 tcal.row( 0 ) = tcal0 ;
1530 for ( uInt ipol = 1 ; ipol < npol ; ipol++ )
1531 tcal.row( ipol ) = col( id[ipol] ) ;
1532 return tcal ;
1533 }
1534};
1535
1536class MSSysCalVisitor : public BaseMSSysCalVisitor
1537{
1538public:
1539 MSSysCalVisitor( const Table &from, Table &to )
1540 : BaseMSSysCalVisitor( from ),
1541 sctab( to ),
1542 rowidx( 0 )
1543 {
1544 scrow = TableRow( sctab ) ;
1545
1546 lastTcalId.resize() ;
1547 theTcalId.resize() ;
1548 startTime = 0.0 ;
1549 endTime = 0.0 ;
1550
1551 const TableRecord &keys = table.keywordSet() ;
1552 Table tcalTable = keys.asTable( "TCAL" ) ;
1553 tcalCol.attach( tcalTable, "TCAL" ) ;
1554 tsysCol.attach( table, "TSYS" ) ;
1555 tcalIdCol.attach( table, "TCAL_ID" ) ;
1556 intervalCol.attach( table, "INTERVAL" ) ;
1557 effectiveTcal.resize( tcalTable.nrow() ) ;
1558 for ( uInt irow = 0 ; irow < tcalTable.nrow() ; irow++ ) {
1559 if ( allEQ( tcalCol( irow ), (Float)1.0 ) )
1560 effectiveTcal[irow] = False ;
1561 else
1562 effectiveTcal[irow] = True ;
1563 }
1564
1565 TableRecord &r = scrow.record() ;
1566 RecordFieldPtr<Int> antennaIdRF( r, "ANTENNA_ID" ) ;
1567 *antennaIdRF = 0 ;
1568 feedIdRF.attachToRecord( r, "FEED_ID" ) ;
1569 specWinIdRF.attachToRecord( r, "SPECTRAL_WINDOW_ID" ) ;
1570 timeRF.attachToRecord( r, "TIME" ) ;
1571 intervalRF.attachToRecord( r, "INTERVAL" ) ;
1572 if ( r.isDefined( "TCAL" ) ) {
1573 tcalRF.attachToRecord( r, "TCAL" ) ;
1574 tcalProcessor = new TcalProcessor( tcalCol ) ;
1575 }
1576 else if ( r.isDefined( "TCAL_SPECTRUM" ) ) {
1577 tcalRF.attachToRecord( r, "TCAL_SPECTRUM" ) ;
1578 tcalProcessor = new TcalSpectrumProcessor( tcalCol ) ;
1579 }
1580 if ( r.isDefined( "TSYS" ) ) {
1581 tsysRF.attachToRecord( r, "TSYS" ) ;
1582 theTsys = new TsysHolder( tsysCol ) ;
1583 lastTsys = new TsysHolder( tsysCol ) ;
1584 }
1585 else {
1586 tsysRF.attachToRecord( r, "TSYS_SPECTRUM" ) ;
1587 theTsys = new TsysSpectrumHolder( tsysCol ) ;
1588 lastTsys = new TsysSpectrumHolder( tsysCol ) ;
1589 }
1590
1591 }
1592
1593 virtual void enterBeamNo(const uInt recordNo, uInt columnValue)
1594 {
1595 *feedIdRF = (Int)columnValue ;
1596 }
1597 virtual void leaveBeamNo(const uInt recordNo, uInt columnValue)
1598 {
1599 }
1600 virtual void enterIfNo(const uInt recordNo, uInt columnValue)
1601 {
1602 //cout << "enterIfNo" << endl ;
1603 ROArrayColumn<Float> sp( table, "SPECTRA" ) ;
1604 uInt nchan = sp( recordNo ).nelements() ;
1605 theTsys->setNchan( nchan ) ;
1606 lastTsys->setNchan( nchan ) ;
1607
1608 *specWinIdRF = (Int)columnValue ;
1609 }
1610 virtual void leaveIfNo(const uInt recordNo, uInt columnValue)
1611 {
1612 //cout << "leaveIfNo" << endl ;
1613 post() ;
1614 lastTsys->reset() ;
1615 lastTcalId.resize() ;
1616 theTsys->reset() ;
1617 theTcalId.resize() ;
1618 startTime = 0.0 ;
1619 endTime = 0.0 ;
1620 }
1621 virtual void enterTime(const uInt recordNo, Double columnValue)
1622 {
1623 //cout << "enterTime" << endl ;
1624 interval = intervalCol.asdouble( recordNo ) ;
1625 // start time and end time
1626 if ( startTime == 0.0 ) {
1627 startTime = columnValue * 86400.0 - 0.5 * interval ;
1628 endTime = columnValue * 86400.0 + 0.5 * interval ;
1629 }
1630 }
1631 virtual void leaveTime(const uInt recordNo, Double columnValue)
1632 {
1633 //cout << "leaveTime" << endl ;
1634 if ( isUpdated() ) {
1635 post() ;
1636 *lastTsys = *theTsys ;
1637 lastTcalId = theTcalId ;
1638 theTsys->reset() ;
1639 theTcalId.resize() ;
1640 startTime = columnValue * 86400.0 - 0.5 * interval ;
1641 endTime = columnValue * 86400.0 + 0.5 * interval ;
1642 }
1643 else {
1644 endTime = columnValue * 86400.0 + 0.5 * interval ;
1645 }
1646 }
1647 virtual void enterPolNo(const uInt recordNo, uInt columnValue)
1648 {
1649 //cout << "enterPolNo" << endl ;
1650 Vector<Float> tsys = tsysCol( recordNo ) ;
1651 uInt tcalId = tcalIdCol.asuInt( recordNo ) ;
1652 // lastTsys.nrow() must be npol
1653 if ( lastTsys->nrow() == columnValue )
1654 lastTsys->appendTsys( recordNo ) ;
1655 // lastTcalId.nelements() must be npol
1656 if ( lastTcalId.nelements() == columnValue )
1657 appendTcalId( lastTcalId, tcalId, columnValue ) ;
1658 // theTsys.nrow() must be npol
1659 if ( theTsys->nrow() == columnValue )
1660 theTsys->appendTsys( recordNo ) ;
1661 else {
1662 theTsys->setTsys( recordNo, columnValue ) ;
1663 }
1664 if ( theTcalId.nelements() == columnValue )
1665 appendTcalId( theTcalId, tcalId, columnValue ) ;
1666 else
1667 setTcalId( theTcalId, tcalId, columnValue ) ;
1668 }
1669 virtual void leavePolNo( const uInt recordNo, uInt columnValue )
1670 {
1671 }
1672
1673private:
1674 void appendTcalId( Vector<uInt> &v, uInt &elem, uInt &polId )
1675 {
1676 v.resize( polId+1, True ) ;
1677 v[polId] = elem ;
1678 }
1679 void setTcalId( Vector<uInt> &v, uInt &elem, uInt &polId )
1680 {
1681 v[polId] = elem ;
1682 }
1683 void post()
1684 {
1685 // check if given Tcal and Tsys is effective
1686 Bool isEffective = False ;
1687 for ( uInt ipol = 0 ; ipol < lastTcalId.nelements() ; ipol++ ) {
1688 if ( effectiveTcal[lastTcalId[ipol]] ) {
1689 isEffective = True ;
1690 break ;
1691 }
1692 }
1693 if ( !isEffective ) {
1694 if ( !(lastTsys->isEffective()) )
1695 return ;
1696 }
1697
1698 //cout << " interval: " << (endTime-startTime) << " lastTcalId = " << lastTcalId << endl ;
1699 Double midTime = 0.5 * ( startTime + endTime ) ;
1700 Double interval = endTime - startTime ;
1701 *timeRF = midTime ;
1702 *intervalRF = interval ;
1703 tcalProcessor->setTcalId( lastTcalId ) ;
1704 Array<Float> tcal = tcalProcessor->getTcal() ;
1705 tcalRF.define( tcal ) ;
1706 tsysRF.define( lastTsys->getTsys() ) ;
1707 sctab.addRow( 1, True ) ;
1708 scrow.put( rowidx ) ;
1709 rowidx++ ;
1710 }
1711
1712 Bool isUpdated()
1713 {
1714 Bool ret = False ;
1715 ret = anyNE( theTcalId, lastTcalId ) ;
1716 if ( !ret )
1717 ret = anyNE( theTsys->getTsys(), lastTsys->getTsys() ) ;
1718 return ret ;
1719 }
1720
1721 Table &sctab;
1722 TableRow scrow;
1723 uInt rowidx;
1724
1725 Double startTime,endTime,interval;
1726
1727 CountedPtr<BaseTsysHolder> lastTsys,theTsys;
1728 Vector<uInt> lastTcalId,theTcalId;
1729 CountedPtr<BaseTcalProcessor> tcalProcessor ;
1730 Vector<Bool> effectiveTcal;
1731
1732 RecordFieldPtr<Int> feedIdRF,specWinIdRF;
1733 RecordFieldPtr<Double> timeRF,intervalRF;
1734 RecordFieldPtr< Array<Float> > tcalRF,tsysRF;
1735
1736 ROArrayColumn<Float> tsysCol,tcalCol;
1737 ROTableColumn tcalIdCol,intervalCol;
1738};
1739
1740MSWriter::MSWriter(CountedPtr<Scantable> stable)
1741 : table_(stable),
1742 isWeather_(False),
1743 tcalSpec_(False),
1744 tsysSpec_(False),
1745 ptTabName_("")
1746{
1747 os_ = LogIO() ;
1748 os_.origin( LogOrigin( "MSWriter", "MSWriter()", WHERE ) ) ;
1749// os_ << "MSWriter::MSWriter()" << LogIO::POST ;
1750
1751 // initialize writer
1752 init() ;
1753}
1754
1755MSWriter::~MSWriter()
1756{
1757 os_.origin( LogOrigin( "MSWriter", "~MSWriter()", WHERE ) ) ;
1758// os_ << "MSWriter::~MSWriter()" << LogIO::POST ;
1759
1760 if ( mstable_ != 0 )
1761 delete mstable_ ;
1762}
1763
1764bool MSWriter::write(const string& filename, const Record& rec)
1765{
1766 os_.origin( LogOrigin( "MSWriter", "write()", WHERE ) ) ;
1767 //double startSec = mathutil::gettimeofday_sec() ;
1768 //os_ << "start MSWriter::write() startSec=" << startSec << LogIO::POST ;
1769
1770 filename_ = filename ;
1771
1772 // parsing MS options
1773 Bool overwrite = False ;
1774 if ( rec.isDefined( "ms" ) ) {
1775 Record msrec = rec.asRecord( "ms" ) ;
1776 if ( msrec.isDefined( "overwrite" ) ) {
1777 overwrite = msrec.asBool( "overwrite" ) ;
1778 }
1779 }
1780
1781 os_ << "Parsing MS options" << endl ;
1782 os_ << " overwrite = " << overwrite << LogIO::POST ;
1783
1784 File file( filename_ ) ;
1785 if ( file.exists() ) {
1786 if ( overwrite ) {
1787 os_ << filename_ << " exists. Overwrite existing data... " << LogIO::POST ;
1788 if ( file.isRegular() ) RegularFile(file).remove() ;
1789 else if ( file.isDirectory() ) Directory(file).removeRecursive() ;
1790 else SymLink(file).remove() ;
1791 }
1792 else {
1793 os_ << LogIO::SEVERE << "ERROR: " << filename_ << " exists..." << LogIO::POST ;
1794 return False ;
1795 }
1796 }
1797
1798 // set up MS
1799 setupMS() ;
1800
1801 // subtables
1802 // OBSERVATION
1803 fillObservation() ;
1804
1805 // ANTENNA
1806 fillAntenna() ;
1807
1808 // PROCESSOR
1809 fillProcessor() ;
1810
1811 // SOURCE
1812 fillSource() ;
1813
1814 // WEATHER
1815 if ( isWeather_ )
1816 fillWeather() ;
1817
1818 // SYSCAL
1819 fillSysCal() ;
1820
1821 /***
1822 * Start iteration using TableVisitor
1823 ***/
1824 {
1825 static const char *cols[] = {
1826 "FIELDNAME", "BEAMNO", "SCANNO", "IFNO", "SRCTYPE", "CYCLENO", "TIME",
1827 "POLNO",
1828 NULL
1829 };
1830 static const TypeManagerImpl<uInt> tmUInt;
1831 static const TypeManagerImpl<Int> tmInt;
1832 static const TypeManagerImpl<Double> tmDouble;
1833 static const TypeManagerImpl<String> tmString;
1834 static const TypeManager *const tms[] = {
1835 &tmString, &tmUInt, &tmUInt, &tmUInt, &tmInt, &tmUInt, &tmDouble, &tmUInt, NULL
1836 };
1837 //double t0 = mathutil::gettimeofday_sec() ;
1838 MSWriterVisitor myVisitor(table_->table(),*mstable_);
1839 //double t1 = mathutil::gettimeofday_sec() ;
1840 //cout << "MSWriterVisitor(): elapsed time " << t1-t0 << " sec" << endl ;
1841 String dataColName = "FLOAT_DATA" ;
1842 if ( useData_ )
1843 dataColName = "DATA" ;
1844 myVisitor.dataColumnName( dataColName ) ;
1845 myVisitor.pointingTableName( ptTabName_ ) ;
1846 myVisitor.setSourceRecord( srcRec_ ) ;
1847 //double t2 = mathutil::gettimeofday_sec() ;
1848 traverseTable(table_->table(), cols, tms, &myVisitor);
1849 //double t3 = mathutil::gettimeofday_sec() ;
1850 //cout << "traverseTable(): elapsed time " << t3-t2 << " sec" << endl ;
1851 }
1852 /***
1853 * End iteration using TableVisitor
1854 ***/
1855
1856 // ASDM tables
1857 const TableRecord &stKeys = table_->table().keywordSet() ;
1858 TableRecord &msKeys = mstable_->rwKeywordSet() ;
1859 uInt nfields = stKeys.nfields() ;
1860 for ( uInt ifield = 0 ; ifield < nfields ; ifield++ ) {
1861 String kname = stKeys.name( ifield ) ;
1862 if ( kname.find( "ASDM" ) != String::npos ) {
1863 String asdmpath = stKeys.asString( ifield ) ;
1864 os_ << "found ASDM table: " << asdmpath << LogIO::POST ;
1865 if ( Table::isReadable( asdmpath ) ) {
1866 Table newAsdmTab( asdmpath, Table::Old ) ;
1867 newAsdmTab.copy( filename_+"/"+kname, Table::New ) ;
1868 os_ << "add subtable: " << kname << LogIO::POST ;
1869 msKeys.defineTable( kname, Table( filename_+"/"+kname, Table::Old ) ) ;
1870 }
1871 }
1872 }
1873
1874 // replace POINTING table with original one if exists
1875 if ( ptTabName_ != "" ) {
1876 delete mstable_ ;
1877 mstable_ = 0 ;
1878 Table newPtTab( ptTabName_, Table::Old ) ;
1879 newPtTab.copy( filename_+"/POINTING", Table::New ) ;
1880 }
1881
1882 //double endSec = mathutil::gettimeofday_sec() ;
1883 //os_ << "end MSWriter::write() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1884
1885 return True ;
1886}
1887
1888void MSWriter::init()
1889{
1890// os_.origin( LogOrigin( "MSWriter", "init()", WHERE ) ) ;
1891// double startSec = mathutil::gettimeofday_sec() ;
1892// os_ << "start MSWriter::init() startSec=" << startSec << LogIO::POST ;
1893
1894 // access to scantable
1895 header_ = table_->getHeader() ;
1896
1897 // FLOAT_DATA? or DATA?
1898 if ( header_.npol > 2 ) {
1899 useFloatData_ = False ;
1900 useData_ = True ;
1901 }
1902 else {
1903 useFloatData_ = True ;
1904 useData_ = False ;
1905 }
1906
1907 // polarization type
1908 polType_ = header_.poltype ;
1909 if ( polType_ == "" )
1910 polType_ = "stokes" ;
1911 else if ( polType_.find( "linear" ) != String::npos )
1912 polType_ = "linear" ;
1913 else if ( polType_.find( "circular" ) != String::npos )
1914 polType_ = "circular" ;
1915 else if ( polType_.find( "stokes" ) != String::npos )
1916 polType_ = "stokes" ;
1917 else if ( polType_.find( "linpol" ) != String::npos )
1918 polType_ = "linpol" ;
1919 else
1920 polType_ = "notype" ;
1921
1922 // Check if some subtables are exists
1923 Bool isTcal = False ;
1924 if ( table_->tcal().table().nrow() != 0 ) {
1925 ROTableColumn col( table_->tcal().table(), "TCAL" ) ;
1926 if ( col.isDefined( 0 ) ) {
1927 os_ << "TCAL table exists: nrow=" << table_->tcal().table().nrow() << LogIO::POST ;
1928 isTcal = True ;
1929 }
1930 else {
1931 os_ << "No TCAL rows" << LogIO::POST ;
1932 }
1933 }
1934 else {
1935 os_ << "No TCAL rows" << LogIO::POST ;
1936 }
1937 if ( table_->weather().table().nrow() != 0 ) {
1938 ROTableColumn col( table_->weather().table(), "TEMPERATURE" ) ;
1939 if ( col.isDefined( 0 ) ) {
1940 os_ << "WEATHER table exists: nrow=" << table_->weather().table().nrow() << LogIO::POST ;
1941 isWeather_ =True ;
1942 }
1943 else {
1944 os_ << "No WEATHER rows" << LogIO::POST ;
1945 }
1946 }
1947 else {
1948 os_ << "No WEATHER rows" << LogIO::POST ;
1949 }
1950
1951 // Are TCAL_SPECTRUM and TSYS_SPECTRUM necessary?
1952 if ( header_.nchan != 1 ) {
1953 if ( isTcal ) {
1954 // examine TCAL subtable
1955 Table tcaltab = table_->tcal().table() ;
1956 ROArrayColumn<Float> tcalCol( tcaltab, "TCAL" ) ;
1957 for ( uInt irow = 0 ; irow < tcaltab.nrow() ; irow++ ) {
1958 if ( tcalCol( irow ).size() != 1 )
1959 tcalSpec_ = True ;
1960 }
1961 }
1962 // examine spectral data
1963 TableIterator iter0( table_->table(), "IFNO" ) ;
1964 while( !iter0.pastEnd() ) {
1965 Table t0( iter0.table() ) ;
1966 ROArrayColumn<Float> sharedFloatArrCol( t0, "SPECTRA" ) ;
1967 uInt len = sharedFloatArrCol( 0 ).size() ;
1968 if ( len != 1 ) {
1969 sharedFloatArrCol.attach( t0, "TSYS" ) ;
1970 if ( sharedFloatArrCol( 0 ).size() != 1 )
1971 tsysSpec_ = True ;
1972 }
1973 iter0.next() ;
1974 }
1975 }
1976
1977 // check if reference for POINTING table exists
1978 const TableRecord &rec = table_->table().keywordSet() ;
1979 if ( rec.isDefined( "POINTING" ) ) {
1980 ptTabName_ = rec.asString( "POINTING" ) ;
1981 if ( !Table::isReadable( ptTabName_ ) ) {
1982 ptTabName_ = "" ;
1983 }
1984 }
1985
1986// double endSec = mathutil::gettimeofday_sec() ;
1987// os_ << "end MSWriter::init() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1988}
1989
1990void MSWriter::setupMS()
1991{
1992// os_.origin( LogOrigin( "MSWriter", "setupMS()", WHERE ) ) ;
1993// double startSec = mathutil::gettimeofday_sec() ;
1994// os_ << "start MSWriter::setupMS() startSec=" << startSec << LogIO::POST ;
1995
1996 String dunit = table_->getHeader().fluxunit ;
1997
1998 TableDesc msDesc = MeasurementSet::requiredTableDesc() ;
1999 if ( useFloatData_ )
2000 MeasurementSet::addColumnToDesc( msDesc, MSMainEnums::FLOAT_DATA, 2 ) ;
2001 else if ( useData_ )
2002 MeasurementSet::addColumnToDesc( msDesc, MSMainEnums::DATA, 2 ) ;
2003
2004 SetupNewTable newtab( filename_, msDesc, Table::New ) ;
2005
2006 mstable_ = new MeasurementSet( newtab ) ;
2007
2008 TableColumn col ;
2009 if ( useFloatData_ )
2010 col.attach( *mstable_, "FLOAT_DATA" ) ;
2011 else if ( useData_ )
2012 col.attach( *mstable_, "DATA" ) ;
2013 col.rwKeywordSet().define( "UNIT", dunit ) ;
2014
2015 // create subtables
2016 TableDesc antennaDesc = MSAntenna::requiredTableDesc() ;
2017 SetupNewTable antennaTab( mstable_->antennaTableName(), antennaDesc, Table::New ) ;
2018 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::ANTENNA ), Table( antennaTab ) ) ;
2019
2020 TableDesc dataDescDesc = MSDataDescription::requiredTableDesc() ;
2021 SetupNewTable dataDescTab( mstable_->dataDescriptionTableName(), dataDescDesc, Table::New ) ;
2022 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::DATA_DESCRIPTION ), Table( dataDescTab ) ) ;
2023
2024 TableDesc dopplerDesc = MSDoppler::requiredTableDesc() ;
2025 SetupNewTable dopplerTab( mstable_->dopplerTableName(), dopplerDesc, Table::New ) ;
2026 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::DOPPLER ), Table( dopplerTab ) ) ;
2027
2028 TableDesc feedDesc = MSFeed::requiredTableDesc() ;
2029 SetupNewTable feedTab( mstable_->feedTableName(), feedDesc, Table::New ) ;
2030 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FEED ), Table( feedTab ) ) ;
2031
2032 TableDesc fieldDesc = MSField::requiredTableDesc() ;
2033 SetupNewTable fieldTab( mstable_->fieldTableName(), fieldDesc, Table::New ) ;
2034 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FIELD ), Table( fieldTab ) ) ;
2035
2036 TableDesc flagCmdDesc = MSFlagCmd::requiredTableDesc() ;
2037 SetupNewTable flagCmdTab( mstable_->flagCmdTableName(), flagCmdDesc, Table::New ) ;
2038 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FLAG_CMD ), Table( flagCmdTab ) ) ;
2039
2040 TableDesc freqOffsetDesc = MSFreqOffset::requiredTableDesc() ;
2041 SetupNewTable freqOffsetTab( mstable_->freqOffsetTableName(), freqOffsetDesc, Table::New ) ;
2042 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FREQ_OFFSET ), Table( freqOffsetTab ) ) ;
2043
2044 TableDesc historyDesc = MSHistory::requiredTableDesc() ;
2045 SetupNewTable historyTab( mstable_->historyTableName(), historyDesc, Table::New ) ;
2046 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::HISTORY ), Table( historyTab ) ) ;
2047
2048 TableDesc observationDesc = MSObservation::requiredTableDesc() ;
2049 SetupNewTable observationTab( mstable_->observationTableName(), observationDesc, Table::New ) ;
2050 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::OBSERVATION ), Table( observationTab ) ) ;
2051
2052 TableDesc pointingDesc = MSPointing::requiredTableDesc() ;
2053 SetupNewTable pointingTab( mstable_->pointingTableName(), pointingDesc, Table::New ) ;
2054 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::POINTING ), Table( pointingTab ) ) ;
2055
2056 TableDesc polarizationDesc = MSPolarization::requiredTableDesc() ;
2057 SetupNewTable polarizationTab( mstable_->polarizationTableName(), polarizationDesc, Table::New ) ;
2058 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::POLARIZATION ), Table( polarizationTab ) ) ;
2059
2060 TableDesc processorDesc = MSProcessor::requiredTableDesc() ;
2061 SetupNewTable processorTab( mstable_->processorTableName(), processorDesc, Table::New ) ;
2062 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::PROCESSOR ), Table( processorTab ) ) ;
2063
2064 TableDesc sourceDesc = MSSource::requiredTableDesc() ;
2065 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::TRANSITION, 1 ) ;
2066 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::REST_FREQUENCY, 1 ) ;
2067 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::SYSVEL, 1 ) ;
2068 SetupNewTable sourceTab( mstable_->sourceTableName(), sourceDesc, Table::New ) ;
2069 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SOURCE ), Table( sourceTab ) ) ;
2070
2071 TableDesc spwDesc = MSSpectralWindow::requiredTableDesc() ;
2072 SetupNewTable spwTab( mstable_->spectralWindowTableName(), spwDesc, Table::New ) ;
2073 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SPECTRAL_WINDOW ), Table( spwTab ) ) ;
2074
2075 TableDesc stateDesc = MSState::requiredTableDesc() ;
2076 SetupNewTable stateTab( mstable_->stateTableName(), stateDesc, Table::New ) ;
2077 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::STATE ), Table( stateTab ) ) ;
2078
2079 TableDesc sysCalDesc = MSSysCal::requiredTableDesc() ;
2080 if ( tcalSpec_ )
2081 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TCAL_SPECTRUM, 2 ) ;
2082 else
2083 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TCAL, 1 ) ;
2084 if ( tsysSpec_ )
2085 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TSYS_SPECTRUM, 2 ) ;
2086 else
2087 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TSYS, 1 ) ;
2088 SetupNewTable sysCalTab( mstable_->sysCalTableName(), sysCalDesc, Table::New ) ;
2089 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SYSCAL ), Table( sysCalTab ) ) ;
2090
2091 TableDesc weatherDesc = MSWeather::requiredTableDesc() ;
2092 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::TEMPERATURE ) ;
2093 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::PRESSURE ) ;
2094 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::REL_HUMIDITY ) ;
2095 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::WIND_SPEED ) ;
2096 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::WIND_DIRECTION ) ;
2097 SetupNewTable weatherTab( mstable_->weatherTableName(), weatherDesc, Table::New ) ;
2098 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::WEATHER ), Table( weatherTab ) ) ;
2099
2100 mstable_->initRefs() ;
2101
2102// double endSec = mathutil::gettimeofday_sec() ;
2103// os_ << "end MSWriter::setupMS() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2104}
2105
2106void MSWriter::fillObservation()
2107{
2108 //double startSec = mathutil::gettimeofday_sec() ;
2109 //os_ << "start MSWriter::fillObservation() startSec=" << startSec << LogIO::POST ;
2110
2111 // only 1 row
2112 mstable_->observation().addRow( 1, True ) ;
2113 MSObservationColumns msObsCols( mstable_->observation() ) ;
2114 msObsCols.observer().put( 0, header_.observer ) ;
2115 // tentatively put antennaname (from ANTENNA subtable)
2116 String hAntennaName = header_.antennaname ;
2117 String::size_type pos = hAntennaName.find( "//" ) ;
2118 String telescopeName ;
2119 if ( pos != String::npos ) {
2120 telescopeName = hAntennaName.substr( 0, pos ) ;
2121 }
2122 else {
2123 pos = hAntennaName.find( "@" ) ;
2124 telescopeName = hAntennaName.substr( 0, pos ) ;
2125 }
2126// os_ << "telescopeName = " << telescopeName << LogIO::POST ;
2127 msObsCols.telescopeName().put( 0, telescopeName ) ;
2128 msObsCols.project().put( 0, header_.project ) ;
2129 //ScalarMeasColumn<MEpoch> timeCol( table_->table().sort("TIME"), "TIME" ) ;
2130 Table sortedtable = table_->table().sort("TIME") ;
2131 ScalarMeasColumn<MEpoch> timeCol( sortedtable, "TIME" ) ;
2132 Vector<MEpoch> trange( 2 ) ;
2133 trange[0] = timeCol( 0 ) ;
2134 trange[1] = timeCol( table_->nrow()-1 ) ;
2135 msObsCols.timeRangeMeas().put( 0, trange ) ;
2136
2137 //double endSec = mathutil::gettimeofday_sec() ;
2138 //os_ << "end MSWriter::fillObservation() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2139}
2140
2141void MSWriter::antennaProperty( String &name, String &m, String &t, Double &d )
2142{
2143 name.upcase() ;
2144
2145 m = "ALT-AZ" ;
2146 t = "GROUND-BASED" ;
2147 if ( name.matches( Regex( "DV[0-9]+$" ) )
2148 || name.matches( Regex( "DA[0-9]+$" ) )
2149 || name.matches( Regex( "PM[0-9]+$" ) ) )
2150 d = 12.0 ;
2151 else if ( name.matches( Regex( "CM[0-9]+$" ) ) )
2152 d = 7.0 ;
2153 else if ( name.contains( "GBT" ) )
2154 d = 104.9 ;
2155 else if ( name.contains( "MOPRA" ) )
2156 d = 22.0 ;
2157 else if ( name.contains( "PKS" ) || name.contains( "PARKS" ) )
2158 d = 64.0 ;
2159 else if ( name.contains( "TIDBINBILLA" ) )
2160 d = 70.0 ;
2161 else if ( name.contains( "CEDUNA" ) )
2162 d = 30.0 ;
2163 else if ( name.contains( "HOBART" ) )
2164 d = 26.0 ;
2165 else if ( name.contains( "APEX" ) )
2166 d = 12.0 ;
2167 else if ( name.contains( "ASTE" ) )
2168 d = 10.0 ;
2169 else if ( name.contains( "NRO" ) )
2170 d = 45.0 ;
2171 else
2172 d = 1.0 ;
2173}
2174
2175void MSWriter::fillAntenna()
2176{
2177 //double startSec = mathutil::gettimeofday_sec() ;
2178 //os_ << "start MSWriter::fillAntenna() startSec=" << startSec << LogIO::POST ;
2179
2180 // only 1 row
2181 Table anttab = mstable_->antenna() ;
2182 anttab.addRow( 1, True ) ;
2183
2184 Table &table = table_->table() ;
2185 const TableRecord &keys = table.keywordSet() ;
2186 String hAntName = keys.asString( "AntennaName" ) ;
2187 String::size_type pos = hAntName.find( "//" ) ;
2188 String antennaName ;
2189 String stationName ;
2190 if ( pos != String::npos ) {
2191 stationName = hAntName.substr( 0, pos ) ;
2192 hAntName = hAntName.substr( pos+2 ) ;
2193 }
2194 pos = hAntName.find( "@" ) ;
2195 if ( pos != String::npos ) {
2196 antennaName = hAntName.substr( 0, pos ) ;
2197 stationName = hAntName.substr( pos+1 ) ;
2198 }
2199 else {
2200 antennaName = hAntName ;
2201 }
2202 Vector<Double> antpos = keys.asArrayDouble( "AntennaPosition" ) ;
2203
2204 String mount, atype ;
2205 Double diameter ;
2206 antennaProperty( antennaName, mount, atype, diameter ) ;
2207
2208 TableRow tr( anttab ) ;
2209 TableRecord &r = tr.record() ;
2210 RecordFieldPtr<String> nameRF( r, "NAME" ) ;
2211 RecordFieldPtr<String> stationRF( r, "STATION" ) ;
2212 RecordFieldPtr<String> mountRF( r, "NAME" ) ;
2213 RecordFieldPtr<String> typeRF( r, "TYPE" ) ;
2214 RecordFieldPtr<Double> dishDiameterRF( r, "DISH_DIAMETER" ) ;
2215 RecordFieldPtr< Vector<Double> > positionRF( r, "POSITION" ) ;
2216 *nameRF = antennaName ;
2217 *mountRF = mount ;
2218 *typeRF = atype ;
2219 *dishDiameterRF = diameter ;
2220 *positionRF = antpos ;
2221
2222 tr.put( 0 ) ;
2223
2224 //double endSec = mathutil::gettimeofday_sec() ;
2225 //os_ << "end MSWriter::fillAntenna() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2226}
2227
2228void MSWriter::fillProcessor()
2229{
2230// double startSec = mathutil::gettimeofday_sec() ;
2231// os_ << "start MSWriter::fillProcessor() startSec=" << startSec << LogIO::POST ;
2232
2233 // only add empty 1 row
2234 MSProcessor msProc = mstable_->processor() ;
2235 msProc.addRow( 1, True ) ;
2236
2237// double endSec = mathutil::gettimeofday_sec() ;
2238// os_ << "end MSWriter::fillProcessor() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2239}
2240
2241void MSWriter::fillSource()
2242{
2243// double startSec = mathutil::gettimeofday_sec() ;
2244// os_ << "start MSWriter::fillSource() startSec=" << startSec << LogIO::POST ;
2245
2246 // access to MS SOURCE subtable
2247 MSSource msSrc = mstable_->source() ;
2248
2249 // access to MOLECULE subtable
2250 STMolecules stm = table_->molecules() ;
2251
2252 Int srcId = 0 ;
2253 Vector<Double> restFreq ;
2254 Vector<String> molName ;
2255 Vector<String> fMolName ;
2256
2257 // row based
2258 TableRow row( msSrc ) ;
2259 TableRecord &rec = row.record() ;
2260 RecordFieldPtr<Int> srcidRF( rec, "SOURCE_ID" ) ;
2261 RecordFieldPtr<String> nameRF( rec, "NAME" ) ;
2262 RecordFieldPtr< Array<Double> > srcpmRF( rec, "PROPER_MOTION" ) ;
2263 RecordFieldPtr< Array<Double> > srcdirRF( rec, "DIRECTION" ) ;
2264 RecordFieldPtr<Int> numlineRF( rec, "NUM_LINES" ) ;
2265 RecordFieldPtr< Array<Double> > restfreqRF( rec, "REST_FREQUENCY" ) ;
2266 RecordFieldPtr< Array<Double> > sysvelRF( rec, "SYSVEL" ) ;
2267 RecordFieldPtr< Array<String> > transitionRF( rec, "TRANSITION" ) ;
2268 RecordFieldPtr<Double> timeRF( rec, "TIME" ) ;
2269 RecordFieldPtr<Double> intervalRF( rec, "INTERVAL" ) ;
2270 RecordFieldPtr<Int> spwidRF( rec, "SPECTRAL_WINDOW_ID" ) ;
2271
2272 //
2273 // ITERATION: SRCNAME
2274 //
2275 TableIterator iter0( table_->table(), "SRCNAME" ) ;
2276 while( !iter0.pastEnd() ) {
2277 //Table t0( iter0.table() ) ;
2278 Table t0 = iter0.table() ;
2279
2280 // get necessary information
2281 ROScalarColumn<String> srcNameCol( t0, "SRCNAME" ) ;
2282 String srcName = srcNameCol( 0 ) ;
2283 ROArrayColumn<Double> sharedDArrRCol( t0, "SRCPROPERMOTION" ) ;
2284 Vector<Double> srcPM = sharedDArrRCol( 0 ) ;
2285 sharedDArrRCol.attach( t0, "SRCDIRECTION" ) ;
2286 Vector<Double> srcDir = sharedDArrRCol( 0 ) ;
2287 ROScalarColumn<Double> srcVelCol( t0, "SRCVELOCITY" ) ;
2288 Double srcVel = srcVelCol( 0 ) ;
2289 srcRec_.define( srcName, srcId ) ;
2290
2291 // NAME
2292 *nameRF = srcName ;
2293
2294 // SOURCE_ID
2295 *srcidRF = srcId ;
2296
2297 // PROPER_MOTION
2298 *srcpmRF = srcPM ;
2299
2300 // DIRECTION
2301 *srcdirRF = srcDir ;
2302
2303 //
2304 // ITERATION: MOLECULE_ID
2305 //
2306 TableIterator iter1( t0, "MOLECULE_ID" ) ;
2307 while( !iter1.pastEnd() ) {
2308 //Table t1( iter1.table() ) ;
2309 Table t1 = iter1.table() ;
2310
2311 // get necessary information
2312 ROScalarColumn<uInt> molIdCol( t1, "MOLECULE_ID" ) ;
2313 uInt molId = molIdCol( 0 ) ;
2314 stm.getEntry( restFreq, molName, fMolName, molId ) ;
2315
2316 uInt numFreq = restFreq.size() ;
2317
2318 // NUM_LINES
2319 *numlineRF = numFreq ;
2320
2321 // REST_FREQUENCY
2322 *restfreqRF = restFreq ;
2323
2324 // TRANSITION
2325 //*transitionRF = fMolName ;
2326 Vector<String> transition ;
2327 if ( fMolName.size() != 0 ) {
2328 transition = fMolName ;
2329 }
2330 else if ( molName.size() != 0 ) {
2331 transition = molName ;
2332 }
2333 else {
2334 transition.resize( numFreq ) ;
2335 transition = "" ;
2336 }
2337 *transitionRF = transition ;
2338
2339 // SYSVEL
2340 Vector<Double> sysvelArr( numFreq, srcVel ) ;
2341 *sysvelRF = sysvelArr ;
2342
2343 //
2344 // ITERATION: IFNO
2345 //
2346 TableIterator iter2( t1, "IFNO" ) ;
2347 while( !iter2.pastEnd() ) {
2348 //Table t2( iter2.table() ) ;
2349 Table t2 = iter2.table() ;
2350 uInt nrow = msSrc.nrow() ;
2351
2352 // get necessary information
2353 ROScalarColumn<uInt> ifNoCol( t2, "IFNO" ) ;
2354 uInt ifno = ifNoCol( 0 ) ; // IFNO = SPECTRAL_WINDOW_ID
2355 Double midTime ;
2356 Double interval ;
2357 getValidTimeRange( midTime, interval, t2 ) ;
2358
2359 // fill SPECTRAL_WINDOW_ID
2360 *spwidRF = ifno ;
2361
2362 // fill TIME, INTERVAL
2363 *timeRF = midTime ;
2364 *intervalRF = interval ;
2365
2366 // add row
2367 msSrc.addRow( 1, True ) ;
2368 row.put( nrow ) ;
2369
2370 iter2.next() ;
2371 }
2372
2373 iter1.next() ;
2374 }
2375
2376 // increment srcId if SRCNAME changed
2377 srcId++ ;
2378
2379 iter0.next() ;
2380 }
2381
2382// double endSec = mathutil::gettimeofday_sec() ;
2383// os_ << "end MSWriter::fillSource() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2384}
2385
2386void MSWriter::fillWeather()
2387{
2388// double startSec = mathutil::gettimeofday_sec() ;
2389// os_ << "start MSWriter::fillWeather() startSec=" << startSec << LogIO::POST ;
2390
2391 // access to MS WEATHER subtable
2392 MSWeather msw = mstable_->weather() ;
2393
2394 // access to WEATHER subtable
2395 Table stw = table_->weather().table() ;
2396 uInt nrow = stw.nrow() ;
2397
2398 if ( nrow == 0 )
2399 return ;
2400
2401 msw.addRow( nrow, True ) ;
2402 MSWeatherColumns mswCols( msw ) ;
2403
2404 // ANTENNA_ID is always 0
2405 Vector<Int> antIdArr( nrow, 0 ) ;
2406 mswCols.antennaId().putColumn( antIdArr ) ;
2407
2408 // fill weather status
2409 ROScalarColumn<Float> sharedFloatCol( stw, "TEMPERATURE" ) ;
2410 mswCols.temperature().putColumn( sharedFloatCol ) ;
2411 sharedFloatCol.attach( stw, "PRESSURE" ) ;
2412 mswCols.pressure().putColumn( sharedFloatCol ) ;
2413 sharedFloatCol.attach( stw, "HUMIDITY" ) ;
2414 mswCols.relHumidity().putColumn( sharedFloatCol ) ;
2415 sharedFloatCol.attach( stw, "WINDSPEED" ) ;
2416 mswCols.windSpeed().putColumn( sharedFloatCol ) ;
2417 sharedFloatCol.attach( stw, "WINDAZ" ) ;
2418 mswCols.windDirection().putColumn( sharedFloatCol ) ;
2419
2420 // fill TIME and INTERVAL
2421 Double midTime ;
2422 Double interval ;
2423 Vector<Double> intervalArr( nrow, 0.0 ) ;
2424 TableIterator iter( table_->table(), "WEATHER_ID" ) ;
2425 while( !iter.pastEnd() ) {
2426 //Table tab( iter.table() ) ;
2427 Table tab = iter.table() ;
2428
2429 ROScalarColumn<uInt> widCol( tab, "WEATHER_ID" ) ;
2430 uInt wid = widCol( 0 ) ;
2431
2432 getValidTimeRange( midTime, interval, tab ) ;
2433 mswCols.time().put( wid, midTime ) ;
2434 intervalArr[wid] = interval ;
2435
2436 iter.next() ;
2437 }
2438 mswCols.interval().putColumn( intervalArr ) ;
2439
2440// double endSec = mathutil::gettimeofday_sec() ;
2441// os_ << "end MSWriter::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2442}
2443
2444void MSWriter::fillSysCal()
2445{
2446 Table mssc = mstable_->sysCal() ;
2447
2448 {
2449 static const char *cols[] = {
2450 "BEAMNO", "IFNO", "TIME", "POLNO",
2451 NULL
2452 };
2453 static const TypeManagerImpl<uInt> tmUInt;
2454 static const TypeManagerImpl<Double> tmDouble;
2455 static const TypeManager *const tms[] = {
2456 &tmUInt, &tmUInt, &tmDouble, &tmUInt, NULL
2457 };
2458 //double t0 = mathutil::gettimeofday_sec() ;
2459 MSSysCalVisitor myVisitor(table_->table(),mssc);
2460 //double t1 = mathutil::gettimeofday_sec() ;
2461 //cout << "MSWriterVisitor(): elapsed time " << t1-t0 << " sec" << endl ;
2462 traverseTable(table_->table(), cols, tms, &myVisitor);
2463 //double t3 = mathutil::gettimeofday_sec() ;
2464 //cout << "traverseTable(): elapsed time " << t3-t2 << " sec" << endl ;
2465 }
2466
2467}
2468
2469void MSWriter::getValidTimeRange( Double &me, Double &interval, Table &tab )
2470{
2471// double startSec = mathutil::gettimeofday_sec() ;
2472// os_ << "start MSWriter::getVaridTimeRange() startSec=" << startSec << LogIO::POST ;
2473
2474 // sort table
2475 //Table stab = tab.sort( "TIME" ) ;
2476
2477 ROScalarColumn<Double> timeCol( tab, "TIME" ) ;
2478 Vector<Double> timeArr = timeCol.getColumn() ;
2479 Double minTime ;
2480 Double maxTime ;
2481 minMax( minTime, maxTime, timeArr ) ;
2482 Double midTime = 0.5 * ( minTime + maxTime ) * 86400.0 ;
2483 // unit for TIME
2484 // Scantable: "d"
2485 // MS: "s"
2486 me = midTime ;
2487 interval = ( maxTime - minTime ) * 86400.0 ;
2488
2489// double endSec = mathutil::gettimeofday_sec() ;
2490// os_ << "end MSWriter::getValidTimeRange() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2491}
2492
2493void MSWriter::getValidTimeRange( Double &me, Double &interval, Vector<Double> &atime, Vector<Double> &ainterval )
2494{
2495// double startSec = mathutil::gettimeofday_sec() ;
2496// os_ << "start MSWriter::getVaridTimeRange() startSec=" << startSec << LogIO::POST ;
2497
2498 // sort table
2499 //Table stab = tab.sort( "TIME" ) ;
2500
2501 Double minTime ;
2502 Double maxTime ;
2503 minMax( minTime, maxTime, atime ) ;
2504 Double midTime = 0.5 * ( minTime + maxTime ) * 86400.0 ;
2505 // unit for TIME
2506 // Scantable: "d"
2507 // MS: "s"
2508 me = midTime ;
2509 interval = ( maxTime - minTime ) * 86400.0 + mean( ainterval ) ;
2510
2511// double endSec = mathutil::gettimeofday_sec() ;
2512// os_ << "end MSWriter::getValidTimeRange() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2513}
2514
2515}
Note: See TracBrowser for help on using the repository browser.