source: trunk/src/MSWriter.cpp@ 2292

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

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: MSFillerUtils and MSWriterUtils added

Test Programs: sd regressions, test_sdsave

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Rewrote MSFiller and MSWriter based on TableTraverse.
TableTraverse is changed a bit since it doesn't work on plain table.


File size: 92.7 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;
183 Int lastSrcType;
184 uInt lastCycleNo;
185 Double lastTime;
186 Int lastPolNo;
187protected:
188 const Table &table;
189 uInt count;
190public:
191 BaseMSWriterVisitor(const Table &table)
192 : table(table)
193 {
194 static const String dummy;
195 lastFieldName = &dummy;
196 count = 0;
197 }
198
199 virtual void enterFieldName(const uInt recordNo, const String &columnValue) {
200 }
201 virtual void leaveFieldName(const uInt recordNo, const String &columnValue) {
202 }
203 virtual void enterBeamNo(const uInt recordNo, uInt columnValue) { }
204 virtual void leaveBeamNo(const uInt recordNo, uInt columnValue) { }
205 virtual void enterScanNo(const uInt recordNo, uInt columnValue) { }
206 virtual void leaveScanNo(const uInt recordNo, uInt columnValue) { }
207 virtual void enterIfNo(const uInt recordNo, uInt columnValue) { }
208 virtual void leaveIfNo(const uInt recordNo, uInt columnValue) { }
209 virtual void enterSrcType(const uInt recordNo, Int columnValue) { }
210 virtual void leaveSrcType(const uInt recordNo, Int columnValue) { }
211 virtual void enterCycleNo(const uInt recordNo, uInt columnValue) { }
212 virtual void leaveCycleNo(const uInt recordNo, uInt columnValue) { }
213 virtual void enterTime(const uInt recordNo, Double columnValue) { }
214 virtual void leaveTime(const uInt recordNo, Double columnValue) { }
215 virtual void enterPolNo(const uInt recordNo, uInt columnValue) { }
216 virtual void leavePolNo(const uInt recordNo, uInt columnValue) { }
217
218 virtual Bool visitRecord(const uInt recordNo,
219 const String &fieldName,
220 const uInt beamNo,
221 const uInt scanNo,
222 const uInt ifNo,
223 const Int srcType,
224 const uInt cycleNo,
225 const Double time,
226 const uInt polNo) { return True ;}
227
228 virtual Bool visit(Bool isFirst, const uInt recordNo,
229 const uInt nCols, void const *const colValues[]) {
230 const String *fieldName = NULL;
231 uInt beamNo, scanNo, ifNo;
232 Int srcType;
233 uInt cycleNo;
234 Double time;
235 Int polNo;
236 { // prologue
237 uInt i = 0;
238 {
239 const String *col = (const String*)colValues[i++];
240 fieldName = &col[recordNo];
241 }
242 {
243 const uInt *col = (const uInt *)colValues[i++];
244 beamNo = col[recordNo];
245 }
246 {
247 const uInt *col = (const uInt *)colValues[i++];
248 scanNo = col[recordNo];
249 }
250 {
251 const uInt *col = (const uInt *)colValues[i++];
252 ifNo = col[recordNo];
253 }
254 {
255 const Int *col = (const Int *)colValues[i++];
256 srcType = col[recordNo];
257 }
258 {
259 const uInt *col = (const uInt *)colValues[i++];
260 cycleNo = col[recordNo];
261 }
262 {
263 const Double *col = (const Double *)colValues[i++];
264 time = col[recordNo];
265 }
266 {
267 const Int *col = (const Int *)colValues[i++];
268 polNo = col[recordNo];
269 }
270 assert(nCols == i);
271 }
272
273 if (isFirst) {
274 enterFieldName(recordNo, *fieldName);
275 enterBeamNo(recordNo, beamNo);
276 enterScanNo(recordNo, scanNo);
277 enterIfNo(recordNo, ifNo);
278 enterSrcType(recordNo, srcType);
279 enterCycleNo(recordNo, cycleNo);
280 enterTime(recordNo, time);
281 enterPolNo(recordNo, polNo);
282 } else {
283 if (lastFieldName->compare(*fieldName) != 0) {
284 leavePolNo(lastRecordNo, lastPolNo);
285 leaveTime(lastRecordNo, lastTime);
286 leaveCycleNo(lastRecordNo, lastCycleNo);
287 leaveSrcType(lastRecordNo, lastSrcType);
288 leaveIfNo(lastRecordNo, lastIfNo);
289 leaveScanNo(lastRecordNo, lastScanNo);
290 leaveBeamNo(lastRecordNo, lastBeamNo);
291 leaveFieldName(lastRecordNo, *lastFieldName);
292
293 enterFieldName(recordNo, *fieldName);
294 enterBeamNo(recordNo, beamNo);
295 enterScanNo(recordNo, scanNo);
296 enterIfNo(recordNo, ifNo);
297 enterSrcType(recordNo, srcType);
298 enterCycleNo(recordNo, cycleNo);
299 enterTime(recordNo, time);
300 enterPolNo(recordNo, polNo);
301 } else if (lastBeamNo != beamNo) {
302 leavePolNo(lastRecordNo, lastPolNo);
303 leaveTime(lastRecordNo, lastTime);
304 leaveCycleNo(lastRecordNo, lastCycleNo);
305 leaveSrcType(lastRecordNo, lastSrcType);
306 leaveIfNo(lastRecordNo, lastIfNo);
307 leaveScanNo(lastRecordNo, lastScanNo);
308 leaveBeamNo(lastRecordNo, lastBeamNo);
309
310 enterBeamNo(recordNo, beamNo);
311 enterScanNo(recordNo, scanNo);
312 enterIfNo(recordNo, ifNo);
313 enterSrcType(recordNo, srcType);
314 enterCycleNo(recordNo, cycleNo);
315 enterTime(recordNo, time);
316 enterPolNo(recordNo, polNo);
317 } else if (lastScanNo != scanNo) {
318 leavePolNo(lastRecordNo, lastPolNo);
319 leaveTime(lastRecordNo, lastTime);
320 leaveCycleNo(lastRecordNo, lastCycleNo);
321 leaveSrcType(lastRecordNo, lastSrcType);
322 leaveIfNo(lastRecordNo, lastIfNo);
323 leaveScanNo(lastRecordNo, lastScanNo);
324
325 enterScanNo(recordNo, scanNo);
326 enterIfNo(recordNo, ifNo);
327 enterSrcType(recordNo, srcType);
328 enterCycleNo(recordNo, cycleNo);
329 enterTime(recordNo, time);
330 enterPolNo(recordNo, polNo);
331 } else if (lastIfNo != ifNo) {
332 leavePolNo(lastRecordNo, lastPolNo);
333 leaveTime(lastRecordNo, lastTime);
334 leaveCycleNo(lastRecordNo, lastCycleNo);
335 leaveSrcType(lastRecordNo, lastSrcType);
336 leaveIfNo(lastRecordNo, lastIfNo);
337
338 enterIfNo(recordNo, ifNo);
339 enterSrcType(recordNo, srcType);
340 enterCycleNo(recordNo, cycleNo);
341 enterTime(recordNo, time);
342 enterPolNo(recordNo, polNo);
343 } else if (lastSrcType != srcType) {
344 leavePolNo(lastRecordNo, lastPolNo);
345 leaveTime(lastRecordNo, lastTime);
346 leaveCycleNo(lastRecordNo, lastCycleNo);
347 leaveSrcType(lastRecordNo, lastSrcType);
348
349 enterSrcType(recordNo, srcType);
350 enterCycleNo(recordNo, cycleNo);
351 enterTime(recordNo, time);
352 enterPolNo(recordNo, polNo);
353 } else if (lastCycleNo != cycleNo) {
354 leavePolNo(lastRecordNo, lastPolNo);
355 leaveTime(lastRecordNo, lastTime);
356 leaveCycleNo(lastRecordNo, lastCycleNo);
357
358 enterCycleNo(recordNo, cycleNo);
359 enterTime(recordNo, time);
360 enterPolNo(recordNo, polNo);
361 } else if (lastTime != time) {
362 leavePolNo(lastRecordNo, lastPolNo);
363 leaveTime(lastRecordNo, lastTime);
364
365 enterTime(recordNo, time);
366 enterPolNo(recordNo, polNo);
367 } else if (lastPolNo != polNo) {
368 leavePolNo(lastRecordNo, lastPolNo);
369 enterPolNo(recordNo, polNo);
370 }
371 }
372 count++;
373 Bool result = visitRecord(recordNo, *fieldName, beamNo, scanNo, ifNo, srcType,
374 cycleNo, time, polNo);
375
376 { // epilogue
377 lastRecordNo = recordNo;
378
379 lastFieldName = fieldName;
380 lastBeamNo = beamNo;
381 lastScanNo = scanNo;
382 lastIfNo = ifNo;
383 lastSrcType = srcType;
384 lastCycleNo = cycleNo;
385 lastTime = time;
386 lastPolNo = polNo;
387 }
388 return result ;
389 }
390
391 virtual void finish() {
392 if (count > 0) {
393 leavePolNo(lastRecordNo, lastPolNo);
394 leaveTime(lastRecordNo, lastTime);
395 leaveCycleNo(lastRecordNo, lastCycleNo);
396 leaveSrcType(lastRecordNo, lastSrcType);
397 leaveIfNo(lastRecordNo, lastIfNo);
398 leaveScanNo(lastRecordNo, lastScanNo);
399 leaveBeamNo(lastRecordNo, lastBeamNo);
400 leaveFieldName(lastRecordNo, *lastFieldName);
401 }
402 }
403};
404
405class MSWriterVisitor: public BaseMSWriterVisitor, public MSWriterUtils {
406public:
407 MSWriterVisitor(const Table &table, Table &mstable)
408 : BaseMSWriterVisitor(table),
409 ms(mstable)
410 {
411 rowidx = 0 ;
412 fieldName = "" ;
413 defaultFieldId = 0 ;
414 spwId = -1 ;
415 subscan = 1 ;
416 ptName = "" ;
417 srcId = 0 ;
418
419 row = TableRow( ms ) ;
420
421 holder.reset() ;
422
423 makePolMap() ;
424 initFrequencies() ;
425 initTcal() ;
426
427 //
428 // add rows to MS
429 //
430 uInt addrow = table.nrow() ;
431 ms.addRow( addrow ) ;
432
433 // attach to Scantable columns
434 spectraCol.attach( table, "SPECTRA" ) ;
435 flagtraCol.attach( table, "FLAGTRA" ) ;
436 flagRowCol.attach( table, "FLAGROW" ) ;
437 tcalIdCol.attach( table, "TCAL_ID" ) ;
438 intervalCol.attach( table, "INTERVAL" ) ;
439 directionCol.attach( table, "DIRECTION" ) ;
440 scanRateCol.attach( table, "SCANRATE" ) ;
441 timeCol.attach( table, "TIME" ) ;
442 freqIdCol.attach( table, "FREQ_ID" ) ;
443 sourceNameCol.attach( table, "SRCNAME" ) ;
444 sourceDirectionCol.attach( table, "SRCDIRECTION" ) ;
445 fieldNameCol.attach( table, "FIELDNAME" ) ;
446
447 // MS subtables
448 attachSubtables() ;
449
450 // attach to MS columns
451 attachMain() ;
452 attachPointing() ;
453 }
454
455 virtual void enterFieldName(const uInt recordNo, const String &columnValue) {
456 //printf("%u: FieldName: %s\n", recordNo, columnValue.c_str());
457 fieldName = fieldNameCol.asString( recordNo ) ;
458 String::size_type pos = fieldName.find( "__" ) ;
459 if ( pos != String::npos ) {
460 fieldId = String::toInt( fieldName.substr( pos+2 ) ) ;
461 fieldName = fieldName.substr( 0, pos ) ;
462 }
463 else {
464 fieldId = defaultFieldId ;
465 defaultFieldId++ ;
466 }
467 Double tSec = timeCol.asdouble( recordNo ) * 86400.0 ;
468 Vector<Double> srcDir = sourceDirectionCol( recordNo ) ;
469 Vector<Double> srate = scanRateCol( recordNo ) ;
470 String srcName = sourceNameCol.asString( recordNo ) ;
471
472 addField( fieldId, fieldName, srcName, srcDir, srate, tSec ) ;
473
474 // put value
475 *fieldIdRF = fieldId ;
476 }
477 virtual void leaveFieldName(const uInt recordNo, const String &columnValue) {
478 }
479 virtual void enterBeamNo(const uInt recordNo, uInt columnValue) {
480 //printf("%u: BeamNo: %u\n", recordNo, columnValue);
481
482 feedId = (Int)columnValue ;
483
484 // put value
485 *feed1RF = feedId ;
486 *feed2RF = feedId ;
487 }
488 virtual void leaveBeamNo(const uInt recordNo, uInt columnValue) {
489 }
490 virtual void enterScanNo(const uInt recordNo, uInt columnValue) {
491 //printf("%u: ScanNo: %u\n", recordNo, columnValue);
492
493 // put value
494 // SCAN_NUMBER is 0-based in Scantable while 1-based in MS
495 *scanNumberRF = (Int)columnValue + 1 ;
496 }
497 virtual void leaveScanNo(const uInt recordNo, uInt columnValue) {
498 subscan = 1 ;
499 }
500 virtual void enterIfNo(const uInt recordNo, uInt columnValue) {
501 //printf("%u: IfNo: %u\n", recordNo, columnValue);
502
503 spwId = (Int)columnValue ;
504 uInt freqId = freqIdCol.asuInt( recordNo ) ;
505
506 Vector<Float> sp = spectraCol( recordNo ) ;
507 uInt nchan = sp.nelements() ;
508 holder.setNchan( nchan ) ;
509
510 addSpectralWindow( spwId, freqId ) ;
511
512 addFeed( feedId, spwId ) ;
513 }
514 virtual void leaveIfNo(const uInt recordNo, uInt columnValue) {
515 }
516 virtual void enterSrcType(const uInt recordNo, Int columnValue) {
517 //printf("%u: SrcType: %d\n", recordNo, columnValue);
518
519 Int stateId = addState( columnValue ) ;
520
521 // put value
522 *stateIdRF = stateId ;
523 }
524 virtual void leaveSrcType(const uInt recordNo, Int columnValue) {
525 }
526 virtual void enterCycleNo(const uInt recordNo, uInt columnValue) {
527 //printf("%u: CycleNo: %u\n", recordNo, columnValue);
528 }
529 virtual void leaveCycleNo(const uInt recordNo, uInt columnValue) {
530 }
531 virtual void enterTime(const uInt recordNo, Double columnValue) {
532 //printf("%u: Time: %f\n", recordNo, columnValue);
533
534 Double timeSec = columnValue * 86400.0 ;
535 Double interval = intervalCol.asdouble( recordNo ) ;
536
537 if ( ptName.empty() ) {
538 Vector<Double> dir = directionCol( recordNo ) ;
539 Vector<Double> rate = scanRateCol( recordNo ) ;
540 if ( anyNE( rate, 0.0 ) ) {
541 Matrix<Double> msdir( 2, 2 ) ;
542 msdir.column( 0 ) = dir ;
543 msdir.column( 1 ) = rate ;
544 addPointing( timeSec, interval, msdir ) ;
545 }
546 else {
547 Matrix<Double> msdir( 2, 1 ) ;
548 msdir.column( 0 ) = dir ;
549 addPointing( timeSec, interval, msdir ) ;
550 }
551 }
552
553 // put value
554 *timeRF = timeSec ;
555 *timeCentroidRF = timeSec ;
556 *intervalRF = interval ;
557 *exposureRF = interval ;
558 }
559 virtual void leaveTime(const uInt recordNo, Double columnValue) {
560 if ( holder.nPol() > 0 ) {
561 Vector<Float> w = holder.getWeight() ;
562 Cube<Bool> c = holder.getFlagCategory() ;
563 Bool flr = holder.getFlagRow() ;
564 Matrix<Bool> fl = holder.getFlag() ;
565 Vector<uInt> polnos = holder.polNos() ;
566 Int polId = addPolarization( polnos ) ;
567 Int ddId = addDataDescription( polId, spwId ) ;
568
569 // put field
570 *dataDescIdRF = ddId ;
571 *flagRowRF = flr ;
572 weightRF.define( w ) ;
573 sigmaRF.define( w ) ;
574 flagCategoryRF.define( c ) ;
575 flagRF.define( fl ) ;
576 if ( useFloat ) {
577 Matrix<Float> sp = holder.getData() ;
578 floatDataRF.define( sp ) ;
579 }
580 else {
581 Matrix<Complex> sp = holder.getComplexData() ;
582 dataRF.define( sp ) ;
583 }
584
585 // commit row
586 row.put( rowidx ) ;
587 rowidx++ ;
588
589 // reset holder
590 holder.reset() ;
591 }
592 if ( tcalKey != -1 ) {
593 tcalNotYet[tcalKey] = False ;
594 tcalKey = -1 ;
595 }
596 }
597 virtual void enterPolNo(const uInt recordNo, uInt columnValue) {
598 //printf("%u: PolNo: %d\n", recordNo, columnValue);
599 uInt tcalId = tcalIdCol.asuInt( recordNo ) ;
600 if ( tcalKey == -1 ) {
601 tcalKey = tcalId ;
602 }
603 if ( tcalNotYet[tcalKey] ) {
604 map< Int,Vector<uInt> >::iterator itr = tcalIdRec.find( tcalKey ) ;
605 if ( itr != tcalIdRec.end() ) {
606 Vector<uInt> ids = itr->second ;
607 uInt nrow = ids.nelements() ;
608 ids.resize( nrow+1, True ) ;
609 ids[nrow] = tcalId ;
610 tcalIdRec.erase( tcalKey ) ;
611 tcalIdRec[tcalKey] = ids ;
612 }
613 else {
614 Vector<uInt> rows( 1, tcalId ) ;
615 tcalIdRec[tcalKey] = rows ;
616 }
617 }
618 map< Int,Vector<uInt> >::iterator itr = tcalRowRec.find( tcalKey ) ;
619 if ( itr != tcalRowRec.end() ) {
620 Vector<uInt> rows = itr->second ;
621 uInt nrow = rows.nelements() ;
622 rows.resize( nrow+1, True ) ;
623 rows[nrow] = recordNo ;
624 tcalRowRec.erase( tcalKey ) ;
625 tcalRowRec[tcalKey] = rows ;
626 }
627 else {
628 Vector<uInt> rows( 1, recordNo ) ;
629 tcalRowRec[tcalKey] = rows ;
630 }
631 }
632 virtual void leavePolNo(const uInt recordNo, uInt columnValue) {
633 }
634
635 virtual Bool visitRecord(const uInt recordNo,
636 const String &fieldName,
637 const uInt beamNo,
638 const uInt scanNo,
639 const uInt ifNo,
640 const Int srcType,
641 const uInt cycleNo,
642 const Double time,
643 const uInt polNo) {
644 //printf("%u: %s, %u, %u, %u, %d, %u, %f, %d\n", recordNo,
645 // fieldName.c_str(), beamNo, scanNo, ifNo, srcType, cycleNo, time, polNo);
646
647 Vector<Float> sp = spectraCol( recordNo ) ;
648 Vector<uChar> tmp = flagtraCol( recordNo ) ;
649 Vector<Bool> fl( tmp.shape() ) ;
650 convertArray( fl, tmp ) ;
651 Bool flr = (Bool)flagRowCol.asuInt( recordNo ) ;
652 holder.accumulate( polNo, sp, fl, flr ) ;
653
654 return True ;
655 }
656
657 virtual void finish() {
658 BaseMSWriterVisitor::finish();
659 //printf("Total: %u\n", count);
660
661 // remove rows
662 if ( ms.nrow() > rowidx ) {
663 uInt numRemove = ms.nrow() - rowidx ;
664 //cout << "numRemove = " << numRemove << endl ;
665 Vector<uInt> rows( numRemove ) ;
666 indgen( rows, rowidx ) ;
667 ms.removeRow( rows ) ;
668 }
669
670 // fill empty SPECTRAL_WINDOW rows
671 infillSpectralWindow() ;
672 }
673
674 void dataColumnName( String name )
675 {
676 TableRecord &r = row.record() ;
677 if ( name == "DATA" ) {
678 useFloat = False ;
679 dataRF.attachToRecord( r, name ) ;
680 }
681 else if ( name == "FLOAT_DATA" ) {
682 useFloat = True ;
683 floatDataRF.attachToRecord( r, name ) ;
684 }
685 }
686 void pointingTableName( String name ) {
687 ptName = name ;
688 }
689 void setSourceRecord( Record &r ) {
690 srcRec = r ;
691 }
692 map< Int,Vector<uInt> > &getTcalIdRecord() { return tcalIdRec ; }
693 map< Int,Vector<uInt> > &getTcalRowRecord() { return tcalRowRec ; }
694private:
695 void addField( Int &fid, String &fname, String &srcName,
696 Vector<Double> &sdir, Vector<Double> &srate,
697 Double &tSec )
698 {
699 uInt nrow = fieldtab.nrow() ;
700 while( (Int)nrow <= fid ) {
701 fieldtab.addRow( 1, True ) ;
702 nrow++ ;
703 }
704
705 Matrix<Double> dir ;
706 Int numPoly = 0 ;
707 if ( anyNE( srate, 0.0 ) ) {
708 dir.resize( 2, 2 ) ;
709 dir.column( 0 ) = sdir ;
710 dir.column( 1 ) = srate ;
711 numPoly = 1 ;
712 }
713 else {
714 dir.resize( 2, 1 ) ;
715 dir.column( 0 ) = sdir ;
716 }
717 srcId = srcRec.asInt( srcName ) ;
718
719 TableRow tr( fieldtab ) ;
720 TableRecord &r = tr.record() ;
721 putField( "NAME", r, fname ) ;
722 putField( "NUM_POLY", r, numPoly ) ;
723 putField( "TIME", r, tSec ) ;
724 putField( "SOURCE_ID", r, srcId ) ;
725 defineField( "DELAY_DIR", r, dir ) ;
726 defineField( "REFERENCE_DIR", r, dir ) ;
727 defineField( "PHASE_DIR", r, dir ) ;
728 tr.put( fid ) ;
729
730 // for POINTING table
731 *poNameRF = fname ;
732 }
733 Int addState( Int &id )
734 {
735 String obsMode ;
736 Bool isSignal ;
737 Double tnoise ;
738 Double tload ;
739 queryType( id, obsMode, isSignal, tnoise, tload ) ;
740
741 String key = obsMode+"_"+String::toString( subscan ) ;
742 Int idx = -1 ;
743 uInt nEntry = stateEntry.nelements() ;
744 for ( uInt i = 0 ; i < nEntry ; i++ ) {
745 if ( stateEntry[i] == key ) {
746 idx = i ;
747 break ;
748 }
749 }
750 if ( idx == -1 ) {
751 uInt nrow = statetab.nrow() ;
752 statetab.addRow( 1, True ) ;
753 TableRow tr( statetab ) ;
754 TableRecord &r = tr.record() ;
755 putField( "OBS_MODE", r, obsMode ) ;
756 putField( "SIG", r, isSignal ) ;
757 isSignal = !isSignal ;
758 putField( "REF", r, isSignal ) ;
759 putField( "CAL", r, tnoise ) ;
760 putField( "LOAD", r, tload ) ;
761 tr.put( nrow ) ;
762 idx = nrow ;
763
764 stateEntry.resize( nEntry+1, True ) ;
765 stateEntry[nEntry] = key ;
766 }
767 subscan++ ;
768
769 return idx ;
770 }
771 void addPointing( Double &tSec, Double &interval, Matrix<Double> &dir )
772 {
773 uInt nrow = potab.nrow() ;
774 potab.addRow( 1, True ) ;
775
776 *poNumPolyRF = dir.ncolumn() - 1 ;
777 *poTimeRF = tSec ;
778 *poTimeOriginRF = tSec ;
779 *poIntervalRF = interval ;
780 poDirectionRF.define( dir ) ;
781 poTargetRF.define( dir ) ;
782 porow.put( nrow ) ;
783 }
784 Int addPolarization( Vector<uInt> &nos )
785 {
786 Int idx = -1 ;
787 uInt nEntry = polEntry.size() ;
788 for ( uInt i = 0 ; i < nEntry ; i++ ) {
789 if ( polEntry[i].conform( nos ) && allEQ( polEntry[i], nos ) ) {
790 idx = i ;
791 break ;
792 }
793 }
794
795 Int numCorr ;
796 Vector<Int> corrType ;
797 Matrix<Int> corrProduct ;
798 polProperty( nos, numCorr, corrType, corrProduct ) ;
799
800 if ( idx == -1 ) {
801 uInt nrow = poltab.nrow() ;
802 poltab.addRow( 1, True ) ;
803 TableRow tr( poltab ) ;
804 TableRecord &r = tr.record() ;
805 putField( "NUM_CORR", r, numCorr ) ;
806 defineField( "CORR_TYPE", r, corrType ) ;
807 defineField( "CORR_PRODUCT", r, corrProduct ) ;
808 tr.put( nrow ) ;
809 idx = nrow ;
810
811 polEntry.resize( nEntry+1 ) ;
812 polEntry[nEntry] = nos ;
813 }
814
815 return idx ;
816 }
817 Int addDataDescription( Int pid, Int sid )
818 {
819 Int idx = -1 ;
820 uInt nEntry = ddEntry.nrow() ;
821 Vector<Int> key( 2 ) ;
822 key[0] = pid ;
823 key[1] = sid ;
824 for ( uInt i = 0 ; i < nEntry ; i++ ) {
825 if ( allEQ( ddEntry.row(i), key ) ) {
826 idx = i ;
827 break ;
828 }
829 }
830
831 if ( idx == -1 ) {
832 uInt nrow = ddtab.nrow() ;
833 ddtab.addRow( 1, True ) ;
834 TableRow tr( ddtab ) ;
835 TableRecord &r = tr.record() ;
836 putField( "POLARIZATION_ID", r, pid ) ;
837 putField( "SPECTRAL_WINDOW_ID", r, sid ) ;
838 tr.put( nrow ) ;
839 idx = nrow ;
840
841 ddEntry.resize( nEntry+1, 2, True ) ;
842 ddEntry.row(nEntry) = key ;
843 }
844
845 return idx ;
846 }
847 void infillSpectralWindow()
848 {
849 ROScalarColumn<Int> nchanCol( spwtab, "NUM_CHAN" ) ;
850 Vector<Int> nchan = nchanCol.getColumn() ;
851 TableRow tr( spwtab ) ;
852 TableRecord &r = tr.record() ;
853 Int mfr = 1 ;
854 Vector<Double> dummy( 1, 0.0 ) ;
855 putField( "MEAS_FREQ_REF", r, mfr ) ;
856 defineField( "CHAN_FREQ", r, dummy ) ;
857 defineField( "CHAN_WIDTH", r, dummy ) ;
858 defineField( "EFFECTIVE_BW", r, dummy ) ;
859 defineField( "RESOLUTION", r, dummy ) ;
860
861 for ( uInt i = 0 ; i < spwtab.nrow() ; i++ ) {
862 if ( nchan[i] == 0 )
863 tr.put( i ) ;
864 }
865 }
866 void addSpectralWindow( Int sid, uInt fid )
867 {
868 if ( !processedFreqId[fid] ) {
869 uInt nrow = spwtab.nrow() ;
870 while( (Int)nrow <= sid ) {
871 spwtab.addRow( 1, True ) ;
872 nrow++ ;
873 }
874 processedFreqId[fid] = True ;
875 }
876
877 Double rp = refpix[fid] ;
878 Double rv = refval[fid] ;
879 Double ic = increment[fid] ;
880
881 Int mfrInt = (Int)freqframe ;
882 Int nchan = holder.nChan() ;
883 Double bw = nchan * abs( ic ) ;
884 Double reffreq = rv - rp * ic ;
885 Int netsb = 0 ; // USB->0, LSB->1
886 if ( ic < 0 )
887 netsb = 1 ;
888 Vector<Double> res( nchan, abs(ic) ) ;
889 Vector<Double> chanf( nchan ) ;
890 indgen( chanf, reffreq, ic ) ;
891
892 TableRow tr( spwtab ) ;
893 TableRecord &r = tr.record() ;
894 putField( "MEAS_FREQ_REF", r, mfrInt ) ;
895 putField( "NUM_CHAN", r, nchan ) ;
896 putField( "TOTAL_BANDWIDTH", r, bw ) ;
897 putField( "REF_FREQUENCY", r, reffreq ) ;
898 putField( "NET_SIDEBAND", r, netsb ) ;
899 defineField( "RESOLUTION", r, res ) ;
900 defineField( "CHAN_WIDTH", r, res ) ;
901 defineField( "EFFECTIVE_BW", r, res ) ;
902 defineField( "CHAN_FREQ", r, chanf ) ;
903 tr.put( sid ) ;
904 }
905 void addFeed( Int fid, Int sid )
906 {
907 Int idx = -1 ;
908 uInt nEntry = feedEntry.nrow() ;
909 Vector<Int> key( 2 ) ;
910 key[0] = fid ;
911 key[1] = sid ;
912 for ( uInt i = 0 ; i < nEntry ; i++ ) {
913 if ( allEQ( feedEntry.row(i), key ) ) {
914 idx = i ;
915 break ;
916 }
917 }
918
919 if ( idx == -1 ) {
920 uInt nrow = feedtab.nrow() ;
921 feedtab.addRow( 1, True ) ;
922 Int numReceptors = 2 ;
923 Vector<String> polType( numReceptors ) ;
924 Matrix<Double> beamOffset( 2, numReceptors, 0.0 ) ;
925 Vector<Double> receptorAngle( numReceptors, 0.0 ) ;
926 if ( poltype == "linear" ) {
927 polType[0] = "X" ;
928 polType[1] = "Y" ;
929 }
930 else if ( poltype == "circular" ) {
931 polType[0] = "R" ;
932 polType[1] = "L" ;
933 }
934 else {
935 polType[0] = "X" ;
936 polType[1] = "Y" ;
937 }
938 Matrix<Complex> polResponse( numReceptors, numReceptors, 0.0 ) ;
939
940 TableRow tr( feedtab ) ;
941 TableRecord &r = tr.record() ;
942 putField( "FEED_ID", r, fid ) ;
943 putField( "BEAM_ID", r, fid ) ;
944 Int tmp = 0 ;
945 putField( "ANTENNA_ID", r, tmp ) ;
946 putField( "SPECTRAL_WINDOW_ID", r, sid ) ;
947 putField( "NUM_RECEPTORS", r, numReceptors ) ;
948 defineField( "POLARIZATION_TYPE", r, polType ) ;
949 defineField( "BEAM_OFFSET", r, beamOffset ) ;
950 defineField( "RECEPTOR_ANGLE", r, receptorAngle ) ;
951 defineField( "POL_RESPONSE", r, polResponse ) ;
952 tr.put( nrow ) ;
953
954 feedEntry.resize( nEntry+1, 2, True ) ;
955 feedEntry.row( nEntry ) = key ;
956 }
957 }
958 void makePolMap()
959 {
960 const TableRecord &keys = table.keywordSet() ;
961 poltype = keys.asString( "POLTYPE" ) ;
962
963 if ( poltype == "stokes" ) {
964 polmap.resize( 4 ) ;
965 polmap[0] = Stokes::I ;
966 polmap[1] = Stokes::Q ;
967 polmap[2] = Stokes::U ;
968 polmap[3] = Stokes::V ;
969 }
970 else if ( poltype == "linear" ) {
971 polmap.resize( 4 ) ;
972 polmap[0] = Stokes::XX ;
973 polmap[1] = Stokes::YY ;
974 polmap[2] = Stokes::XY ;
975 polmap[3] = Stokes::YX ;
976 }
977 else if ( poltype == "circular" ) {
978 polmap.resize( 4 ) ;
979 polmap[0] = Stokes::RR ;
980 polmap[1] = Stokes::LL ;
981 polmap[2] = Stokes::RL ;
982 polmap[3] = Stokes::LR ;
983 }
984 else if ( poltype == "linpol" ) {
985 polmap.resize( 2 ) ;
986 polmap[0] = Stokes::Plinear ;
987 polmap[1] = Stokes::Pangle ;
988 }
989 else {
990 polmap.resize( 0 ) ;
991 }
992 }
993 void initFrequencies()
994 {
995 const TableRecord &keys = table.keywordSet() ;
996 Table tab = keys.asTable( "FREQUENCIES" ) ;
997 ROScalarColumn<uInt> idcol( tab, "ID" ) ;
998 ROScalarColumn<Double> rpcol( tab, "REFPIX" ) ;
999 ROScalarColumn<Double> rvcol( tab, "REFVAL" ) ;
1000 ROScalarColumn<Double> iccol( tab, "INCREMENT" ) ;
1001 Vector<uInt> id = idcol.getColumn() ;
1002 Vector<Double> rp = rpcol.getColumn() ;
1003 Vector<Double> rv = rvcol.getColumn() ;
1004 Vector<Double> ic = iccol.getColumn() ;
1005 for ( uInt i = 0 ; i < id.nelements() ; i++ ) {
1006 processedFreqId.insert( pair<uInt,Bool>( id[i], False ) ) ;
1007 refpix.insert( pair<uInt,Double>( id[i], rp[i] ) ) ;
1008 refval.insert( pair<uInt,Double>( id[i], rv[i] ) ) ;
1009 increment.insert( pair<uInt,Double>( id[i], ic[i] ) ) ;
1010 }
1011 String frameStr = tab.keywordSet().asString( "BASEFRAME" ) ;
1012 MFrequency::getType( freqframe, frameStr ) ;
1013 }
1014 void attachSubtables()
1015 {
1016 const TableRecord &keys = table.keywordSet() ;
1017 TableRecord &mskeys = ms.rwKeywordSet() ;
1018
1019 // FIELD table
1020 fieldtab = mskeys.asTable( "FIELD" ) ;
1021
1022 // SPECTRAL_WINDOW table
1023 spwtab = mskeys.asTable( "SPECTRAL_WINDOW" ) ;
1024
1025 // POINTING table
1026 potab = mskeys.asTable( "POINTING" ) ;
1027
1028 // POLARIZATION table
1029 poltab = mskeys.asTable( "POLARIZATION" ) ;
1030
1031 // DATA_DESCRIPTION table
1032 ddtab = mskeys.asTable( "DATA_DESCRIPTION" ) ;
1033
1034 // STATE table
1035 statetab = mskeys.asTable( "STATE" ) ;
1036
1037 // FEED table
1038 feedtab = mskeys.asTable( "FEED" ) ;
1039 }
1040 void attachMain()
1041 {
1042 TableRecord &r = row.record() ;
1043 dataDescIdRF.attachToRecord( r, "DATA_DESC_ID" ) ;
1044 flagRowRF.attachToRecord( r, "FLAG_ROW" ) ;
1045 weightRF.attachToRecord( r, "WEIGHT" ) ;
1046 sigmaRF.attachToRecord( r, "SIGMA" ) ;
1047 flagCategoryRF.attachToRecord( r, "FLAG_CATEGORY" ) ;
1048 flagRF.attachToRecord( r, "FLAG" ) ;
1049 timeRF.attachToRecord( r, "TIME" ) ;
1050 timeCentroidRF.attachToRecord( r, "TIME_CENTROID" ) ;
1051 intervalRF.attachToRecord( r, "INTERVAL" ) ;
1052 exposureRF.attachToRecord( r, "EXPOSURE" ) ;
1053 fieldIdRF.attachToRecord( r, "FIELD_ID" ) ;
1054 feed1RF.attachToRecord( r, "FEED1" ) ;
1055 feed2RF.attachToRecord( r, "FEED2" ) ;
1056 scanNumberRF.attachToRecord( r, "SCAN_NUMBER" ) ;
1057 stateIdRF.attachToRecord( r, "STATE_ID" ) ;
1058
1059 // constant values
1060 Int id = 0 ;
1061 RecordFieldPtr<Int> intRF( r, "OBSERVATION_ID" ) ;
1062 *intRF = 0 ;
1063 intRF.attachToRecord( r, "ANTENNA1" ) ;
1064 *intRF = 0 ;
1065 intRF.attachToRecord( r, "ANTENNA2" ) ;
1066 *intRF = 0 ;
1067 intRF.attachToRecord( r, "ARRAY_ID" ) ;
1068 *intRF = 0 ;
1069 intRF.attachToRecord( r, "PROCESSOR_ID" ) ;
1070 *intRF = 0 ;
1071 RecordFieldPtr< Vector<Double> > arrayRF( r, "UVW" ) ;
1072 arrayRF.define( Vector<Double>( 3, 0.0 ) ) ;
1073 }
1074 void attachPointing()
1075 {
1076 porow = TableRow( potab ) ;
1077 TableRecord &r = porow.record() ;
1078 poNumPolyRF.attachToRecord( r, "NUM_POLY" ) ;
1079 poTimeRF.attachToRecord( r, "TIME" ) ;
1080 poTimeOriginRF.attachToRecord( r, "TIME_ORIGIN" ) ;
1081 poIntervalRF.attachToRecord( r, "INTERVAL" ) ;
1082 poNameRF.attachToRecord( r, "NAME" ) ;
1083 poDirectionRF.attachToRecord( r, "DIRECTION" ) ;
1084 poTargetRF.attachToRecord( r, "TARGET" ) ;
1085
1086 // constant values
1087 RecordFieldPtr<Int> antIdRF( r, "ANTENNA_ID" ) ;
1088 *antIdRF = 0 ;
1089 RecordFieldPtr<Bool> trackingRF( r, "TRACKING" ) ;
1090 *trackingRF = True ;
1091 }
1092 void queryType( Int type, String &stype, Bool &b, Double &t, Double &l )
1093 {
1094 t = 0.0 ;
1095 l = 0.0 ;
1096
1097 String sep1="#" ;
1098 String sep2="," ;
1099 String target="OBSERVE_TARGET" ;
1100 String atmcal="CALIBRATE_TEMPERATURE" ;
1101 String onstr="ON_SOURCE" ;
1102 String offstr="OFF_SOURCE" ;
1103 String pswitch="POSITION_SWITCH" ;
1104 String nod="NOD" ;
1105 String fswitch="FREQUENCY_SWITCH" ;
1106 String sigstr="SIG" ;
1107 String refstr="REF" ;
1108 String unspecified="UNSPECIFIED" ;
1109 String ftlow="LOWER" ;
1110 String fthigh="HIGHER" ;
1111 switch ( type ) {
1112 case SrcType::PSON:
1113 stype = target+sep1+onstr+sep2+pswitch ;
1114 b = True ;
1115 break ;
1116 case SrcType::PSOFF:
1117 stype = target+sep1+offstr+sep2+pswitch ;
1118 b = False ;
1119 break ;
1120 case SrcType::NOD:
1121 stype = target+sep1+onstr+sep2+nod ;
1122 b = True ;
1123 break ;
1124 case SrcType::FSON:
1125 stype = target+sep1+onstr+sep2+fswitch+sep1+sigstr ;
1126 b = True ;
1127 break ;
1128 case SrcType::FSOFF:
1129 stype = target+sep1+onstr+sep2+fswitch+sep1+refstr ;
1130 b = False ;
1131 break ;
1132 case SrcType::SKY:
1133 stype = atmcal+sep1+offstr+sep2+unspecified ;
1134 b = False ;
1135 break ;
1136 case SrcType::HOT:
1137 stype = atmcal+sep1+offstr+sep2+unspecified ;
1138 b = False ;
1139 break ;
1140 case SrcType::WARM:
1141 stype = atmcal+sep1+offstr+sep2+unspecified ;
1142 b = False ;
1143 break ;
1144 case SrcType::COLD:
1145 stype = atmcal+sep1+offstr+sep2+unspecified ;
1146 b = False ;
1147 break ;
1148 case SrcType::PONCAL:
1149 stype = atmcal+sep1+onstr+sep2+pswitch ;
1150 b = True ;
1151 break ;
1152 case SrcType::POFFCAL:
1153 stype = atmcal+sep1+offstr+sep2+pswitch ;
1154 b = False ;
1155 break ;
1156 case SrcType::NODCAL:
1157 stype = atmcal+sep1+onstr+sep2+nod ;
1158 b = True ;
1159 break ;
1160 case SrcType::FONCAL:
1161 stype = atmcal+sep1+onstr+sep2+fswitch+sep1+sigstr ;
1162 b = True ;
1163 break ;
1164 case SrcType::FOFFCAL:
1165 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+refstr ;
1166 b = False ;
1167 break ;
1168 case SrcType::FSLO:
1169 stype = target+sep1+onstr+sep2+fswitch+sep1+ftlow ;
1170 b = True ;
1171 break ;
1172 case SrcType::FLOOFF:
1173 stype = target+sep1+offstr+sep2+fswitch+sep1+ftlow ;
1174 b = False ;
1175 break ;
1176 case SrcType::FLOSKY:
1177 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
1178 b = False ;
1179 break ;
1180 case SrcType::FLOHOT:
1181 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
1182 b = False ;
1183 break ;
1184 case SrcType::FLOWARM:
1185 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
1186 b = False ;
1187 break ;
1188 case SrcType::FLOCOLD:
1189 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
1190 b = False ;
1191 break ;
1192 case SrcType::FSHI:
1193 stype = target+sep1+onstr+sep2+fswitch+sep1+fthigh ;
1194 b = True ;
1195 break ;
1196 case SrcType::FHIOFF:
1197 stype = target+sep1+offstr+sep2+fswitch+sep1+fthigh ;
1198 b = False ;
1199 break ;
1200 case SrcType::FHISKY:
1201 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
1202 b = False ;
1203 break ;
1204 case SrcType::FHIHOT:
1205 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
1206 b = False ;
1207 break ;
1208 case SrcType::FHIWARM:
1209 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
1210 b = False ;
1211 break ;
1212 case SrcType::FHICOLD:
1213 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
1214 b = False ;
1215 break ;
1216 case SrcType::SIG:
1217 stype = target+sep1+onstr+sep2+unspecified ;
1218 b = True ;
1219 break ;
1220 case SrcType::REF:
1221 stype = target+sep1+offstr+sep2+unspecified ;
1222 b = False ;
1223 break ;
1224 default:
1225 stype = unspecified ;
1226 b = True ;
1227 break ;
1228 }
1229 }
1230 void polProperty( Vector<uInt> &nos, Int &n, Vector<Int> &c, Matrix<Int> &cp )
1231 {
1232 n = nos.nelements() ;
1233 c.resize( n ) ;
1234 for ( Int i = 0 ; i < n ; i++ )
1235 c[i] = (Int)polmap[nos[i]] ;
1236 cp.resize( 2, n ) ;
1237 if ( n == 1 )
1238 cp = 0 ;
1239 else if ( n == 2 ) {
1240 cp.column( 0 ) = 0 ;
1241 cp.column( 1 ) = 1 ;
1242 }
1243 else {
1244 cp.column( 0 ) = 0 ;
1245 cp.column( 1 ) = 1 ;
1246 cp( 0, 1 ) = 0 ;
1247 cp( 1, 1 ) = 1 ;
1248 cp( 0, 2 ) = 1 ;
1249 cp( 1, 2 ) = 0 ;
1250 }
1251 }
1252 void initTcal()
1253 {
1254 const TableRecord &rec = table.keywordSet() ;
1255 Table tcalTable = rec.asTable( "TCAL" ) ;
1256 ROScalarColumn<uInt> idCol( tcalTable, "ID" ) ;
1257 Vector<uInt> id = idCol.getColumn() ;
1258 uInt maxId = max( id ) ;
1259 tcalNotYet.resize( maxId+1 ) ;
1260 tcalNotYet = True ;
1261 tcalKey = -1 ;
1262 }
1263
1264 Table &ms;
1265 TableRow row;
1266 uInt rowidx;
1267 String fieldName;
1268 Int fieldId;
1269 Int srcId;
1270 Int defaultFieldId;
1271 Int spwId;
1272 Int feedId;
1273 Int subscan;
1274 PolarizedComponentHolder holder;
1275 String ptName;
1276 Bool useFloat;
1277 String poltype;
1278 Vector<Stokes::StokesTypes> polmap;
1279
1280 // MS subtables
1281 Table spwtab;
1282 Table statetab;
1283 Table ddtab;
1284 Table poltab;
1285 Table fieldtab;
1286 Table feedtab;
1287 Table potab;
1288
1289 // Scantable MAIN columns
1290 ROArrayColumn<Float> spectraCol;
1291 ROArrayColumn<Double> directionCol,scanRateCol,sourceDirectionCol;
1292 ROArrayColumn<uChar> flagtraCol;
1293 ROTableColumn tcalIdCol,intervalCol,flagRowCol,timeCol,freqIdCol,
1294 sourceNameCol,fieldNameCol;
1295
1296 // MS MAIN columns
1297 RecordFieldPtr<Int> dataDescIdRF,fieldIdRF,feed1RF,feed2RF,
1298 scanNumberRF,stateIdRF;
1299 RecordFieldPtr<Bool> flagRowRF;
1300 RecordFieldPtr<Double> timeRF,timeCentroidRF,intervalRF,exposureRF;
1301 RecordFieldPtr< Vector<Float> > weightRF,sigmaRF;
1302 RecordFieldPtr< Cube<Bool> > flagCategoryRF;
1303 RecordFieldPtr< Matrix<Bool> > flagRF;
1304 RecordFieldPtr< Matrix<Float> > floatDataRF;
1305 RecordFieldPtr< Matrix<Complex> > dataRF;
1306
1307 // MS POINTING columns
1308 TableRow porow;
1309 RecordFieldPtr<Int> poNumPolyRF ;
1310 RecordFieldPtr<Double> poTimeRF,
1311 poTimeOriginRF,
1312 poIntervalRF ;
1313 RecordFieldPtr<String> poNameRF ;
1314 RecordFieldPtr< Matrix<Double> > poDirectionRF,
1315 poTargetRF ;
1316
1317 Vector<String> stateEntry;
1318 Matrix<Int> ddEntry;
1319 Matrix<Int> feedEntry;
1320 vector< Vector<uInt> > polEntry;
1321 map<uInt,Bool> processedFreqId;
1322 map<uInt,Double> refpix;
1323 map<uInt,Double> refval;
1324 map<uInt,Double> increment;
1325 MFrequency::Types freqframe;
1326 Record srcRec;
1327 map< Int,Vector<uInt> > tcalIdRec;
1328 map< Int,Vector<uInt> > tcalRowRec;
1329 Int tcalKey;
1330 Vector<Bool> tcalNotYet;
1331};
1332
1333MSWriter::MSWriter(CountedPtr<Scantable> stable)
1334 : table_(stable),
1335 isTcal_(False),
1336 isWeather_(False),
1337 tcalSpec_(False),
1338 tsysSpec_(False),
1339 ptTabName_("")
1340{
1341 os_ = LogIO() ;
1342 os_.origin( LogOrigin( "MSWriter", "MSWriter()", WHERE ) ) ;
1343// os_ << "MSWriter::MSWriter()" << LogIO::POST ;
1344
1345 // initialize writer
1346 init() ;
1347}
1348
1349MSWriter::~MSWriter()
1350{
1351 os_.origin( LogOrigin( "MSWriter", "~MSWriter()", WHERE ) ) ;
1352// os_ << "MSWriter::~MSWriter()" << LogIO::POST ;
1353
1354 if ( mstable_ != 0 )
1355 delete mstable_ ;
1356}
1357
1358bool MSWriter::write(const string& filename, const Record& rec)
1359{
1360 os_.origin( LogOrigin( "MSWriter", "write()", WHERE ) ) ;
1361 //double startSec = mathutil::gettimeofday_sec() ;
1362 //os_ << "start MSWriter::write() startSec=" << startSec << LogIO::POST ;
1363
1364 filename_ = filename ;
1365
1366 // parsing MS options
1367 Bool overwrite = False ;
1368 if ( rec.isDefined( "ms" ) ) {
1369 Record msrec = rec.asRecord( "ms" ) ;
1370 if ( msrec.isDefined( "overwrite" ) ) {
1371 overwrite = msrec.asBool( "overwrite" ) ;
1372 }
1373 }
1374
1375 os_ << "Parsing MS options" << endl ;
1376 os_ << " overwrite = " << overwrite << LogIO::POST ;
1377
1378 File file( filename_ ) ;
1379 if ( file.exists() ) {
1380 if ( overwrite ) {
1381 os_ << filename_ << " exists. Overwrite existing data... " << LogIO::POST ;
1382 if ( file.isRegular() ) RegularFile(file).remove() ;
1383 else if ( file.isDirectory() ) Directory(file).removeRecursive() ;
1384 else SymLink(file).remove() ;
1385 }
1386 else {
1387 os_ << LogIO::SEVERE << "ERROR: " << filename_ << " exists..." << LogIO::POST ;
1388 return False ;
1389 }
1390 }
1391
1392 // set up MS
1393 setupMS() ;
1394
1395 // subtables
1396 // OBSERVATION
1397 fillObservation() ;
1398
1399 // ANTENNA
1400 fillAntenna() ;
1401
1402 // PROCESSOR
1403 fillProcessor() ;
1404
1405 // SOURCE
1406 fillSource() ;
1407
1408 // WEATHER
1409 if ( isWeather_ )
1410 fillWeather() ;
1411
1412 /***
1413 * Start iteration using TableVisitor
1414 ***/
1415 {
1416 static const char *cols[] = {
1417 "FIELDNAME", "BEAMNO", "SCANNO", "IFNO", "SRCTYPE", "CYCLENO", "TIME",
1418 "POLNO",
1419 NULL
1420 };
1421 static const TypeManagerImpl<uInt> tmUInt;
1422 static const TypeManagerImpl<Int> tmInt;
1423 static const TypeManagerImpl<Double> tmDouble;
1424 static const TypeManagerImpl<String> tmString;
1425 static const TypeManager *const tms[] = {
1426 &tmString, &tmUInt, &tmUInt, &tmUInt, &tmInt, &tmUInt, &tmDouble, &tmInt, NULL
1427 };
1428 //double t0 = mathutil::gettimeofday_sec() ;
1429 MSWriterVisitor myVisitor(table_->table(),*mstable_);
1430 //double t1 = mathutil::gettimeofday_sec() ;
1431 //cout << "MSWriterVisitor(): elapsed time " << t1-t0 << " sec" << endl ;
1432 String dataColName = "FLOAT_DATA" ;
1433 if ( useData_ )
1434 dataColName = "DATA" ;
1435 myVisitor.dataColumnName( dataColName ) ;
1436 myVisitor.pointingTableName( ptTabName_ ) ;
1437 myVisitor.setSourceRecord( srcRec_ ) ;
1438 //double t2 = mathutil::gettimeofday_sec() ;
1439 traverseTable(table_->table(), cols, tms, &myVisitor);
1440 //double t3 = mathutil::gettimeofday_sec() ;
1441 //cout << "traverseTable(): elapsed time " << t3-t2 << " sec" << endl ;
1442 map< Int,Vector<uInt> > &idRec = myVisitor.getTcalIdRecord() ;
1443 map< Int,Vector<uInt> > &rowRec = myVisitor.getTcalRowRecord() ;
1444
1445 // SYSCAL
1446 if ( isTcal_ )
1447 fillSysCal( idRec, rowRec ) ;
1448 }
1449 /***
1450 * End iteration using TableVisitor
1451 ***/
1452
1453 // ASDM tables
1454 const TableRecord &stKeys = table_->table().keywordSet() ;
1455 TableRecord &msKeys = mstable_->rwKeywordSet() ;
1456 uInt nfields = stKeys.nfields() ;
1457 for ( uInt ifield = 0 ; ifield < nfields ; ifield++ ) {
1458 String kname = stKeys.name( ifield ) ;
1459 if ( kname.find( "ASDM" ) != String::npos ) {
1460 String asdmpath = stKeys.asString( ifield ) ;
1461 os_ << "found ASDM table: " << asdmpath << LogIO::POST ;
1462 if ( Table::isReadable( asdmpath ) ) {
1463 Table newAsdmTab( asdmpath, Table::Old ) ;
1464 newAsdmTab.copy( filename_+"/"+kname, Table::New ) ;
1465 os_ << "add subtable: " << kname << LogIO::POST ;
1466 msKeys.defineTable( kname, Table( filename_+"/"+kname, Table::Old ) ) ;
1467 }
1468 }
1469 }
1470
1471 // replace POINTING table with original one if exists
1472 if ( ptTabName_ != "" ) {
1473 delete mstable_ ;
1474 mstable_ = 0 ;
1475 Table newPtTab( ptTabName_, Table::Old ) ;
1476 newPtTab.copy( filename_+"/POINTING", Table::New ) ;
1477 }
1478
1479 //double endSec = mathutil::gettimeofday_sec() ;
1480 //os_ << "end MSWriter::write() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1481
1482 return True ;
1483}
1484
1485void MSWriter::init()
1486{
1487// os_.origin( LogOrigin( "MSWriter", "init()", WHERE ) ) ;
1488// double startSec = mathutil::gettimeofday_sec() ;
1489// os_ << "start MSWriter::init() startSec=" << startSec << LogIO::POST ;
1490
1491 // access to scantable
1492 header_ = table_->getHeader() ;
1493
1494 // FLOAT_DATA? or DATA?
1495 if ( header_.npol > 2 ) {
1496 useFloatData_ = False ;
1497 useData_ = True ;
1498 }
1499 else {
1500 useFloatData_ = True ;
1501 useData_ = False ;
1502 }
1503
1504 // polarization type
1505 polType_ = header_.poltype ;
1506 if ( polType_ == "" )
1507 polType_ = "stokes" ;
1508 else if ( polType_.find( "linear" ) != String::npos )
1509 polType_ = "linear" ;
1510 else if ( polType_.find( "circular" ) != String::npos )
1511 polType_ = "circular" ;
1512 else if ( polType_.find( "stokes" ) != String::npos )
1513 polType_ = "stokes" ;
1514 else if ( polType_.find( "linpol" ) != String::npos )
1515 polType_ = "linpol" ;
1516 else
1517 polType_ = "notype" ;
1518
1519 // Check if some subtables are exists
1520 if ( table_->tcal().table().nrow() != 0 ) {
1521 ROTableColumn col( table_->tcal().table(), "TCAL" ) ;
1522 if ( col.isDefined( 0 ) ) {
1523 os_ << "TCAL table exists: nrow=" << table_->tcal().table().nrow() << LogIO::POST ;
1524 isTcal_ = True ;
1525 }
1526 else {
1527 os_ << "No TCAL rows" << LogIO::POST ;
1528 }
1529 }
1530 else {
1531 os_ << "No TCAL rows" << LogIO::POST ;
1532 }
1533 if ( table_->weather().table().nrow() != 0 ) {
1534 ROTableColumn col( table_->weather().table(), "TEMPERATURE" ) ;
1535 if ( col.isDefined( 0 ) ) {
1536 os_ << "WEATHER table exists: nrow=" << table_->weather().table().nrow() << LogIO::POST ;
1537 isWeather_ =True ;
1538 }
1539 else {
1540 os_ << "No WEATHER rows" << LogIO::POST ;
1541 }
1542 }
1543 else {
1544 os_ << "No WEATHER rows" << LogIO::POST ;
1545 }
1546
1547 // Are TCAL_SPECTRUM and TSYS_SPECTRUM necessary?
1548 if ( isTcal_ && header_.nchan != 1 ) {
1549 // examine TCAL subtable
1550 Table tcaltab = table_->tcal().table() ;
1551 ROArrayColumn<Float> tcalCol( tcaltab, "TCAL" ) ;
1552 for ( uInt irow = 0 ; irow < tcaltab.nrow() ; irow++ ) {
1553 if ( tcalCol( irow ).size() != 1 )
1554 tcalSpec_ = True ;
1555 }
1556 // examine spectral data
1557 TableIterator iter0( table_->table(), "IFNO" ) ;
1558 while( !iter0.pastEnd() ) {
1559 Table t0( iter0.table() ) ;
1560 ROArrayColumn<Float> sharedFloatArrCol( t0, "SPECTRA" ) ;
1561 uInt len = sharedFloatArrCol( 0 ).size() ;
1562 if ( len != 1 ) {
1563 sharedFloatArrCol.attach( t0, "TSYS" ) ;
1564 if ( sharedFloatArrCol( 0 ).size() != 1 )
1565 tsysSpec_ = True ;
1566 }
1567 iter0.next() ;
1568 }
1569 }
1570
1571 // check if reference for POINTING table exists
1572 const TableRecord &rec = table_->table().keywordSet() ;
1573 if ( rec.isDefined( "POINTING" ) ) {
1574 ptTabName_ = rec.asString( "POINTING" ) ;
1575 if ( !Table::isReadable( ptTabName_ ) ) {
1576 ptTabName_ = "" ;
1577 }
1578 }
1579
1580// double endSec = mathutil::gettimeofday_sec() ;
1581// os_ << "end MSWriter::init() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1582}
1583
1584void MSWriter::setupMS()
1585{
1586// os_.origin( LogOrigin( "MSWriter", "setupMS()", WHERE ) ) ;
1587// double startSec = mathutil::gettimeofday_sec() ;
1588// os_ << "start MSWriter::setupMS() startSec=" << startSec << LogIO::POST ;
1589
1590 String dunit = table_->getHeader().fluxunit ;
1591
1592 TableDesc msDesc = MeasurementSet::requiredTableDesc() ;
1593 if ( useFloatData_ )
1594 MeasurementSet::addColumnToDesc( msDesc, MSMainEnums::FLOAT_DATA, 2 ) ;
1595 else if ( useData_ )
1596 MeasurementSet::addColumnToDesc( msDesc, MSMainEnums::DATA, 2 ) ;
1597
1598 SetupNewTable newtab( filename_, msDesc, Table::New ) ;
1599
1600 mstable_ = new MeasurementSet( newtab ) ;
1601
1602 TableColumn col ;
1603 if ( useFloatData_ )
1604 col.attach( *mstable_, "FLOAT_DATA" ) ;
1605 else if ( useData_ )
1606 col.attach( *mstable_, "DATA" ) ;
1607 col.rwKeywordSet().define( "UNIT", dunit ) ;
1608
1609 // create subtables
1610 TableDesc antennaDesc = MSAntenna::requiredTableDesc() ;
1611 SetupNewTable antennaTab( mstable_->antennaTableName(), antennaDesc, Table::New ) ;
1612 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::ANTENNA ), Table( antennaTab ) ) ;
1613
1614 TableDesc dataDescDesc = MSDataDescription::requiredTableDesc() ;
1615 SetupNewTable dataDescTab( mstable_->dataDescriptionTableName(), dataDescDesc, Table::New ) ;
1616 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::DATA_DESCRIPTION ), Table( dataDescTab ) ) ;
1617
1618 TableDesc dopplerDesc = MSDoppler::requiredTableDesc() ;
1619 SetupNewTable dopplerTab( mstable_->dopplerTableName(), dopplerDesc, Table::New ) ;
1620 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::DOPPLER ), Table( dopplerTab ) ) ;
1621
1622 TableDesc feedDesc = MSFeed::requiredTableDesc() ;
1623 SetupNewTable feedTab( mstable_->feedTableName(), feedDesc, Table::New ) ;
1624 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FEED ), Table( feedTab ) ) ;
1625
1626 TableDesc fieldDesc = MSField::requiredTableDesc() ;
1627 SetupNewTable fieldTab( mstable_->fieldTableName(), fieldDesc, Table::New ) ;
1628 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FIELD ), Table( fieldTab ) ) ;
1629
1630 TableDesc flagCmdDesc = MSFlagCmd::requiredTableDesc() ;
1631 SetupNewTable flagCmdTab( mstable_->flagCmdTableName(), flagCmdDesc, Table::New ) ;
1632 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FLAG_CMD ), Table( flagCmdTab ) ) ;
1633
1634 TableDesc freqOffsetDesc = MSFreqOffset::requiredTableDesc() ;
1635 SetupNewTable freqOffsetTab( mstable_->freqOffsetTableName(), freqOffsetDesc, Table::New ) ;
1636 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FREQ_OFFSET ), Table( freqOffsetTab ) ) ;
1637
1638 TableDesc historyDesc = MSHistory::requiredTableDesc() ;
1639 SetupNewTable historyTab( mstable_->historyTableName(), historyDesc, Table::New ) ;
1640 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::HISTORY ), Table( historyTab ) ) ;
1641
1642 TableDesc observationDesc = MSObservation::requiredTableDesc() ;
1643 SetupNewTable observationTab( mstable_->observationTableName(), observationDesc, Table::New ) ;
1644 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::OBSERVATION ), Table( observationTab ) ) ;
1645
1646 TableDesc pointingDesc = MSPointing::requiredTableDesc() ;
1647 SetupNewTable pointingTab( mstable_->pointingTableName(), pointingDesc, Table::New ) ;
1648 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::POINTING ), Table( pointingTab ) ) ;
1649
1650 TableDesc polarizationDesc = MSPolarization::requiredTableDesc() ;
1651 SetupNewTable polarizationTab( mstable_->polarizationTableName(), polarizationDesc, Table::New ) ;
1652 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::POLARIZATION ), Table( polarizationTab ) ) ;
1653
1654 TableDesc processorDesc = MSProcessor::requiredTableDesc() ;
1655 SetupNewTable processorTab( mstable_->processorTableName(), processorDesc, Table::New ) ;
1656 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::PROCESSOR ), Table( processorTab ) ) ;
1657
1658 TableDesc sourceDesc = MSSource::requiredTableDesc() ;
1659 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::TRANSITION, 1 ) ;
1660 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::REST_FREQUENCY, 1 ) ;
1661 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::SYSVEL, 1 ) ;
1662 SetupNewTable sourceTab( mstable_->sourceTableName(), sourceDesc, Table::New ) ;
1663 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SOURCE ), Table( sourceTab ) ) ;
1664
1665 TableDesc spwDesc = MSSpectralWindow::requiredTableDesc() ;
1666 SetupNewTable spwTab( mstable_->spectralWindowTableName(), spwDesc, Table::New ) ;
1667 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SPECTRAL_WINDOW ), Table( spwTab ) ) ;
1668
1669 TableDesc stateDesc = MSState::requiredTableDesc() ;
1670 SetupNewTable stateTab( mstable_->stateTableName(), stateDesc, Table::New ) ;
1671 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::STATE ), Table( stateTab ) ) ;
1672
1673 // TODO: add TCAL_SPECTRUM and TSYS_SPECTRUM if necessary
1674 TableDesc sysCalDesc = MSSysCal::requiredTableDesc() ;
1675 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TCAL, 2 ) ;
1676 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TSYS, 2 ) ;
1677 if ( tcalSpec_ )
1678 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TCAL_SPECTRUM, 2 ) ;
1679 if ( tsysSpec_ )
1680 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TSYS_SPECTRUM, 2 ) ;
1681 SetupNewTable sysCalTab( mstable_->sysCalTableName(), sysCalDesc, Table::New ) ;
1682 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SYSCAL ), Table( sysCalTab ) ) ;
1683
1684 TableDesc weatherDesc = MSWeather::requiredTableDesc() ;
1685 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::TEMPERATURE ) ;
1686 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::PRESSURE ) ;
1687 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::REL_HUMIDITY ) ;
1688 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::WIND_SPEED ) ;
1689 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::WIND_DIRECTION ) ;
1690 SetupNewTable weatherTab( mstable_->weatherTableName(), weatherDesc, Table::New ) ;
1691 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::WEATHER ), Table( weatherTab ) ) ;
1692
1693 mstable_->initRefs() ;
1694
1695// double endSec = mathutil::gettimeofday_sec() ;
1696// os_ << "end MSWriter::setupMS() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1697}
1698
1699void MSWriter::fillObservation()
1700{
1701 //double startSec = mathutil::gettimeofday_sec() ;
1702 //os_ << "start MSWriter::fillObservation() startSec=" << startSec << LogIO::POST ;
1703
1704 // only 1 row
1705 mstable_->observation().addRow( 1, True ) ;
1706 MSObservationColumns msObsCols( mstable_->observation() ) ;
1707 msObsCols.observer().put( 0, header_.observer ) ;
1708 // tentatively put antennaname (from ANTENNA subtable)
1709 String hAntennaName = header_.antennaname ;
1710 String::size_type pos = hAntennaName.find( "//" ) ;
1711 String telescopeName ;
1712 if ( pos != String::npos ) {
1713 telescopeName = hAntennaName.substr( 0, pos ) ;
1714 }
1715 else {
1716 pos = hAntennaName.find( "@" ) ;
1717 telescopeName = hAntennaName.substr( 0, pos ) ;
1718 }
1719// os_ << "telescopeName = " << telescopeName << LogIO::POST ;
1720 msObsCols.telescopeName().put( 0, telescopeName ) ;
1721 msObsCols.project().put( 0, header_.project ) ;
1722 //ScalarMeasColumn<MEpoch> timeCol( table_->table().sort("TIME"), "TIME" ) ;
1723 Table sortedtable = table_->table().sort("TIME") ;
1724 ScalarMeasColumn<MEpoch> timeCol( sortedtable, "TIME" ) ;
1725 Vector<MEpoch> trange( 2 ) ;
1726 trange[0] = timeCol( 0 ) ;
1727 trange[1] = timeCol( table_->nrow()-1 ) ;
1728 msObsCols.timeRangeMeas().put( 0, trange ) ;
1729
1730 //double endSec = mathutil::gettimeofday_sec() ;
1731 //os_ << "end MSWriter::fillObservation() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1732}
1733
1734void MSWriter::antennaProperty( String &name, String &m, String &t, Double &d )
1735{
1736 name.upcase() ;
1737
1738 m = "ALT-AZ" ;
1739 t = "GROUND-BASED" ;
1740 if ( name.matches( Regex( "DV[0-9]+$" ) )
1741 || name.matches( Regex( "DA[0-9]+$" ) )
1742 || name.matches( Regex( "PM[0-9]+$" ) ) )
1743 d = 12.0 ;
1744 else if ( name.matches( Regex( "CM[0-9]+$" ) ) )
1745 d = 7.0 ;
1746 else if ( name.contains( "GBT" ) )
1747 d = 104.9 ;
1748 else if ( name.contains( "MOPRA" ) )
1749 d = 22.0 ;
1750 else if ( name.contains( "PKS" ) || name.contains( "PARKS" ) )
1751 d = 64.0 ;
1752 else if ( name.contains( "TIDBINBILLA" ) )
1753 d = 70.0 ;
1754 else if ( name.contains( "CEDUNA" ) )
1755 d = 30.0 ;
1756 else if ( name.contains( "HOBART" ) )
1757 d = 26.0 ;
1758 else if ( name.contains( "APEX" ) )
1759 d = 12.0 ;
1760 else if ( name.contains( "ASTE" ) )
1761 d = 10.0 ;
1762 else if ( name.contains( "NRO" ) )
1763 d = 45.0 ;
1764 else
1765 d = 1.0 ;
1766}
1767
1768void MSWriter::fillAntenna()
1769{
1770 //double startSec = mathutil::gettimeofday_sec() ;
1771 //os_ << "start MSWriter::fillAntenna() startSec=" << startSec << LogIO::POST ;
1772
1773 // only 1 row
1774 Table anttab = mstable_->antenna() ;
1775 anttab.addRow( 1, True ) ;
1776
1777 Table &table = table_->table() ;
1778 const TableRecord &keys = table.keywordSet() ;
1779 String hAntName = keys.asString( "AntennaName" ) ;
1780 String::size_type pos = hAntName.find( "//" ) ;
1781 String antennaName ;
1782 String stationName ;
1783 if ( pos != String::npos ) {
1784 stationName = hAntName.substr( 0, pos ) ;
1785 hAntName = hAntName.substr( pos+2 ) ;
1786 }
1787 pos = hAntName.find( "@" ) ;
1788 if ( pos != String::npos ) {
1789 antennaName = hAntName.substr( 0, pos ) ;
1790 stationName = hAntName.substr( pos+1 ) ;
1791 }
1792 else {
1793 antennaName = hAntName ;
1794 }
1795 Vector<Double> antpos = keys.asArrayDouble( "AntennaPosition" ) ;
1796
1797 String mount, atype ;
1798 Double diameter ;
1799 antennaProperty( antennaName, mount, atype, diameter ) ;
1800
1801 TableRow tr( anttab ) ;
1802 TableRecord &r = tr.record() ;
1803 RecordFieldPtr<String> nameRF( r, "NAME" ) ;
1804 RecordFieldPtr<String> stationRF( r, "STATION" ) ;
1805 RecordFieldPtr<String> mountRF( r, "NAME" ) ;
1806 RecordFieldPtr<String> typeRF( r, "TYPE" ) ;
1807 RecordFieldPtr<Double> dishDiameterRF( r, "DISH_DIAMETER" ) ;
1808 RecordFieldPtr< Vector<Double> > positionRF( r, "POSITION" ) ;
1809 *nameRF = antennaName ;
1810 *mountRF = mount ;
1811 *typeRF = atype ;
1812 *dishDiameterRF = diameter ;
1813 *positionRF = antpos ;
1814
1815 tr.put( 0 ) ;
1816
1817 //double endSec = mathutil::gettimeofday_sec() ;
1818 //os_ << "end MSWriter::fillAntenna() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1819}
1820
1821// void MSWriter::fillAntenna()
1822// {
1823// // double startSec = mathutil::gettimeofday_sec() ;
1824// // os_ << "start MSWriter::fillAntenna() startSec=" << startSec << LogIO::POST ;
1825
1826// // only 1 row
1827// mstable_->antenna().addRow( 1, True ) ;
1828// MSAntennaColumns msAntCols( mstable_->antenna() ) ;
1829
1830// String hAntennaName = header_.antennaname ;
1831// String::size_type pos = hAntennaName.find( "//" ) ;
1832// String antennaName ;
1833// String stationName ;
1834// if ( pos != String::npos ) {
1835// hAntennaName = hAntennaName.substr( pos+2 ) ;
1836// }
1837// pos = hAntennaName.find( "@" ) ;
1838// if ( pos != String::npos ) {
1839// antennaName = hAntennaName.substr( 0, pos ) ;
1840// stationName = hAntennaName.substr( pos+1 ) ;
1841// }
1842// else {
1843// antennaName = hAntennaName ;
1844// stationName = hAntennaName ;
1845// }
1846// // os_ << "antennaName = " << antennaName << LogIO::POST ;
1847// // os_ << "stationName = " << stationName << LogIO::POST ;
1848
1849// msAntCols.name().put( 0, antennaName ) ;
1850// msAntCols.station().put( 0, stationName ) ;
1851
1852// // os_ << "antennaPosition = " << header_.antennaposition << LogIO::POST ;
1853
1854// msAntCols.position().put( 0, header_.antennaposition ) ;
1855
1856// // MOUNT is set to "ALT-AZ"
1857// msAntCols.mount().put( 0, "ALT-AZ" ) ;
1858
1859// Double diameter = getDishDiameter( antennaName ) ;
1860// msAntCols.dishDiameterQuant().put( 0, Quantity( diameter, "m" ) ) ;
1861
1862// // double endSec = mathutil::gettimeofday_sec() ;
1863// // os_ << "end MSWriter::fillAntenna() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1864// }
1865
1866void MSWriter::fillProcessor()
1867{
1868// double startSec = mathutil::gettimeofday_sec() ;
1869// os_ << "start MSWriter::fillProcessor() startSec=" << startSec << LogIO::POST ;
1870
1871 // only add empty 1 row
1872 MSProcessor msProc = mstable_->processor() ;
1873 msProc.addRow( 1, True ) ;
1874
1875// double endSec = mathutil::gettimeofday_sec() ;
1876// os_ << "end MSWriter::fillProcessor() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1877}
1878
1879void MSWriter::fillSource()
1880{
1881// double startSec = mathutil::gettimeofday_sec() ;
1882// os_ << "start MSWriter::fillSource() startSec=" << startSec << LogIO::POST ;
1883
1884 // access to MS SOURCE subtable
1885 MSSource msSrc = mstable_->source() ;
1886
1887 // access to MOLECULE subtable
1888 STMolecules stm = table_->molecules() ;
1889
1890 Int srcId = 0 ;
1891 Vector<Double> restFreq ;
1892 Vector<String> molName ;
1893 Vector<String> fMolName ;
1894
1895 // row based
1896 TableRow row( msSrc ) ;
1897 TableRecord &rec = row.record() ;
1898 RecordFieldPtr<Int> srcidRF( rec, "SOURCE_ID" ) ;
1899 RecordFieldPtr<String> nameRF( rec, "NAME" ) ;
1900 RecordFieldPtr< Array<Double> > srcpmRF( rec, "PROPER_MOTION" ) ;
1901 RecordFieldPtr< Array<Double> > srcdirRF( rec, "DIRECTION" ) ;
1902 RecordFieldPtr<Int> numlineRF( rec, "NUM_LINES" ) ;
1903 RecordFieldPtr< Array<Double> > restfreqRF( rec, "REST_FREQUENCY" ) ;
1904 RecordFieldPtr< Array<Double> > sysvelRF( rec, "SYSVEL" ) ;
1905 RecordFieldPtr< Array<String> > transitionRF( rec, "TRANSITION" ) ;
1906 RecordFieldPtr<Double> timeRF( rec, "TIME" ) ;
1907 RecordFieldPtr<Double> intervalRF( rec, "INTERVAL" ) ;
1908 RecordFieldPtr<Int> spwidRF( rec, "SPECTRAL_WINDOW_ID" ) ;
1909
1910 //
1911 // ITERATION: SRCNAME
1912 //
1913 TableIterator iter0( table_->table(), "SRCNAME" ) ;
1914 while( !iter0.pastEnd() ) {
1915 //Table t0( iter0.table() ) ;
1916 Table t0 = iter0.table() ;
1917
1918 // get necessary information
1919 ROScalarColumn<String> srcNameCol( t0, "SRCNAME" ) ;
1920 String srcName = srcNameCol( 0 ) ;
1921 ROArrayColumn<Double> sharedDArrRCol( t0, "SRCPROPERMOTION" ) ;
1922 Vector<Double> srcPM = sharedDArrRCol( 0 ) ;
1923 sharedDArrRCol.attach( t0, "SRCDIRECTION" ) ;
1924 Vector<Double> srcDir = sharedDArrRCol( 0 ) ;
1925 ROScalarColumn<Double> srcVelCol( t0, "SRCVELOCITY" ) ;
1926 Double srcVel = srcVelCol( 0 ) ;
1927 srcRec_.define( srcName, srcId ) ;
1928
1929 // NAME
1930 *nameRF = srcName ;
1931
1932 // SOURCE_ID
1933 *srcidRF = srcId ;
1934
1935 // PROPER_MOTION
1936 *srcpmRF = srcPM ;
1937
1938 // DIRECTION
1939 *srcdirRF = srcDir ;
1940
1941 //
1942 // ITERATION: MOLECULE_ID
1943 //
1944 TableIterator iter1( t0, "MOLECULE_ID" ) ;
1945 while( !iter1.pastEnd() ) {
1946 //Table t1( iter1.table() ) ;
1947 Table t1 = iter1.table() ;
1948
1949 // get necessary information
1950 ROScalarColumn<uInt> molIdCol( t1, "MOLECULE_ID" ) ;
1951 uInt molId = molIdCol( 0 ) ;
1952 stm.getEntry( restFreq, molName, fMolName, molId ) ;
1953
1954 uInt numFreq = restFreq.size() ;
1955
1956 // NUM_LINES
1957 *numlineRF = numFreq ;
1958
1959 // REST_FREQUENCY
1960 *restfreqRF = restFreq ;
1961
1962 // TRANSITION
1963 //*transitionRF = fMolName ;
1964 Vector<String> transition ;
1965 if ( fMolName.size() != 0 ) {
1966 transition = fMolName ;
1967 }
1968 else if ( molName.size() != 0 ) {
1969 transition = molName ;
1970 }
1971 else {
1972 transition.resize( numFreq ) ;
1973 transition = "" ;
1974 }
1975 *transitionRF = transition ;
1976
1977 // SYSVEL
1978 Vector<Double> sysvelArr( numFreq, srcVel ) ;
1979 *sysvelRF = sysvelArr ;
1980
1981 //
1982 // ITERATION: IFNO
1983 //
1984 TableIterator iter2( t1, "IFNO" ) ;
1985 while( !iter2.pastEnd() ) {
1986 //Table t2( iter2.table() ) ;
1987 Table t2 = iter2.table() ;
1988 uInt nrow = msSrc.nrow() ;
1989
1990 // get necessary information
1991 ROScalarColumn<uInt> ifNoCol( t2, "IFNO" ) ;
1992 uInt ifno = ifNoCol( 0 ) ; // IFNO = SPECTRAL_WINDOW_ID
1993 Double midTime ;
1994 Double interval ;
1995 getValidTimeRange( midTime, interval, t2 ) ;
1996
1997 // fill SPECTRAL_WINDOW_ID
1998 *spwidRF = ifno ;
1999
2000 // fill TIME, INTERVAL
2001 *timeRF = midTime ;
2002 *intervalRF = interval ;
2003
2004 // add row
2005 msSrc.addRow( 1, True ) ;
2006 row.put( nrow ) ;
2007
2008 iter2.next() ;
2009 }
2010
2011 iter1.next() ;
2012 }
2013
2014 // increment srcId if SRCNAME changed
2015 srcId++ ;
2016
2017 iter0.next() ;
2018 }
2019
2020// double endSec = mathutil::gettimeofday_sec() ;
2021// os_ << "end MSWriter::fillSource() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2022}
2023
2024void MSWriter::fillWeather()
2025{
2026// double startSec = mathutil::gettimeofday_sec() ;
2027// os_ << "start MSWriter::fillWeather() startSec=" << startSec << LogIO::POST ;
2028
2029 // access to MS WEATHER subtable
2030 MSWeather msw = mstable_->weather() ;
2031
2032 // access to WEATHER subtable
2033 Table stw = table_->weather().table() ;
2034 uInt nrow = stw.nrow() ;
2035
2036 if ( nrow == 0 )
2037 return ;
2038
2039 msw.addRow( nrow, True ) ;
2040 MSWeatherColumns mswCols( msw ) ;
2041
2042 // ANTENNA_ID is always 0
2043 Vector<Int> antIdArr( nrow, 0 ) ;
2044 mswCols.antennaId().putColumn( antIdArr ) ;
2045
2046 // fill weather status
2047 ROScalarColumn<Float> sharedFloatCol( stw, "TEMPERATURE" ) ;
2048 mswCols.temperature().putColumn( sharedFloatCol ) ;
2049 sharedFloatCol.attach( stw, "PRESSURE" ) ;
2050 mswCols.pressure().putColumn( sharedFloatCol ) ;
2051 sharedFloatCol.attach( stw, "HUMIDITY" ) ;
2052 mswCols.relHumidity().putColumn( sharedFloatCol ) ;
2053 sharedFloatCol.attach( stw, "WINDSPEED" ) ;
2054 mswCols.windSpeed().putColumn( sharedFloatCol ) ;
2055 sharedFloatCol.attach( stw, "WINDAZ" ) ;
2056 mswCols.windDirection().putColumn( sharedFloatCol ) ;
2057
2058 // fill TIME and INTERVAL
2059 Double midTime ;
2060 Double interval ;
2061 Vector<Double> intervalArr( nrow, 0.0 ) ;
2062 TableIterator iter( table_->table(), "WEATHER_ID" ) ;
2063 while( !iter.pastEnd() ) {
2064 //Table tab( iter.table() ) ;
2065 Table tab = iter.table() ;
2066
2067 ROScalarColumn<uInt> widCol( tab, "WEATHER_ID" ) ;
2068 uInt wid = widCol( 0 ) ;
2069
2070 getValidTimeRange( midTime, interval, tab ) ;
2071 mswCols.time().put( wid, midTime ) ;
2072 intervalArr[wid] = interval ;
2073
2074 iter.next() ;
2075 }
2076 mswCols.interval().putColumn( intervalArr ) ;
2077
2078// double endSec = mathutil::gettimeofday_sec() ;
2079// os_ << "end MSWriter::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2080}
2081
2082void MSWriter::fillSysCal( map< Int,Vector<uInt> > &idRec,
2083 map< Int,Vector<uInt> > &rowRec )
2084{
2085 //double startSec = mathutil::gettimeofday_sec() ;
2086 //os_ << "start MSWriter::fillSysCal() startSec=" << startSec << LogIO::POST ;
2087
2088 //idRec.print( cout ) ;
2089
2090 // access to MS SYSCAL subtable
2091 MSSysCal mssc = mstable_->sysCal() ;
2092
2093 // access to TCAL subtable
2094 Table stt = table_->tcal().table() ;
2095 uInt nrow = stt.nrow() ;
2096
2097 // access to MAIN table
2098 Block<String> cols( 6 ) ;
2099 cols[0] = "TIME" ;
2100 cols[1] = "TCAL_ID" ;
2101 cols[2] = "TSYS" ;
2102 cols[3] = "BEAMNO" ;
2103 cols[4] = "IFNO" ;
2104 cols[5] = "INTERVAL" ;
2105 Table tab = table_->table().project( cols ) ;
2106
2107 if ( nrow == 0 )
2108 return ;
2109
2110 //nrow = idRec.nfields() ;
2111 nrow = idRec.size() ;
2112
2113 Double midTime ;
2114 Double interval ;
2115 String timeStr ;
2116
2117 // row base
2118 TableRow row( mssc ) ;
2119 TableRecord &trec = row.record() ;
2120 RecordFieldPtr<Int> antennaRF( trec, "ANTENNA_ID" ) ;
2121 RecordFieldPtr<Int> feedRF( trec, "FEED_ID" ) ;
2122 RecordFieldPtr<Int> spwRF( trec, "SPECTRAL_WINDOW_ID" ) ;
2123 RecordFieldPtr<Double> timeRF( trec, "TIME" ) ;
2124 RecordFieldPtr<Double> intervalRF( trec, "INTERVAL" ) ;
2125 RecordFieldPtr< Array<Float> > tsysRF( trec, "TSYS" ) ;
2126 RecordFieldPtr< Array<Float> > tcalRF( trec, "TCAL" ) ;
2127 RecordFieldPtr< Array<Float> > tsysspRF ;
2128 RecordFieldPtr< Array<Float> > tcalspRF ;
2129 if ( tsysSpec_ )
2130 tsysspRF.attachToRecord( trec, "TSYS_SPECTRUM" ) ;
2131 if ( tcalSpec_ )
2132 tcalspRF.attachToRecord( trec, "TCAL_SPECTRUM" ) ;
2133
2134 // ANTENNA_ID is always 0
2135 *antennaRF = 0 ;
2136
2137 Table sortedstt = stt.sort( "ID" ) ;
2138 ROArrayColumn<Float> tcalCol( sortedstt, "TCAL" ) ;
2139 ROTableColumn idCol( sortedstt, "ID" ) ;
2140 ROArrayColumn<Float> tsysCol( tab, "TSYS" ) ;
2141 ROTableColumn tcalidCol( tab, "TCAL_ID" ) ;
2142 ROTableColumn timeCol( tab, "TIME" ) ;
2143 ROTableColumn intervalCol( tab, "INTERVAL" ) ;
2144 ROTableColumn beamnoCol( tab, "BEAMNO" ) ;
2145 ROTableColumn ifnoCol( tab, "IFNO" ) ;
2146 map< Int,Vector<uInt> >::iterator itr0 = idRec.begin() ;
2147 map< Int,Vector<uInt> >::iterator itr1 = rowRec.begin() ;
2148 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
2149// double t1 = mathutil::gettimeofday_sec() ;
2150 Vector<uInt> ids = itr0->second ;
2151 itr0++ ;
2152// os_ << "ids = " << ids << LogIO::POST ;
2153 uInt npol = ids.size() ;
2154 Vector<uInt> rows = itr1->second ;
2155 itr1++ ;
2156// os_ << "rows = " << rows << LogIO::POST ;
2157 Vector<Double> atime( rows.nelements() ) ;
2158 Vector<Double> ainterval( rows.nelements() ) ;
2159 Vector<uInt> atcalid( rows.nelements() ) ;
2160 for( uInt jrow = 0 ; jrow < rows.nelements() ; jrow++ ) {
2161 atime[jrow] = (Double)timeCol.asdouble( rows[jrow] ) ;
2162 ainterval[jrow] = (Double)intervalCol.asdouble( rows[jrow] ) ;
2163 atcalid[jrow] = tcalidCol.asuInt( rows[jrow] ) ;
2164 }
2165 Vector<Float> dummy = tsysCol( rows[0] ) ;
2166 Matrix<Float> tsys( npol,dummy.nelements() ) ;
2167 tsys.row( 0 ) = dummy ;
2168 for ( uInt jrow = 1 ; jrow < npol ; jrow++ )
2169 tsys.row( jrow ) = tsysCol( rows[jrow] ) ;
2170
2171 // FEED_ID
2172 *feedRF = beamnoCol.asuInt( rows[0] ) ;
2173
2174 // SPECTRAL_WINDOW_ID
2175 *spwRF = ifnoCol.asuInt( rows[0] ) ;
2176
2177 // TIME and INTERVAL
2178 getValidTimeRange( midTime, interval, atime, ainterval ) ;
2179 *timeRF = midTime ;
2180 *intervalRF = interval ;
2181
2182 // TCAL and TSYS
2183 Matrix<Float> tcal ;
2184 Table t ;
2185 if ( idCol.asuInt( ids[0] ) == ids[0] ) {
2186// os_ << "sorted at irow=" << irow << " ids[0]=" << ids[0] << LogIO::POST ;
2187 Vector<Float> dummyC = tcalCol( ids[0] ) ;
2188 tcal.resize( npol, dummyC.size() ) ;
2189 tcal.row( 0 ) = dummyC ;
2190 }
2191 else {
2192// os_ << "NOT sorted at irow=" << irow << " ids[0]=" << ids[0] << LogIO::POST ;
2193 t = stt( stt.col("ID") == ids[0], 1 ) ;
2194 Vector<Float> dummyC = tcalCol( 0 ) ;
2195 tcal.resize( npol, dummyC.size(), True ) ;
2196 tcal.row( 0 ) = dummyC ;
2197 }
2198 if ( npol == 2 ) {
2199 if ( idCol.asuInt( ids[1] ) == ids[1] ) {
2200// os_ << "sorted at irow=" << irow << " ids[1]=" << ids[1] << LogIO::POST ;
2201 tcal.row( 1 ) = tcalCol( ids[1] ) ;
2202 }
2203 else {
2204// os_ << "NOT sorted at irow=" << irow << " ids[1]=" << ids[1] << LogIO::POST ;
2205 t = stt( stt.col("ID") == ids[1], 1 ) ;
2206 tcalCol.attach( t, "TCAL" ) ;
2207 tcal.row( 1 ) = tcalCol( 0 ) ;
2208 }
2209 }
2210 else if ( npol == 3 ) {
2211 if ( idCol.asuInt( ids[2] ) == ids[2] )
2212 tcal.row( 1 ) = tcalCol( ids[2] ) ;
2213 else {
2214 t = stt( stt.col("ID") == ids[2], 1 ) ;
2215 tcalCol.attach( t, "TCAL" ) ;
2216 tcal.row( 1 ) = tcalCol( 0 ) ;
2217 }
2218 if ( idCol.asuInt( ids[1] ) == ids[1] )
2219 tcal.row( 2 ) = tcalCol( ids[1] ) ;
2220 else {
2221 t = stt( stt.col("ID") == ids[1], 1 ) ;
2222 tcalCol.attach( t, "TCAL" ) ;
2223 tcal.row( 2 ) = tcalCol( 0 ) ;
2224 }
2225 }
2226 else if ( npol == 4 ) {
2227 if ( idCol.asuInt( ids[2] ) == ids[2] )
2228 tcal.row( 1 ) = tcalCol( ids[2] ) ;
2229 else {
2230 t = stt( stt.col("ID") == ids[2], 1 ) ;
2231 tcalCol.attach( t, "TCAL" ) ;
2232 tcal.row( 1 ) = tcalCol( 0 ) ;
2233 }
2234 if ( idCol.asuInt( ids[3] ) == ids[3] )
2235 tcal.row( 2 ) = tcalCol( ids[3] ) ;
2236 else {
2237 t = stt( stt.col("ID") == ids[3], 1 ) ;
2238 tcalCol.attach( t, "TCAL" ) ;
2239 tcal.row( 2 ) = tcalCol( 0 ) ;
2240 }
2241 if ( idCol.asuInt( ids[1] ) == ids[1] )
2242 tcal.row( 2 ) = tcalCol( ids[1] ) ;
2243 else {
2244 t = stt( stt.col("ID") == ids[1], 1 ) ;
2245 tcalCol.attach( t, "TCAL" ) ;
2246 tcal.row( 3 ) = tcalCol( 0 ) ;
2247 }
2248 }
2249 if ( tcalSpec_ ) {
2250 // put TCAL_SPECTRUM
2251 tcalspRF.define( tcal ) ;
2252 // set TCAL (mean of TCAL_SPECTRUM)
2253 Matrix<Float> tcalMean( npol, 1 ) ;
2254 for ( uInt iid = 0 ; iid < npol ; iid++ ) {
2255 tcalMean( iid, 0 ) = mean( tcal.row(iid) ) ;
2256 }
2257 // put TCAL
2258 tcalRF.define( tcalMean ) ;
2259 }
2260 else {
2261 // put TCAL
2262 tcalRF.define( tcal ) ;
2263 }
2264
2265 if ( tsysSpec_ ) {
2266 // put TSYS_SPECTRUM
2267 tsysspRF.define( tsys ) ;
2268 // set TSYS (mean of TSYS_SPECTRUM)
2269 Matrix<Float> tsysMean( npol, 1 ) ;
2270 for ( uInt iid = 0 ; iid < npol ; iid++ ) {
2271 tsysMean( iid, 0 ) = mean( tsys.row(iid) ) ;
2272 }
2273 // put TSYS
2274 tsysRF.define( tsysMean ) ;
2275 }
2276 else {
2277 // put TSYS
2278 tsysRF.define( tsys ) ;
2279 }
2280
2281 // add row
2282 mssc.addRow( 1, True ) ;
2283 row.put( mssc.nrow()-1 ) ;
2284
2285// double t2 = mathutil::gettimeofday_sec() ;
2286// os_ << irow << "th loop elapsed time = " << t2-t1 << "sec" << LogIO::POST ;
2287 }
2288
2289 //double endSec = mathutil::gettimeofday_sec() ;
2290 //os_ << "end MSWriter::fillSysCal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2291}
2292
2293void MSWriter::addFeed( Int id )
2294{
2295// double startSec = mathutil::gettimeofday_sec() ;
2296// os_ << "start MSWriter::addFeed() startSec=" << startSec << LogIO::POST ;
2297
2298 // add row
2299 MSFeed msFeed = mstable_->feed() ;
2300 msFeed.addRow( 1, True ) ;
2301 Int nrow = msFeed.nrow() ;
2302 Int numReceptors = 2 ;
2303 Vector<String> polType( numReceptors ) ;
2304 Matrix<Double> beamOffset( 2, numReceptors ) ;
2305 beamOffset = 0.0 ;
2306 Vector<Double> receptorAngle( numReceptors, 0.0 ) ;
2307 if ( polType_ == "linear" ) {
2308 polType[0] = "X" ;
2309 polType[1] = "Y" ;
2310 }
2311 else if ( polType_ == "circular" ) {
2312 polType[0] = "R" ;
2313 polType[1] = "L" ;
2314 }
2315 else {
2316 polType[0] = "X" ;
2317 polType[1] = "Y" ;
2318 }
2319 Matrix<Complex> polResponse( numReceptors, numReceptors, 0.0 ) ;
2320 for ( Int i = 0 ; i < numReceptors ; i++ )
2321 polResponse( i, i ) = 0.0 ;
2322
2323 MSFeedColumns msFeedCols( mstable_->feed() ) ;
2324
2325 msFeedCols.feedId().put( nrow-1, id ) ;
2326 msFeedCols.antennaId().put( nrow-1, 0 ) ;
2327 msFeedCols.numReceptors().put( nrow-1, numReceptors ) ;
2328 msFeedCols.polarizationType().put( nrow-1, polType ) ;
2329 msFeedCols.beamOffset().put( nrow-1, beamOffset ) ;
2330 msFeedCols.receptorAngle().put( nrow-1, receptorAngle ) ;
2331 msFeedCols.polResponse().put( nrow-1, polResponse ) ;
2332
2333// double endSec = mathutil::gettimeofday_sec() ;
2334// os_ << "end MSWriter::addFeed() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2335}
2336
2337void MSWriter::addSpectralWindow( Int spwid, Int freqid )
2338{
2339// double startSec = mathutil::gettimeofday_sec() ;
2340// os_ << "start MSWriter::addSpectralWindow() startSec=" << startSec << LogIO::POST ;
2341
2342 // add row
2343 MSSpectralWindow msSpw = mstable_->spectralWindow() ;
2344 while( (Int)msSpw.nrow() <= spwid ) {
2345 msSpw.addRow( 1, True ) ;
2346 }
2347
2348 MSSpWindowColumns msSpwCols( msSpw ) ;
2349
2350 STFrequencies stf = table_->frequencies() ;
2351
2352 // MEAS_FREQ_REF
2353 msSpwCols.measFreqRef().put( spwid, stf.getFrame( True ) ) ;
2354
2355 Double refpix ;
2356 Double refval ;
2357 Double inc ;
2358 stf.getEntry( refpix, refval, inc, (uInt)freqid ) ;
2359
2360 // NUM_CHAN
2361 //Int nchan = (Int)(refpix * 2) + 1 ;
2362 //if ( nchan == 0 )
2363 //nchan = 1 ;
2364 Int nchan = table_->nchan( spwid ) ;
2365 msSpwCols.numChan().put( spwid, nchan ) ;
2366
2367 // TOTAL_BANDWIDTH
2368 Double bw = nchan * abs( inc ) ;
2369 msSpwCols.totalBandwidth().put( spwid, bw ) ;
2370
2371 // REF_FREQUENCY
2372 Double refFreq = refval - refpix * inc ;
2373 msSpwCols.refFrequency().put( spwid, refFreq ) ;
2374
2375 // NET_SIDEBAND
2376 // tentative: USB->0, LSB->1
2377 Int netSideband = 0 ;
2378 if ( inc < 0 )
2379 netSideband = 1 ;
2380 msSpwCols.netSideband().put( spwid, netSideband ) ;
2381
2382 // RESOLUTION, CHAN_WIDTH, EFFECTIVE_BW
2383 Vector<Double> sharedDoubleArr( nchan, abs(inc) ) ;
2384 msSpwCols.resolution().put( spwid, sharedDoubleArr ) ;
2385 msSpwCols.chanWidth().put( spwid, sharedDoubleArr ) ;
2386 msSpwCols.effectiveBW().put( spwid, sharedDoubleArr ) ;
2387
2388 // CHAN_FREQ
2389 indgen( sharedDoubleArr, refFreq, inc ) ;
2390 msSpwCols.chanFreq().put( spwid, sharedDoubleArr ) ;
2391
2392// double endSec = mathutil::gettimeofday_sec() ;
2393// os_ << "end MSWriter::addSpectralWindow() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2394}
2395
2396void MSWriter::addField( Int fid, String fieldname, String srcname, Double t, Vector<Double> rate )
2397{
2398// double startSec = mathutil::gettimeofday_sec() ;
2399// os_ << "start MSWriter::addField() startSec=" << startSec << LogIO::POST ;
2400
2401 MSField msField = mstable_->field() ;
2402 while( (Int)msField.nrow() <= fid ) {
2403 msField.addRow( 1, True ) ;
2404 }
2405 MSFieldColumns msFieldCols( msField ) ;
2406
2407 // Access to SOURCE table
2408 MSSource msSrc = mstable_->source() ;
2409
2410 // fill target row
2411 msFieldCols.name().put( fid, fieldname ) ;
2412 msFieldCols.time().put( fid, t ) ;
2413 Int numPoly = 0 ;
2414 if ( anyNE( rate, 0.0 ) )
2415 numPoly = 1 ;
2416 msFieldCols.numPoly().put( fid, numPoly ) ;
2417 MSSourceIndex msSrcIdx( msSrc ) ;
2418 Int srcId = -1 ;
2419 Vector<Int> srcIdArr = msSrcIdx.matchSourceName( srcname ) ;
2420 if ( srcIdArr.size() != 0 ) {
2421 srcId = srcIdArr[0] ;
2422 MSSource msSrcSel = msSrc( msSrc.col("SOURCE_ID") == srcId, 1 ) ;
2423 ROMSSourceColumns msSrcCols( msSrcSel ) ;
2424 Vector<Double> srcDir = msSrcCols.direction()( 0 ) ;
2425 Matrix<Double> srcDirA( IPosition( 2, 2, 1+numPoly ) ) ;
2426// os_ << "srcDirA = " << srcDirA << LogIO::POST ;
2427// os_ << "sliced srcDirA = " << srcDirA.column( 0 ) << LogIO::POST ;
2428 srcDirA.column( 0 ) = srcDir ;
2429// os_ << "srcDirA = " << srcDirA << LogIO::POST ;
2430 if ( numPoly != 0 )
2431 srcDirA.column( 1 ) = rate ;
2432 msFieldCols.phaseDir().put( fid, srcDirA ) ;
2433 msFieldCols.referenceDir().put( fid, srcDirA ) ;
2434 msFieldCols.delayDir().put( fid, srcDirA ) ;
2435 }
2436 msFieldCols.sourceId().put( fid, srcId ) ;
2437
2438// double endSec = mathutil::gettimeofday_sec() ;
2439// os_ << "end MSWriter::addField() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2440}
2441
2442void MSWriter::addPointing( String &name, Double &me, Double &interval, Matrix<Double> &dir )
2443{
2444// double startSec = mathutil::gettimeofday_sec() ;
2445// os_ << "start MSWriter::addPointing() startSec=" << startSec << LogIO::POST ;
2446
2447 // access to POINTING subtable
2448 MSPointing msp = mstable_->pointing() ;
2449 uInt nrow = msp.nrow() ;
2450
2451 // add row
2452 msp.addRow( 1, True ) ;
2453
2454 // fill row
2455 TableRow row( msp ) ;
2456 TableRecord &rec = row.record() ;
2457 RecordFieldPtr<Int> antennaRF( rec, "ANTENNA_ID" ) ;
2458 *antennaRF = 0 ;
2459 RecordFieldPtr<Int> numpolyRF( rec, "NUM_POLY" ) ;
2460 *numpolyRF = dir.ncolumn() - 1 ;
2461 RecordFieldPtr<Double> timeRF( rec, "TIME" ) ;
2462 *timeRF = me ;
2463 RecordFieldPtr<Double> toriginRF( rec, "TIME_ORIGIN" ) ;
2464 *toriginRF = me ;
2465 RecordFieldPtr<Double> intervalRF( rec, "INTERVAL" ) ;
2466 *intervalRF = interval ;
2467 RecordFieldPtr<String> nameRF( rec, "NAME" ) ;
2468 *nameRF = name ;
2469 RecordFieldPtr<Bool> trackRF( rec, "TRACKING" ) ;
2470 *trackRF = True ;
2471 RecordFieldPtr< Array<Double> > dirRF( rec, "DIRECTION" ) ;
2472 *dirRF = dir ;
2473 RecordFieldPtr< Array<Double> > targetRF( rec, "TARGET" ) ;
2474 *targetRF = dir ;
2475 row.put( nrow ) ;
2476
2477// double endSec = mathutil::gettimeofday_sec() ;
2478// os_ << "end MSWriter::addPointing() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2479}
2480
2481Int MSWriter::addPolarization( Vector<Int> polnos )
2482{
2483// double startSec = mathutil::gettimeofday_sec() ;
2484// os_ << "start MSWriter::addPolarization() startSec=" << startSec << LogIO::POST ;
2485
2486// os_ << "polnos = " << polnos << LogIO::POST ;
2487 MSPolarization msPol = mstable_->polarization() ;
2488 uInt nrow = msPol.nrow() ;
2489
2490// // only 1 POLARIZATION row for 1 scantable
2491// if ( nrow > 0 )
2492// return 0 ;
2493
2494 Vector<Int> corrType = toCorrType( polnos ) ;
2495
2496 ROArrayColumn<Int> corrtCol( msPol, "CORR_TYPE" ) ;
2497 //Matrix<Int> corrTypeArr = corrtCol.getColumn() ;
2498 Int polid = -1 ;
2499 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
2500 Vector<Int> corrTypeArr = corrtCol( irow ) ;
2501 if ( corrType.nelements() == corrTypeArr.nelements()
2502 && allEQ( corrType, corrTypeArr ) ) {
2503 polid = irow ;
2504 break ;
2505 }
2506 }
2507
2508 if ( polid == -1 ) {
2509 MSPolarizationColumns msPolCols( msPol ) ;
2510
2511 // add row
2512 msPol.addRow( 1, True ) ;
2513 polid = (Int)nrow ;
2514
2515 // CORR_TYPE
2516 msPolCols.corrType().put( nrow, corrType ) ;
2517
2518 // NUM_CORR
2519 uInt npol = corrType.size() ;
2520 msPolCols.numCorr().put( nrow, npol ) ;
2521
2522 // CORR_PRODUCT
2523 Matrix<Int> corrProd( 2, npol, -1 ) ;
2524 if ( npol == 1 ) {
2525 corrProd = 0 ;
2526 }
2527 else if ( npol == 2 ) {
2528 corrProd.column( 0 ) = 0 ;
2529 corrProd.column( 1 ) = 1 ;
2530 }
2531 else {
2532 corrProd.column( 0 ) = 0 ;
2533 corrProd.column( 3 ) = 1 ;
2534 corrProd( 0,1 ) = 0 ;
2535 corrProd( 1,1 ) = 1 ;
2536 corrProd( 0,2 ) = 1 ;
2537 corrProd( 1,2 ) = 0 ;
2538 }
2539 msPolCols.corrProduct().put( nrow, corrProd ) ;
2540 }
2541
2542// double endSec = mathutil::gettimeofday_sec() ;
2543// os_ << "end MSWriter::addPolarization() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2544
2545 return polid ;
2546}
2547
2548Int MSWriter::addDataDescription( Int polid, Int spwid )
2549{
2550// double startSec = mathutil::gettimeofday_sec() ;
2551// os_ << "start MSWriter::addDataDescription() startSec=" << startSec << LogIO::POST ;
2552
2553 MSDataDescription msDataDesc = mstable_->dataDescription() ;
2554 uInt nrow = msDataDesc.nrow() ;
2555
2556 // only 1 POLARIZATION_ID for 1 scantable
2557 Int ddid = -1 ;
2558 ROScalarColumn<Int> spwCol( msDataDesc, "SPECTRAL_WINDOW_ID" ) ;
2559 Vector<Int> spwIds = spwCol.getColumn() ;
2560 //ROScalarColumn<Int> polCol( msDataDesc, "POLARIZATION_ID" ) ;
2561 //Vector<Int> polIds = polCol.getColumn() ;
2562 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
2563 //if ( spwid == spwIds[irow] && polid == polIds[irow] ) {
2564 if ( spwid == spwIds[irow] ) {
2565 ddid = irow ;
2566 break ;
2567 }
2568 }
2569// os_ << "ddid = " << ddid << LogIO::POST ;
2570
2571
2572 if ( ddid == -1 ) {
2573 msDataDesc.addRow( 1, True ) ;
2574 MSDataDescColumns msDataDescCols( msDataDesc ) ;
2575 msDataDescCols.polarizationId().put( nrow, polid ) ;
2576 msDataDescCols.spectralWindowId().put( nrow, spwid ) ;
2577 ddid = (Int)nrow ;
2578 }
2579
2580// double endSec = mathutil::gettimeofday_sec() ;
2581// os_ << "end MSWriter::addDataDescription() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2582
2583 return ddid ;
2584}
2585
2586Int MSWriter::addState( Int st, Int &subscan )
2587{
2588// double startSec = mathutil::gettimeofday_sec() ;
2589// os_ << "start MSWriter::addState() startSec=" << startSec << LogIO::POST ;
2590
2591 // access to STATE subtable
2592 MSState msState = mstable_->state() ;
2593 uInt nrow = msState.nrow() ;
2594
2595 String obsMode ;
2596 Bool isSignal ;
2597 Double tnoise ;
2598 Double tload ;
2599 queryType( st, obsMode, isSignal, tnoise, tload ) ;
2600// os_ << "obsMode = " << obsMode << " isSignal = " << isSignal << LogIO::POST ;
2601
2602 Int idx = -1 ;
2603 ROScalarColumn<String> obsModeCol( msState, "OBS_MODE" ) ;
2604 ROScalarColumn<Int> subscanCol( msState, "SUB_SCAN" ) ;
2605 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
2606 if ( obsModeCol(irow) == obsMode
2607 //&& sigCol(irow) == isSignal
2608 //&& refCol(irow) != isSignal
2609 && subscanCol(irow) == subscan ) {
2610 idx = irow ;
2611 break ;
2612 }
2613 }
2614 if ( idx == -1 ) {
2615 msState.addRow( 1, True ) ;
2616 TableRow row( msState ) ;
2617 TableRecord &rec = row.record() ;
2618 RecordFieldPtr<String> obsmodeRF( rec, "OBS_MODE" ) ;
2619 *obsmodeRF = obsMode ;
2620 RecordFieldPtr<Bool> sigRF( rec, "SIG" ) ;
2621 *sigRF = isSignal ;
2622 RecordFieldPtr<Bool> refRF( rec, "REF" ) ;
2623 *refRF = !isSignal ;
2624 RecordFieldPtr<Int> subscanRF( rec, "SUB_SCAN" ) ;
2625 *subscanRF = subscan ;
2626 RecordFieldPtr<Double> noiseRF( rec, "CAL" ) ;
2627 *noiseRF = tnoise ;
2628 RecordFieldPtr<Double> loadRF( rec, "LOAD" ) ;
2629 *loadRF = tload ;
2630 row.put( nrow ) ;
2631 idx = nrow ;
2632 }
2633 subscan++ ;
2634
2635// double endSec = mathutil::gettimeofday_sec() ;
2636// os_ << "end MSWriter::addState() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2637
2638 return idx ;
2639}
2640
2641Vector<Int> MSWriter::toCorrType( Vector<Int> polnos )
2642{
2643// double startSec = mathutil::gettimeofday_sec() ;
2644// os_ << "start MSWriter::toCorrType() startSec=" << startSec << LogIO::POST ;
2645
2646 uInt npol = polnos.size() ;
2647 Vector<Int> corrType( npol, Stokes::Undefined ) ;
2648
2649 if ( npol == 4 ) {
2650 if ( polType_ == "linear" ) {
2651 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
2652 if ( polnos[ipol] == 0 )
2653 corrType[ipol] = Stokes::XX ;
2654 else if ( polnos[ipol] == 1 )
2655 corrType[ipol] = Stokes::XY ;
2656 else if ( polnos[ipol] == 2 )
2657 corrType[ipol] = Stokes::YX ;
2658 else if ( polnos[ipol] == 3 )
2659 corrType[ipol] = Stokes::YY ;
2660 }
2661 }
2662 else if ( polType_ == "circular" ) {
2663 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
2664 if ( polnos[ipol] == 0 )
2665 corrType[ipol] = Stokes::RR ;
2666 else if ( polnos[ipol] == 1 )
2667 corrType[ipol] = Stokes::RL ;
2668 else if ( polnos[ipol] == 2 )
2669 corrType[ipol] = Stokes::LR ;
2670 else if ( polnos[ipol] == 3 )
2671 corrType[ipol] = Stokes::LL ;
2672 }
2673 }
2674 else if ( polType_ == "stokes" ) {
2675 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
2676 if ( polnos[ipol] == 0 )
2677 corrType[ipol] = Stokes::I ;
2678 else if ( polnos[ipol] == 1 )
2679 corrType[ipol] = Stokes::Q ;
2680 else if ( polnos[ipol] == 2 )
2681 corrType[ipol] = Stokes::U ;
2682 else if ( polnos[ipol] == 3 )
2683 corrType[ipol] = Stokes::V ;
2684 }
2685 }
2686 }
2687 else if ( npol == 2 ) {
2688 if ( polType_ == "linear" ) {
2689 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
2690 if ( polnos[ipol] == 0 )
2691 corrType[ipol] = Stokes::XX ;
2692 else if ( polnos[ipol] == 1 )
2693 corrType[ipol] = Stokes::YY ;
2694 }
2695 }
2696 else if ( polType_ == "circular" ) {
2697 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
2698 if ( polnos[ipol] == 0 )
2699 corrType[ipol] = Stokes::RR ;
2700 else if ( polnos[ipol] == 1 )
2701 corrType[ipol] = Stokes::LL ;
2702 }
2703 }
2704 else if ( polType_ == "stokes" ) {
2705 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
2706 if ( polnos[ipol] == 0 )
2707 corrType[ipol] = Stokes::I ;
2708 else if ( polnos[ipol] == 1 )
2709 corrType[ipol] = Stokes::V ;
2710 }
2711 }
2712 else if ( polType_ == "linpol" ) {
2713 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
2714 if ( polnos[ipol] == 1 )
2715 corrType[ipol] = Stokes::Plinear ;
2716 else if ( polnos[ipol] == 2 )
2717 corrType[ipol] = Stokes::Pangle ;
2718 }
2719 }
2720 }
2721 else if ( npol == 1 ) {
2722 if ( polType_ == "linear" )
2723 corrType[0] = Stokes::XX ;
2724 else if ( polType_ == "circular" )
2725 corrType[0] = Stokes::RR ;
2726 else if ( polType_ == "stokes" )
2727 corrType[0] = Stokes::I ;
2728 }
2729
2730// double endSec = mathutil::gettimeofday_sec() ;
2731// os_ << "end MSWriter::toCorrType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2732
2733 return corrType ;
2734}
2735
2736void MSWriter::getValidTimeRange( Double &me, Double &interval, Table &tab )
2737{
2738// double startSec = mathutil::gettimeofday_sec() ;
2739// os_ << "start MSWriter::getVaridTimeRange() startSec=" << startSec << LogIO::POST ;
2740
2741 // sort table
2742 //Table stab = tab.sort( "TIME" ) ;
2743
2744 ROScalarColumn<Double> timeCol( tab, "TIME" ) ;
2745 Vector<Double> timeArr = timeCol.getColumn() ;
2746 Double minTime ;
2747 Double maxTime ;
2748 minMax( minTime, maxTime, timeArr ) ;
2749 Double midTime = 0.5 * ( minTime + maxTime ) * 86400.0 ;
2750 // unit for TIME
2751 // Scantable: "d"
2752 // MS: "s"
2753 me = midTime ;
2754 interval = ( maxTime - minTime ) * 86400.0 ;
2755
2756// double endSec = mathutil::gettimeofday_sec() ;
2757// os_ << "end MSWriter::getValidTimeRange() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2758}
2759
2760void MSWriter::getValidTimeRange( Double &me, Double &interval, Vector<Double> &atime, Vector<Double> &ainterval )
2761{
2762// double startSec = mathutil::gettimeofday_sec() ;
2763// os_ << "start MSWriter::getVaridTimeRange() startSec=" << startSec << LogIO::POST ;
2764
2765 // sort table
2766 //Table stab = tab.sort( "TIME" ) ;
2767
2768 Double minTime ;
2769 Double maxTime ;
2770 minMax( minTime, maxTime, atime ) ;
2771 Double midTime = 0.5 * ( minTime + maxTime ) * 86400.0 ;
2772 // unit for TIME
2773 // Scantable: "d"
2774 // MS: "s"
2775 me = midTime ;
2776 interval = ( maxTime - minTime ) * 86400.0 + mean( ainterval ) ;
2777
2778// double endSec = mathutil::gettimeofday_sec() ;
2779// os_ << "end MSWriter::getValidTimeRange() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2780}
2781
2782//void MSWriter::queryType( Int type, String &stype, Bool &b )
2783void MSWriter::queryType( Int type, String &stype, Bool &b, Double &t, Double &l )
2784{
2785// double startSec = mathutil::gettimeofday_sec() ;
2786// os_ << "start MSWriter::queryType() startSec=" << startSec << LogIO::POST ;
2787
2788 // 2011/03/14 TN
2789 // OBS_MODE string of MS created by importasdm task is slightly
2790 // (but critically) changed.
2791 // 2011/05/20 TN
2792 // OBS_MODE string of MS created by importasdm task is changed
2793 // again (separator is now "#" instead of "_"
2794 String sep1="#" ;
2795 String sep2="," ;
2796 String target="OBSERVE_TARGET" ;
2797 String atmcal="CALIBRATE_TEMPERATURE" ;
2798 String onstr="ON_SOURCE" ;
2799 String offstr="OFF_SOURCE" ;
2800 String pswitch="POSITION_SWITCH" ;
2801 String nod="NOD" ;
2802 String fswitch="FREQUENCY_SWITCH" ;
2803 String sigstr="SIG" ;
2804 String refstr="REF" ;
2805 String unspecified="UNSPECIFIED" ;
2806 String ftlow="LOWER" ;
2807 String fthigh="HIGHER" ;
2808 switch ( type ) {
2809 case SrcType::PSON:
2810 //stype = "OBSERVE_TARGET_ON_SOURCE,POSITION_SWITCH" ;
2811 stype = target+sep1+onstr+sep2+pswitch ;
2812 b = True ;
2813 t = 0.0 ;
2814 l = 0.0 ;
2815 break ;
2816 case SrcType::PSOFF:
2817 //stype = "OBSERVE_TARGET_OFF_SOURCE,POSITION_SWITCH" ;
2818 stype = target+sep1+offstr+sep2+pswitch ;
2819 b = False ;
2820 t = 0.0 ;
2821 l = 0.0 ;
2822 break ;
2823 case SrcType::NOD:
2824 //stype = "OBSERVE_TARGET_ON_SOURCE,NOD" ;
2825 stype = target+sep1+onstr+sep2+nod ;
2826 b = True ;
2827 t = 0.0 ;
2828 l = 0.0 ;
2829 break ;
2830 case SrcType::FSON:
2831 //stype = "OBSERVE_TARGET_ON_SOURCE,FREQUENCY_SWITCH_SIG" ;
2832 stype = target+sep1+onstr+sep2+fswitch+sep1+sigstr ;
2833 b = True ;
2834 t = 0.0 ;
2835 l = 0.0 ;
2836 break ;
2837 case SrcType::FSOFF:
2838 //stype = "OBSERVE_TARGET_ON_SOURCE,FREQUENCY_SWITCH_REF" ;
2839 stype = target+sep1+onstr+sep2+fswitch+sep1+refstr ;
2840 b = False ;
2841 t = 0.0 ;
2842 l = 0.0 ;
2843 break ;
2844 case SrcType::SKY:
2845 //stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,UNSPECIFIED" ;
2846 stype = atmcal+sep1+offstr+sep2+unspecified ;
2847 b = False ;
2848 t = 0.0 ;
2849 l = 1.0 ;
2850 break ;
2851 case SrcType::HOT:
2852 //stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,UNSPECIFIED" ;
2853 stype = atmcal+sep1+offstr+sep2+unspecified ;
2854 b = False ;
2855 t = 0.0 ;
2856 l = 1.0 ;
2857 break ;
2858 case SrcType::WARM:
2859 //stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,UNSPECIFIED" ;
2860 stype = atmcal+sep1+offstr+sep2+unspecified ;
2861 t = 0.0 ;
2862 b = False ;
2863 l = 1.0 ;
2864 break ;
2865 case SrcType::COLD:
2866 //stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,UNSPECIFIED" ;
2867 stype = atmcal+sep1+offstr+sep2+unspecified ;
2868 b = False ;
2869 t = 0.0 ;
2870 l = 1.0 ;
2871 break ;
2872 case SrcType::PONCAL:
2873 //stype = "CALIBRATE_TEMPERATURE_ON_SOURCE,POSITION_SWITCH" ;
2874 stype = atmcal+sep1+onstr+sep2+pswitch ;
2875 b = True ;
2876 t = 1.0 ;
2877 l = 0.0 ;
2878 break ;
2879 case SrcType::POFFCAL:
2880 //stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,POSITION_SWITCH" ;
2881 stype = atmcal+sep1+offstr+sep2+pswitch ;
2882 b = False ;
2883 t = 1.0 ;
2884 l = 0.0 ;
2885 break ;
2886 case SrcType::NODCAL:
2887 //stype = "CALIBRATE_TEMPERATURE_ON_SOURCE,NOD" ;
2888 stype = atmcal+sep1+onstr+sep2+nod ;
2889 b = True ;
2890 t = 1.0 ;
2891 l = 0.0 ;
2892 break ;
2893 case SrcType::FONCAL:
2894 //stype = "CALIBRATE_TEMPERATURE_ON_SOURCE,FREQUENCY_SWITCH_SIG" ;
2895 stype = atmcal+sep1+onstr+sep2+fswitch+sep1+sigstr ;
2896 b = True ;
2897 t = 1.0 ;
2898 l = 0.0 ;
2899 break ;
2900 case SrcType::FOFFCAL:
2901 //stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,FREQUENCY_SWITCH_REF" ;
2902 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+refstr ;
2903 b = False ;
2904 t = 1.0 ;
2905 l = 0.0 ;
2906 break ;
2907 case SrcType::FSLO:
2908 //stype = "OBSERVE_TARGET_ON_SOURCE,FREQUENCY_SWITCH_LOWER" ;
2909 stype = target+sep1+onstr+sep2+fswitch+sep1+ftlow ;
2910 b = True ;
2911 t = 0.0 ;
2912 l = 0.0 ;
2913 break ;
2914 case SrcType::FLOOFF:
2915 //stype = "OBSERVE_TARGET_OFF_SOURCE,FREQUENCY_SWITCH_LOWER" ;
2916 stype = target+sep1+offstr+sep2+fswitch+sep1+ftlow ;
2917 b = False ;
2918 t = 0.0 ;
2919 l = 0.0 ;
2920 break ;
2921 case SrcType::FLOSKY:
2922 //stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,FREQUENCY_SWITCH_LOWER" ;
2923 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
2924 b = False ;
2925 t = 0.0 ;
2926 l = 1.0 ;
2927 break ;
2928 case SrcType::FLOHOT:
2929 //stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,FREQUENCY_SWITCH_LOWER" ;
2930 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
2931 b = False ;
2932 t = 0.0 ;
2933 l = 1.0 ;
2934 break ;
2935 case SrcType::FLOWARM:
2936 //stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,FREQUENCY_SWITCH_LOWER" ;
2937 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
2938 b = False ;
2939 t = 0.0 ;
2940 l = 1.0 ;
2941 break ;
2942 case SrcType::FLOCOLD:
2943 //stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,FREQUENCY_SWITCH_LOWER" ;
2944 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
2945 b = False ;
2946 t = 0.0 ;
2947 l = 1.0 ;
2948 break ;
2949 case SrcType::FSHI:
2950 //stype = "OBSERVE_TARGET_ON_SOURCE,FREQUENCY_SWITCH_HIGHER" ;
2951 stype = target+sep1+onstr+sep2+fswitch+sep1+fthigh ;
2952 b = True ;
2953 t = 0.0 ;
2954 l = 0.0 ;
2955 break ;
2956 case SrcType::FHIOFF:
2957 //stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,FREQUENCY_SWITCH_HIGHER" ;
2958 stype = target+sep1+offstr+sep2+fswitch+sep1+fthigh ;
2959 b = False ;
2960 t = 0.0 ;
2961 l = 0.0 ;
2962 break ;
2963 case SrcType::FHISKY:
2964 //stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,FREQUENCY_SWITCH_HIGHER" ;
2965 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
2966 b = False ;
2967 t = 0.0 ;
2968 l = 1.0 ;
2969 break ;
2970 case SrcType::FHIHOT:
2971 //stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,FREQUENCY_SWITCH_HIGHER" ;
2972 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
2973 b = False ;
2974 t = 0.0 ;
2975 l = 1.0 ;
2976 break ;
2977 case SrcType::FHIWARM:
2978 //stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,FREQUENCY_SWITCH_HIGHER" ;
2979 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
2980 b = False ;
2981 t = 0.0 ;
2982 l = 1.0 ;
2983 break ;
2984 case SrcType::FHICOLD:
2985 //stype = "CALIBRATE_TEMPERATURE_OFF_SOURCE,FREQUENCY_SWITCH_HIGHER" ;
2986 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
2987 b = False ;
2988 t = 0.0 ;
2989 l = 1.0 ;
2990 break ;
2991 case SrcType::SIG:
2992 //stype = "OBSERVE_TARGET_ON_SOURCE,UNSPECIFIED" ;
2993 stype = target+sep1+onstr+sep2+unspecified ;
2994 b = True ;
2995 t = 0.0 ;
2996 l = 0.0 ;
2997 break ;
2998 case SrcType::REF:
2999 //stype = "OBSERVE_TARGET_ON_SOURCE,UNSPECIFIED" ;
3000 stype = target+sep1+offstr+sep2+unspecified ;
3001 b = False ;
3002 t = 0.0 ;
3003 l = 0.0 ;
3004 break ;
3005 default:
3006 //stype = "UNSPECIFIED" ;
3007 stype = unspecified ;
3008 b = True ;
3009 t = 0.0 ;
3010 l = 0.0 ;
3011 break ;
3012 }
3013
3014// double endSec = mathutil::gettimeofday_sec() ;
3015// os_ << "end MSWriter::queryType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
3016}
3017
3018Double MSWriter::getDishDiameter( String antname )
3019{
3020 Double diameter = 0.0 ;
3021
3022 antname.upcase() ;
3023
3024 if ( antname.matches( Regex( "DV[0-9]+$" ) )
3025 || antname.matches( Regex( "DA[0-9]+$" ) )
3026 || antname.matches( Regex( "PM[0-9]+$" ) ) )
3027 diameter = 12.0 ;
3028 else if ( antname.matches( Regex( "CM[0-9]+$" ) ) )
3029 diameter = 7.0 ;
3030 else if ( antname.contains( "GBT" ) )
3031 diameter = 104.9 ;
3032 else if ( antname.contains( "MOPRA" ) )
3033 diameter = 22.0 ;
3034 else if ( antname.contains( "PKS" ) || antname.contains( "PARKS" ) )
3035 diameter = 64.0 ;
3036 else if ( antname.contains( "TIDBINBILLA" ) )
3037 diameter = 70.0 ;
3038 else if ( antname.contains( "CEDUNA" ) )
3039 diameter = 30.0 ;
3040 else if ( antname.contains( "HOBART" ) )
3041 diameter = 26.0 ;
3042 else if ( antname.contains( "APEX" ) )
3043 diameter = 12.0 ;
3044 else if ( antname.contains( "ASTE" ) )
3045 diameter = 10.0 ;
3046 else if ( antname.contains( "NRO" ) )
3047 diameter = 45.0 ;
3048 else
3049 diameter = 1.0 ;
3050
3051 return diameter ;
3052}
3053
3054void MSWriter::infillSpectralWindow()
3055{
3056 MSSpectralWindow msSpw = mstable_->spectralWindow() ;
3057 MSSpWindowColumns msSpwCols( msSpw ) ;
3058 uInt nrow = msSpw.nrow() ;
3059
3060 ScalarColumn<Int> measFreqRefCol = msSpwCols.measFreqRef() ;
3061 ArrayColumn<Double> chanFreqCol = msSpwCols.chanFreq() ;
3062 ArrayColumn<Double> chanWidthCol = msSpwCols.chanWidth() ;
3063 ArrayColumn<Double> effectiveBWCol = msSpwCols.effectiveBW() ;
3064 ArrayColumn<Double> resolutionCol = msSpwCols.resolution() ;
3065 Vector<Double> dummy( 1, 0.0 ) ;
3066 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
3067 if ( !(chanFreqCol.isDefined( irow )) ) {
3068 measFreqRefCol.put( irow, 1 ) ;
3069 chanFreqCol.put( irow, dummy ) ;
3070 chanWidthCol.put( irow, dummy ) ;
3071 effectiveBWCol.put( irow, dummy ) ;
3072 resolutionCol.put( irow, dummy ) ;
3073 }
3074 }
3075
3076}
3077
3078}
Note: See TracBrowser for help on using the repository browser.