source: trunk/src/MSWriter.cpp@ 2296

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

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Removed unused methods in MSWriter.
Cleaned up comments.


File size: 69.0 KB
RevLine 
[1974]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//
[2291]13#include <assert.h>
[1974]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>
[2016]20#include <casa/Arrays/Cube.h>
[1974]21
[1975]22#include <tables/Tables/ExprNode.h>
[1974]23#include <tables/Tables/TableDesc.h>
24#include <tables/Tables/SetupNewTab.h>
25#include <tables/Tables/TableIter.h>
26#include <tables/Tables/RefRows.h>
[2016]27#include <tables/Tables/TableRow.h>
[1974]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>
[1975]33#include <ms/MeasurementSets/MSSourceIndex.h>
[1974]34
35#include "MSWriter.h"
36#include "STHeader.h"
37#include "STFrequencies.h"
[1975]38#include "STMolecules.h"
[1977]39#include "STTcal.h"
[2258]40#include "MathUtils.h"
[2291]41#include "TableTraverse.h"
[1974]42
43using namespace casa ;
[2016]44using namespace std ;
[1974]45
[2016]46namespace asap {
[1974]47
[2291]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
[1974]1333MSWriter::MSWriter(CountedPtr<Scantable> stable)
[2016]1334 : table_(stable),
[2019]1335 isTcal_(False),
1336 isWeather_(False),
1337 tcalSpec_(False),
1338 tsysSpec_(False),
[2016]1339 ptTabName_("")
[1974]1340{
1341 os_ = LogIO() ;
1342 os_.origin( LogOrigin( "MSWriter", "MSWriter()", WHERE ) ) ;
[2026]1343// os_ << "MSWriter::MSWriter()" << LogIO::POST ;
[1974]1344
1345 // initialize writer
1346 init() ;
1347}
1348
1349MSWriter::~MSWriter()
1350{
1351 os_.origin( LogOrigin( "MSWriter", "~MSWriter()", WHERE ) ) ;
[2026]1352// os_ << "MSWriter::~MSWriter()" << LogIO::POST ;
[2016]1353
1354 if ( mstable_ != 0 )
1355 delete mstable_ ;
[1974]1356}
[2291]1357
[1974]1358bool MSWriter::write(const string& filename, const Record& rec)
1359{
1360 os_.origin( LogOrigin( "MSWriter", "write()", WHERE ) ) ;
[2291]1361 //double startSec = mathutil::gettimeofday_sec() ;
1362 //os_ << "start MSWriter::write() startSec=" << startSec << LogIO::POST ;
[1974]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
[1975]1402 // PROCESSOR
1403 fillProcessor() ;
1404
1405 // SOURCE
1406 fillSource() ;
1407
[1977]1408 // WEATHER
[2019]1409 if ( isWeather_ )
1410 fillWeather() ;
[1977]1411
[2291]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() ;
[2016]1444
[2291]1445 // SYSCAL
1446 if ( isTcal_ )
1447 fillSysCal( idRec, rowRec ) ;
[1974]1448 }
[2291]1449 /***
1450 * End iteration using TableVisitor
1451 ***/
[1974]1452
[2022]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
[2016]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 }
[1974]1478
[2291]1479 //double endSec = mathutil::gettimeofday_sec() ;
1480 //os_ << "end MSWriter::write() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
[1975]1481
[1974]1482 return True ;
1483}
[2291]1484
[1974]1485void MSWriter::init()
1486{
1487// os_.origin( LogOrigin( "MSWriter", "init()", WHERE ) ) ;
[2258]1488// double startSec = mathutil::gettimeofday_sec() ;
[2026]1489// os_ << "start MSWriter::init() startSec=" << startSec << LogIO::POST ;
[1974]1490
1491 // access to scantable
1492 header_ = table_->getHeader() ;
1493
1494 // FLOAT_DATA? or DATA?
1495 if ( header_.npol > 2 ) {
[1977]1496 useFloatData_ = False ;
1497 useData_ = True ;
[1974]1498 }
1499 else {
[1977]1500 useFloatData_ = True ;
1501 useData_ = False ;
[1974]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
[2019]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 ) ) {
[2026]1523 os_ << "TCAL table exists: nrow=" << table_->tcal().table().nrow() << LogIO::POST ;
[2019]1524 isTcal_ = True ;
1525 }
1526 else {
[2026]1527 os_ << "No TCAL rows" << LogIO::POST ;
[2019]1528 }
1529 }
1530 else {
[2026]1531 os_ << "No TCAL rows" << LogIO::POST ;
[2019]1532 }
1533 if ( table_->weather().table().nrow() != 0 ) {
1534 ROTableColumn col( table_->weather().table(), "TEMPERATURE" ) ;
1535 if ( col.isDefined( 0 ) ) {
[2026]1536 os_ << "WEATHER table exists: nrow=" << table_->weather().table().nrow() << LogIO::POST ;
[2019]1537 isWeather_ =True ;
1538 }
1539 else {
[2026]1540 os_ << "No WEATHER rows" << LogIO::POST ;
[2019]1541 }
1542 }
1543 else {
[2026]1544 os_ << "No WEATHER rows" << LogIO::POST ;
[2019]1545 }
1546
[1977]1547 // Are TCAL_SPECTRUM and TSYS_SPECTRUM necessary?
[2019]1548 if ( isTcal_ && header_.nchan != 1 ) {
[1977]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 }
[2016]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
[2258]1580// double endSec = mathutil::gettimeofday_sec() ;
[2026]1581// os_ << "end MSWriter::init() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
[1974]1582}
1583
1584void MSWriter::setupMS()
1585{
1586// os_.origin( LogOrigin( "MSWriter", "setupMS()", WHERE ) ) ;
[2258]1587// double startSec = mathutil::gettimeofday_sec() ;
[2026]1588// os_ << "start MSWriter::setupMS() startSec=" << startSec << LogIO::POST ;
[2185]1589
1590 String dunit = table_->getHeader().fluxunit ;
[1974]1591
1592 TableDesc msDesc = MeasurementSet::requiredTableDesc() ;
[1977]1593 if ( useFloatData_ )
1594 MeasurementSet::addColumnToDesc( msDesc, MSMainEnums::FLOAT_DATA, 2 ) ;
1595 else if ( useData_ )
1596 MeasurementSet::addColumnToDesc( msDesc, MSMainEnums::DATA, 2 ) ;
[1974]1597
1598 SetupNewTable newtab( filename_, msDesc, Table::New ) ;
1599
1600 mstable_ = new MeasurementSet( newtab ) ;
1601
[2185]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
[1974]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() ;
[1975]1659 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::TRANSITION, 1 ) ;
1660 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::REST_FREQUENCY, 1 ) ;
1661 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::SYSVEL, 1 ) ;
[1974]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() ;
[1977]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 ) ;
[1974]1681 SetupNewTable sysCalTab( mstable_->sysCalTableName(), sysCalDesc, Table::New ) ;
1682 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SYSCAL ), Table( sysCalTab ) ) ;
1683
1684 TableDesc weatherDesc = MSWeather::requiredTableDesc() ;
[1977]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 ) ;
[1974]1690 SetupNewTable weatherTab( mstable_->weatherTableName(), weatherDesc, Table::New ) ;
1691 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::WEATHER ), Table( weatherTab ) ) ;
1692
1693 mstable_->initRefs() ;
1694
[2258]1695// double endSec = mathutil::gettimeofday_sec() ;
[2026]1696// os_ << "end MSWriter::setupMS() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
[1974]1697}
1698
1699void MSWriter::fillObservation()
1700{
[2291]1701 //double startSec = mathutil::gettimeofday_sec() ;
1702 //os_ << "start MSWriter::fillObservation() startSec=" << startSec << LogIO::POST ;
[2016]1703
[1974]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 }
[2026]1719// os_ << "telescopeName = " << telescopeName << LogIO::POST ;
[1974]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 ) ;
[2016]1729
[2291]1730 //double endSec = mathutil::gettimeofday_sec() ;
1731 //os_ << "end MSWriter::fillObservation() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
[1974]1732}
1733
[2291]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
[1974]1768void MSWriter::fillAntenna()
1769{
[2291]1770 //double startSec = mathutil::gettimeofday_sec() ;
1771 //os_ << "start MSWriter::fillAntenna() startSec=" << startSec << LogIO::POST ;
[2016]1772
[1974]1773 // only 1 row
[2291]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( "//" ) ;
[1974]1781 String antennaName ;
1782 String stationName ;
1783 if ( pos != String::npos ) {
[2291]1784 stationName = hAntName.substr( 0, pos ) ;
1785 hAntName = hAntName.substr( pos+2 ) ;
[1974]1786 }
[2291]1787 pos = hAntName.find( "@" ) ;
[1974]1788 if ( pos != String::npos ) {
[2291]1789 antennaName = hAntName.substr( 0, pos ) ;
1790 stationName = hAntName.substr( pos+1 ) ;
[1974]1791 }
1792 else {
[2291]1793 antennaName = hAntName ;
[1974]1794 }
[2291]1795 Vector<Double> antpos = keys.asArrayDouble( "AntennaPosition" ) ;
[1974]1796
[2291]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 ) ;
[1974]1816
[2291]1817 //double endSec = mathutil::gettimeofday_sec() ;
1818 //os_ << "end MSWriter::fillAntenna() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1819}
[1974]1820
[1975]1821void MSWriter::fillProcessor()
1822{
[2258]1823// double startSec = mathutil::gettimeofday_sec() ;
[2026]1824// os_ << "start MSWriter::fillProcessor() startSec=" << startSec << LogIO::POST ;
[1975]1825
1826 // only add empty 1 row
1827 MSProcessor msProc = mstable_->processor() ;
1828 msProc.addRow( 1, True ) ;
[2016]1829
[2258]1830// double endSec = mathutil::gettimeofday_sec() ;
[2026]1831// os_ << "end MSWriter::fillProcessor() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
[1975]1832}
1833
1834void MSWriter::fillSource()
1835{
[2258]1836// double startSec = mathutil::gettimeofday_sec() ;
[2026]1837// os_ << "start MSWriter::fillSource() startSec=" << startSec << LogIO::POST ;
[2016]1838
[1975]1839 // access to MS SOURCE subtable
1840 MSSource msSrc = mstable_->source() ;
1841
1842 // access to MOLECULE subtable
1843 STMolecules stm = table_->molecules() ;
1844
1845 Int srcId = 0 ;
1846 Vector<Double> restFreq ;
1847 Vector<String> molName ;
1848 Vector<String> fMolName ;
[2016]1849
1850 // row based
1851 TableRow row( msSrc ) ;
1852 TableRecord &rec = row.record() ;
1853 RecordFieldPtr<Int> srcidRF( rec, "SOURCE_ID" ) ;
1854 RecordFieldPtr<String> nameRF( rec, "NAME" ) ;
1855 RecordFieldPtr< Array<Double> > srcpmRF( rec, "PROPER_MOTION" ) ;
1856 RecordFieldPtr< Array<Double> > srcdirRF( rec, "DIRECTION" ) ;
1857 RecordFieldPtr<Int> numlineRF( rec, "NUM_LINES" ) ;
1858 RecordFieldPtr< Array<Double> > restfreqRF( rec, "REST_FREQUENCY" ) ;
1859 RecordFieldPtr< Array<Double> > sysvelRF( rec, "SYSVEL" ) ;
1860 RecordFieldPtr< Array<String> > transitionRF( rec, "TRANSITION" ) ;
1861 RecordFieldPtr<Double> timeRF( rec, "TIME" ) ;
1862 RecordFieldPtr<Double> intervalRF( rec, "INTERVAL" ) ;
1863 RecordFieldPtr<Int> spwidRF( rec, "SPECTRAL_WINDOW_ID" ) ;
1864
[1975]1865 //
1866 // ITERATION: SRCNAME
1867 //
1868 TableIterator iter0( table_->table(), "SRCNAME" ) ;
1869 while( !iter0.pastEnd() ) {
[2016]1870 //Table t0( iter0.table() ) ;
1871 Table t0 = iter0.table() ;
[1975]1872
1873 // get necessary information
1874 ROScalarColumn<String> srcNameCol( t0, "SRCNAME" ) ;
1875 String srcName = srcNameCol( 0 ) ;
1876 ROArrayColumn<Double> sharedDArrRCol( t0, "SRCPROPERMOTION" ) ;
1877 Vector<Double> srcPM = sharedDArrRCol( 0 ) ;
1878 sharedDArrRCol.attach( t0, "SRCDIRECTION" ) ;
1879 Vector<Double> srcDir = sharedDArrRCol( 0 ) ;
1880 ROScalarColumn<Double> srcVelCol( t0, "SRCVELOCITY" ) ;
1881 Double srcVel = srcVelCol( 0 ) ;
[2291]1882 srcRec_.define( srcName, srcId ) ;
[1975]1883
[2016]1884 // NAME
1885 *nameRF = srcName ;
1886
1887 // SOURCE_ID
1888 *srcidRF = srcId ;
1889
1890 // PROPER_MOTION
1891 *srcpmRF = srcPM ;
1892
1893 // DIRECTION
1894 *srcdirRF = srcDir ;
1895
[1975]1896 //
1897 // ITERATION: MOLECULE_ID
1898 //
1899 TableIterator iter1( t0, "MOLECULE_ID" ) ;
1900 while( !iter1.pastEnd() ) {
[2016]1901 //Table t1( iter1.table() ) ;
1902 Table t1 = iter1.table() ;
[1975]1903
1904 // get necessary information
1905 ROScalarColumn<uInt> molIdCol( t1, "MOLECULE_ID" ) ;
1906 uInt molId = molIdCol( 0 ) ;
1907 stm.getEntry( restFreq, molName, fMolName, molId ) ;
1908
[2016]1909 uInt numFreq = restFreq.size() ;
1910
1911 // NUM_LINES
1912 *numlineRF = numFreq ;
1913
1914 // REST_FREQUENCY
1915 *restfreqRF = restFreq ;
1916
1917 // TRANSITION
[2019]1918 //*transitionRF = fMolName ;
1919 Vector<String> transition ;
1920 if ( fMolName.size() != 0 ) {
1921 transition = fMolName ;
1922 }
1923 else if ( molName.size() != 0 ) {
1924 transition = molName ;
1925 }
1926 else {
1927 transition.resize( numFreq ) ;
1928 transition = "" ;
1929 }
1930 *transitionRF = transition ;
[2016]1931
1932 // SYSVEL
1933 Vector<Double> sysvelArr( numFreq, srcVel ) ;
1934 *sysvelRF = sysvelArr ;
1935
[1975]1936 //
1937 // ITERATION: IFNO
1938 //
1939 TableIterator iter2( t1, "IFNO" ) ;
1940 while( !iter2.pastEnd() ) {
[2016]1941 //Table t2( iter2.table() ) ;
1942 Table t2 = iter2.table() ;
1943 uInt nrow = msSrc.nrow() ;
[1975]1944
1945 // get necessary information
1946 ROScalarColumn<uInt> ifNoCol( t2, "IFNO" ) ;
1947 uInt ifno = ifNoCol( 0 ) ; // IFNO = SPECTRAL_WINDOW_ID
[2016]1948 Double midTime ;
[1975]1949 Double interval ;
[2016]1950 getValidTimeRange( midTime, interval, t2 ) ;
[1975]1951
1952 // fill SPECTRAL_WINDOW_ID
[2016]1953 *spwidRF = ifno ;
[1975]1954
1955 // fill TIME, INTERVAL
[2016]1956 *timeRF = midTime ;
1957 *intervalRF = interval ;
[1975]1958
[2016]1959 // add row
1960 msSrc.addRow( 1, True ) ;
1961 row.put( nrow ) ;
1962
[1975]1963 iter2.next() ;
1964 }
1965
1966 iter1.next() ;
1967 }
1968
1969 // increment srcId if SRCNAME changed
1970 srcId++ ;
1971
1972 iter0.next() ;
1973 }
[2016]1974
[2258]1975// double endSec = mathutil::gettimeofday_sec() ;
[2026]1976// os_ << "end MSWriter::fillSource() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
[1975]1977}
1978
[1977]1979void MSWriter::fillWeather()
1980{
[2258]1981// double startSec = mathutil::gettimeofday_sec() ;
[2026]1982// os_ << "start MSWriter::fillWeather() startSec=" << startSec << LogIO::POST ;
[1977]1983
1984 // access to MS WEATHER subtable
1985 MSWeather msw = mstable_->weather() ;
1986
1987 // access to WEATHER subtable
1988 Table stw = table_->weather().table() ;
1989 uInt nrow = stw.nrow() ;
1990
1991 if ( nrow == 0 )
1992 return ;
1993
1994 msw.addRow( nrow, True ) ;
1995 MSWeatherColumns mswCols( msw ) ;
1996
1997 // ANTENNA_ID is always 0
1998 Vector<Int> antIdArr( nrow, 0 ) ;
1999 mswCols.antennaId().putColumn( antIdArr ) ;
2000
2001 // fill weather status
2002 ROScalarColumn<Float> sharedFloatCol( stw, "TEMPERATURE" ) ;
2003 mswCols.temperature().putColumn( sharedFloatCol ) ;
2004 sharedFloatCol.attach( stw, "PRESSURE" ) ;
2005 mswCols.pressure().putColumn( sharedFloatCol ) ;
2006 sharedFloatCol.attach( stw, "HUMIDITY" ) ;
2007 mswCols.relHumidity().putColumn( sharedFloatCol ) ;
2008 sharedFloatCol.attach( stw, "WINDSPEED" ) ;
2009 mswCols.windSpeed().putColumn( sharedFloatCol ) ;
2010 sharedFloatCol.attach( stw, "WINDAZ" ) ;
2011 mswCols.windDirection().putColumn( sharedFloatCol ) ;
2012
2013 // fill TIME and INTERVAL
[2016]2014 Double midTime ;
[1977]2015 Double interval ;
2016 Vector<Double> intervalArr( nrow, 0.0 ) ;
2017 TableIterator iter( table_->table(), "WEATHER_ID" ) ;
2018 while( !iter.pastEnd() ) {
[2016]2019 //Table tab( iter.table() ) ;
2020 Table tab = iter.table() ;
[1977]2021
2022 ROScalarColumn<uInt> widCol( tab, "WEATHER_ID" ) ;
2023 uInt wid = widCol( 0 ) ;
2024
[2016]2025 getValidTimeRange( midTime, interval, tab ) ;
2026 mswCols.time().put( wid, midTime ) ;
[1977]2027 intervalArr[wid] = interval ;
2028
2029 iter.next() ;
2030 }
2031 mswCols.interval().putColumn( intervalArr ) ;
[2016]2032
[2258]2033// double endSec = mathutil::gettimeofday_sec() ;
[2026]2034// os_ << "end MSWriter::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
[1977]2035}
2036
[2291]2037void MSWriter::fillSysCal( map< Int,Vector<uInt> > &idRec,
2038 map< Int,Vector<uInt> > &rowRec )
[1977]2039{
[2291]2040 //double startSec = mathutil::gettimeofday_sec() ;
2041 //os_ << "start MSWriter::fillSysCal() startSec=" << startSec << LogIO::POST ;
[1977]2042
[2291]2043 //idRec.print( cout ) ;
[1977]2044
2045 // access to MS SYSCAL subtable
2046 MSSysCal mssc = mstable_->sysCal() ;
2047
2048 // access to TCAL subtable
[2016]2049 Table stt = table_->tcal().table() ;
2050 uInt nrow = stt.nrow() ;
[1977]2051
[2016]2052 // access to MAIN table
[2018]2053 Block<String> cols( 6 ) ;
[2016]2054 cols[0] = "TIME" ;
2055 cols[1] = "TCAL_ID" ;
2056 cols[2] = "TSYS" ;
2057 cols[3] = "BEAMNO" ;
2058 cols[4] = "IFNO" ;
[2018]2059 cols[5] = "INTERVAL" ;
[2016]2060 Table tab = table_->table().project( cols ) ;
2061
[1977]2062 if ( nrow == 0 )
2063 return ;
2064
[2291]2065 //nrow = idRec.nfields() ;
2066 nrow = idRec.size() ;
[1977]2067
[2016]2068 Double midTime ;
[1977]2069 Double interval ;
2070 String timeStr ;
2071
[2016]2072 // row base
2073 TableRow row( mssc ) ;
2074 TableRecord &trec = row.record() ;
2075 RecordFieldPtr<Int> antennaRF( trec, "ANTENNA_ID" ) ;
2076 RecordFieldPtr<Int> feedRF( trec, "FEED_ID" ) ;
2077 RecordFieldPtr<Int> spwRF( trec, "SPECTRAL_WINDOW_ID" ) ;
2078 RecordFieldPtr<Double> timeRF( trec, "TIME" ) ;
2079 RecordFieldPtr<Double> intervalRF( trec, "INTERVAL" ) ;
2080 RecordFieldPtr< Array<Float> > tsysRF( trec, "TSYS" ) ;
2081 RecordFieldPtr< Array<Float> > tcalRF( trec, "TCAL" ) ;
2082 RecordFieldPtr< Array<Float> > tsysspRF ;
2083 RecordFieldPtr< Array<Float> > tcalspRF ;
2084 if ( tsysSpec_ )
2085 tsysspRF.attachToRecord( trec, "TSYS_SPECTRUM" ) ;
2086 if ( tcalSpec_ )
2087 tcalspRF.attachToRecord( trec, "TCAL_SPECTRUM" ) ;
[1977]2088
[2016]2089 // ANTENNA_ID is always 0
2090 *antennaRF = 0 ;
[1977]2091
[2016]2092 Table sortedstt = stt.sort( "ID" ) ;
2093 ROArrayColumn<Float> tcalCol( sortedstt, "TCAL" ) ;
2094 ROTableColumn idCol( sortedstt, "ID" ) ;
[2018]2095 ROArrayColumn<Float> tsysCol( tab, "TSYS" ) ;
2096 ROTableColumn tcalidCol( tab, "TCAL_ID" ) ;
2097 ROTableColumn timeCol( tab, "TIME" ) ;
2098 ROTableColumn intervalCol( tab, "INTERVAL" ) ;
2099 ROTableColumn beamnoCol( tab, "BEAMNO" ) ;
2100 ROTableColumn ifnoCol( tab, "IFNO" ) ;
[2291]2101 map< Int,Vector<uInt> >::iterator itr0 = idRec.begin() ;
2102 map< Int,Vector<uInt> >::iterator itr1 = rowRec.begin() ;
[2016]2103 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
[2258]2104// double t1 = mathutil::gettimeofday_sec() ;
[2291]2105 Vector<uInt> ids = itr0->second ;
2106 itr0++ ;
[2026]2107// os_ << "ids = " << ids << LogIO::POST ;
[2016]2108 uInt npol = ids.size() ;
[2291]2109 Vector<uInt> rows = itr1->second ;
2110 itr1++ ;
[2026]2111// os_ << "rows = " << rows << LogIO::POST ;
[2018]2112 Vector<Double> atime( rows.nelements() ) ;
2113 Vector<Double> ainterval( rows.nelements() ) ;
2114 Vector<uInt> atcalid( rows.nelements() ) ;
2115 for( uInt jrow = 0 ; jrow < rows.nelements() ; jrow++ ) {
2116 atime[jrow] = (Double)timeCol.asdouble( rows[jrow] ) ;
2117 ainterval[jrow] = (Double)intervalCol.asdouble( rows[jrow] ) ;
2118 atcalid[jrow] = tcalidCol.asuInt( rows[jrow] ) ;
2119 }
2120 Vector<Float> dummy = tsysCol( rows[0] ) ;
2121 Matrix<Float> tsys( npol,dummy.nelements() ) ;
2122 tsys.row( 0 ) = dummy ;
2123 for ( uInt jrow = 1 ; jrow < npol ; jrow++ )
2124 tsys.row( jrow ) = tsysCol( rows[jrow] ) ;
[1977]2125
[2016]2126 // FEED_ID
[2018]2127 *feedRF = beamnoCol.asuInt( rows[0] ) ;
[1977]2128
[2016]2129 // SPECTRAL_WINDOW_ID
[2018]2130 *spwRF = ifnoCol.asuInt( rows[0] ) ;
[1977]2131
[2016]2132 // TIME and INTERVAL
[2018]2133 getValidTimeRange( midTime, interval, atime, ainterval ) ;
[2016]2134 *timeRF = midTime ;
2135 *intervalRF = interval ;
2136
2137 // TCAL and TSYS
2138 Matrix<Float> tcal ;
2139 Table t ;
2140 if ( idCol.asuInt( ids[0] ) == ids[0] ) {
[2026]2141// os_ << "sorted at irow=" << irow << " ids[0]=" << ids[0] << LogIO::POST ;
[2016]2142 Vector<Float> dummyC = tcalCol( ids[0] ) ;
2143 tcal.resize( npol, dummyC.size() ) ;
2144 tcal.row( 0 ) = dummyC ;
[1977]2145 }
[2016]2146 else {
[2026]2147// os_ << "NOT sorted at irow=" << irow << " ids[0]=" << ids[0] << LogIO::POST ;
[2244]2148 t = stt( stt.col("ID") == ids[0], 1 ) ;
[2016]2149 Vector<Float> dummyC = tcalCol( 0 ) ;
[2244]2150 tcal.resize( npol, dummyC.size(), True ) ;
[2016]2151 tcal.row( 0 ) = dummyC ;
2152 }
2153 if ( npol == 2 ) {
2154 if ( idCol.asuInt( ids[1] ) == ids[1] ) {
[2026]2155// os_ << "sorted at irow=" << irow << " ids[1]=" << ids[1] << LogIO::POST ;
[2016]2156 tcal.row( 1 ) = tcalCol( ids[1] ) ;
2157 }
2158 else {
[2026]2159// os_ << "NOT sorted at irow=" << irow << " ids[1]=" << ids[1] << LogIO::POST ;
[2244]2160 t = stt( stt.col("ID") == ids[1], 1 ) ;
[2016]2161 tcalCol.attach( t, "TCAL" ) ;
[2244]2162 tcal.row( 1 ) = tcalCol( 0 ) ;
[2016]2163 }
2164 }
2165 else if ( npol == 3 ) {
2166 if ( idCol.asuInt( ids[2] ) == ids[2] )
2167 tcal.row( 1 ) = tcalCol( ids[2] ) ;
2168 else {
[2244]2169 t = stt( stt.col("ID") == ids[2], 1 ) ;
[2016]2170 tcalCol.attach( t, "TCAL" ) ;
2171 tcal.row( 1 ) = tcalCol( 0 ) ;
2172 }
2173 if ( idCol.asuInt( ids[1] ) == ids[1] )
2174 tcal.row( 2 ) = tcalCol( ids[1] ) ;
2175 else {
[2244]2176 t = stt( stt.col("ID") == ids[1], 1 ) ;
[2016]2177 tcalCol.attach( t, "TCAL" ) ;
2178 tcal.row( 2 ) = tcalCol( 0 ) ;
2179 }
2180 }
2181 else if ( npol == 4 ) {
2182 if ( idCol.asuInt( ids[2] ) == ids[2] )
2183 tcal.row( 1 ) = tcalCol( ids[2] ) ;
2184 else {
[2244]2185 t = stt( stt.col("ID") == ids[2], 1 ) ;
[2016]2186 tcalCol.attach( t, "TCAL" ) ;
2187 tcal.row( 1 ) = tcalCol( 0 ) ;
2188 }
2189 if ( idCol.asuInt( ids[3] ) == ids[3] )
2190 tcal.row( 2 ) = tcalCol( ids[3] ) ;
2191 else {
[2244]2192 t = stt( stt.col("ID") == ids[3], 1 ) ;
[2016]2193 tcalCol.attach( t, "TCAL" ) ;
2194 tcal.row( 2 ) = tcalCol( 0 ) ;
2195 }
2196 if ( idCol.asuInt( ids[1] ) == ids[1] )
2197 tcal.row( 2 ) = tcalCol( ids[1] ) ;
2198 else {
[2244]2199 t = stt( stt.col("ID") == ids[1], 1 ) ;
[2016]2200 tcalCol.attach( t, "TCAL" ) ;
2201 tcal.row( 3 ) = tcalCol( 0 ) ;
2202 }
2203 }
2204 if ( tcalSpec_ ) {
2205 // put TCAL_SPECTRUM
[2019]2206 tcalspRF.define( tcal ) ;
[2016]2207 // set TCAL (mean of TCAL_SPECTRUM)
2208 Matrix<Float> tcalMean( npol, 1 ) ;
2209 for ( uInt iid = 0 ; iid < npol ; iid++ ) {
2210 tcalMean( iid, 0 ) = mean( tcal.row(iid) ) ;
2211 }
2212 // put TCAL
[2291]2213 tcalRF.define( tcalMean ) ;
[2016]2214 }
2215 else {
2216 // put TCAL
[2291]2217 tcalRF.define( tcal ) ;
[2016]2218 }
[1977]2219
[2016]2220 if ( tsysSpec_ ) {
2221 // put TSYS_SPECTRUM
[2019]2222 tsysspRF.define( tsys ) ;
[2016]2223 // set TSYS (mean of TSYS_SPECTRUM)
2224 Matrix<Float> tsysMean( npol, 1 ) ;
2225 for ( uInt iid = 0 ; iid < npol ; iid++ ) {
2226 tsysMean( iid, 0 ) = mean( tsys.row(iid) ) ;
2227 }
2228 // put TSYS
[2291]2229 tsysRF.define( tsysMean ) ;
[2016]2230 }
2231 else {
2232 // put TSYS
[2291]2233 tsysRF.define( tsys ) ;
[2016]2234 }
[1977]2235
[2016]2236 // add row
2237 mssc.addRow( 1, True ) ;
2238 row.put( mssc.nrow()-1 ) ;
2239
[2258]2240// double t2 = mathutil::gettimeofday_sec() ;
[2026]2241// os_ << irow << "th loop elapsed time = " << t2-t1 << "sec" << LogIO::POST ;
[1977]2242 }
2243
[2291]2244 //double endSec = mathutil::gettimeofday_sec() ;
2245 //os_ << "end MSWriter::fillSysCal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
[1977]2246}
2247
[2016]2248void MSWriter::getValidTimeRange( Double &me, Double &interval, Table &tab )
[1974]2249{
[2258]2250// double startSec = mathutil::gettimeofday_sec() ;
[2026]2251// os_ << "start MSWriter::getVaridTimeRange() startSec=" << startSec << LogIO::POST ;
[2016]2252
2253 // sort table
2254 //Table stab = tab.sort( "TIME" ) ;
2255
[1975]2256 ROScalarColumn<Double> timeCol( tab, "TIME" ) ;
2257 Vector<Double> timeArr = timeCol.getColumn() ;
2258 Double minTime ;
2259 Double maxTime ;
2260 minMax( minTime, maxTime, timeArr ) ;
[2016]2261 Double midTime = 0.5 * ( minTime + maxTime ) * 86400.0 ;
2262 // unit for TIME
2263 // Scantable: "d"
2264 // MS: "s"
2265 me = midTime ;
2266 interval = ( maxTime - minTime ) * 86400.0 ;
2267
[2258]2268// double endSec = mathutil::gettimeofday_sec() ;
[2026]2269// os_ << "end MSWriter::getValidTimeRange() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
[1975]2270}
[1974]2271
[2018]2272void MSWriter::getValidTimeRange( Double &me, Double &interval, Vector<Double> &atime, Vector<Double> &ainterval )
2273{
[2258]2274// double startSec = mathutil::gettimeofday_sec() ;
[2026]2275// os_ << "start MSWriter::getVaridTimeRange() startSec=" << startSec << LogIO::POST ;
[2018]2276
2277 // sort table
2278 //Table stab = tab.sort( "TIME" ) ;
2279
2280 Double minTime ;
2281 Double maxTime ;
2282 minMax( minTime, maxTime, atime ) ;
2283 Double midTime = 0.5 * ( minTime + maxTime ) * 86400.0 ;
2284 // unit for TIME
2285 // Scantable: "d"
2286 // MS: "s"
2287 me = midTime ;
[2020]2288 interval = ( maxTime - minTime ) * 86400.0 + mean( ainterval ) ;
[2018]2289
[2258]2290// double endSec = mathutil::gettimeofday_sec() ;
[2026]2291// os_ << "end MSWriter::getValidTimeRange() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
[2018]2292}
2293
[1974]2294}
Note: See TracBrowser for help on using the repository browser.