source: trunk/src/MSFiller.cpp@ 2574

Last change on this file since 2574 was 2554, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: Yes CSV-1829

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

Modified the code to be able to handle INCREMENT for spectral axis correctly
regardless of sign of MS/SPECTRA_WINDOW/CHAN_WIDTH frequency.


File size: 65.8 KB
Line 
1
2//
3// C++ Interface: MSFiller
4//
5// Description:
6//
7// This class is specific filler for MS format
8// New version that is implemented using TableVisitor instead of TableIterator
9//
10// Takeshi Nakazato <takeshi.nakazato@nao.ac.jp>, (C) 2011
11//
12// Copyright: See COPYING file that comes with this distribution
13//
14//
15
16#include <assert.h>
17#include <iostream>
18#include <map>
19
20#include <tables/Tables/ExprNode.h>
21#include <tables/Tables/TableIter.h>
22#include <tables/Tables/TableColumn.h>
23#include <tables/Tables/ScalarColumn.h>
24#include <tables/Tables/ArrayColumn.h>
25#include <tables/Tables/TableParse.h>
26#include <tables/Tables/TableRow.h>
27
28#include <casa/Containers/RecordField.h>
29#include <casa/Logging/LogIO.h>
30#include <casa/Arrays/Slicer.h>
31#include <casa/Quanta/MVTime.h>
32#include <casa/OS/Path.h>
33
34#include <measures/Measures/Stokes.h>
35#include <measures/Measures/MEpoch.h>
36#include <measures/Measures/MCEpoch.h>
37#include <measures/Measures/MFrequency.h>
38#include <measures/Measures/MCFrequency.h>
39#include <measures/Measures/MPosition.h>
40#include <measures/Measures/MCPosition.h>
41#include <measures/Measures/MDirection.h>
42#include <measures/Measures/MCDirection.h>
43#include <measures/Measures/MeasConvert.h>
44#include <measures/TableMeasures/ScalarMeasColumn.h>
45#include <measures/TableMeasures/ArrayMeasColumn.h>
46#include <measures/TableMeasures/ScalarQuantColumn.h>
47#include <measures/TableMeasures/ArrayQuantColumn.h>
48
49#include <ms/MeasurementSets/MSAntennaIndex.h>
50
51#include <atnf/PKSIO/SrcType.h>
52
53#include "MSFiller.h"
54#include "STHeader.h"
55
56#include "MathUtils.h"
57
58using namespace casa ;
59using namespace std ;
60
61namespace asap {
62
63class BaseMSFillerVisitor: public TableVisitor {
64 uInt lastRecordNo ;
65 Int lastObservationId ;
66 Int lastFeedId ;
67 Int lastFieldId ;
68 Int lastDataDescId ;
69 Int lastScanNo ;
70 Int lastStateId ;
71 Double lastTime ;
72protected:
73 const Table &table;
74 uInt count;
75public:
76 BaseMSFillerVisitor(const Table &table)
77 : table(table)
78 {
79 count = 0;
80 }
81
82 virtual void enterObservationId(const uInt recordNo, Int columnValue) { }
83 virtual void leaveObservationId(const uInt recordNo, Int columnValue) { }
84 virtual void enterFeedId(const uInt recordNo, Int columnValue) { }
85 virtual void leaveFeedId(const uInt recordNo, Int columnValue) { }
86 virtual void enterFieldId(const uInt recordNo, Int columnValue) { }
87 virtual void leaveFieldId(const uInt recordNo, Int columnValue) { }
88 virtual void enterDataDescId(const uInt recordNo, Int columnValue) { }
89 virtual void leaveDataDescId(const uInt recordNo, Int columnValue) { }
90 virtual void enterScanNo(const uInt recordNo, Int columnValue) { }
91 virtual void leaveScanNo(const uInt recordNo, Int columnValue) { }
92 virtual void enterStateId(const uInt recordNo, Int columnValue) { }
93 virtual void leaveStateId(const uInt recordNo, Int columnValue) { }
94 virtual void enterTime(const uInt recordNo, Double columnValue) { }
95 virtual void leaveTime(const uInt recordNo, Double columnValue) { }
96
97 virtual Bool visitRecord(const uInt recordNo,
98 const Int ObservationId,
99 const Int feedId,
100 const Int fieldId,
101 const Int dataDescId,
102 const Int scanNo,
103 const Int stateId,
104 const Double time) { return True ; }
105
106 virtual Bool visit(Bool isFirst, const uInt recordNo,
107 const uInt nCols, void const *const colValues[]) {
108 Int observationId, feedId, fieldId, dataDescId, scanNo, stateId;
109 Double time;
110 { // prologue
111 uInt i = 0;
112 {
113 const Int *col = (const Int *)colValues[i++];
114 observationId = col[recordNo];
115 }
116 {
117 const Int *col = (const Int *)colValues[i++];
118 feedId = col[recordNo];
119 }
120 {
121 const Int *col = (const Int *)colValues[i++];
122 fieldId = col[recordNo];
123 }
124 {
125 const Int *col = (const Int *)colValues[i++];
126 dataDescId = col[recordNo];
127 }
128 {
129 const Int *col = (const Int *)colValues[i++];
130 scanNo = col[recordNo];
131 }
132 {
133 const Int *col = (const Int *)colValues[i++];
134 stateId = col[recordNo];
135 }
136 {
137 const Double *col = (const Double *)colValues[i++];
138 time = col[recordNo];
139 }
140 assert(nCols == i);
141 }
142
143 if (isFirst) {
144 enterObservationId(recordNo, observationId);
145 enterFeedId(recordNo, feedId);
146 enterFieldId(recordNo, fieldId);
147 enterDataDescId(recordNo, dataDescId);
148 enterScanNo(recordNo, scanNo);
149 enterStateId(recordNo, stateId);
150 enterTime(recordNo, time);
151 } else {
152 if (lastObservationId != observationId) {
153 leaveTime(lastRecordNo, lastTime);
154 leaveStateId(lastRecordNo, lastStateId);
155 leaveScanNo(lastRecordNo, lastScanNo);
156 leaveDataDescId(lastRecordNo, lastDataDescId);
157 leaveFieldId(lastRecordNo, lastFieldId);
158 leaveFeedId(lastRecordNo, lastFeedId);
159 leaveObservationId(lastRecordNo, lastObservationId);
160
161 enterObservationId(recordNo, observationId);
162 enterFeedId(recordNo, feedId);
163 enterFieldId(recordNo, fieldId);
164 enterDataDescId(recordNo, dataDescId);
165 enterScanNo(recordNo, scanNo);
166 enterStateId(recordNo, stateId);
167 enterTime(recordNo, time);
168 } else if (lastFeedId != feedId) {
169 leaveTime(lastRecordNo, lastTime);
170 leaveStateId(lastRecordNo, lastStateId);
171 leaveScanNo(lastRecordNo, lastScanNo);
172 leaveDataDescId(lastRecordNo, lastDataDescId);
173 leaveFieldId(lastRecordNo, lastFieldId);
174 leaveFeedId(lastRecordNo, lastFeedId);
175
176 enterFeedId(recordNo, feedId);
177 enterFieldId(recordNo, fieldId);
178 enterDataDescId(recordNo, dataDescId);
179 enterScanNo(recordNo, scanNo);
180 enterStateId(recordNo, stateId);
181 enterTime(recordNo, time);
182 } else if (lastFieldId != fieldId) {
183 leaveTime(lastRecordNo, lastTime);
184 leaveStateId(lastRecordNo, lastStateId);
185 leaveScanNo(lastRecordNo, lastScanNo);
186 leaveDataDescId(lastRecordNo, lastDataDescId);
187 leaveFieldId(lastRecordNo, lastFieldId);
188
189 enterFieldId(recordNo, fieldId);
190 enterDataDescId(recordNo, dataDescId);
191 enterScanNo(recordNo, scanNo);
192 enterStateId(recordNo, stateId);
193 enterTime(recordNo, time);
194 } else if (lastDataDescId != dataDescId) {
195 leaveTime(lastRecordNo, lastTime);
196 leaveStateId(lastRecordNo, lastStateId);
197 leaveScanNo(lastRecordNo, lastScanNo);
198 leaveDataDescId(lastRecordNo, lastDataDescId);
199
200 enterDataDescId(recordNo, dataDescId);
201 enterScanNo(recordNo, scanNo);
202 enterStateId(recordNo, stateId);
203 enterTime(recordNo, time);
204 } else if (lastScanNo != scanNo) {
205 leaveTime(lastRecordNo, lastTime);
206 leaveStateId(lastRecordNo, lastStateId);
207 leaveScanNo(lastRecordNo, lastScanNo);
208
209 enterScanNo(recordNo, scanNo);
210 enterStateId(recordNo, stateId);
211 enterTime(recordNo, time);
212 } else if (lastStateId != stateId) {
213 leaveTime(lastRecordNo, lastTime);
214 leaveStateId(lastRecordNo, lastStateId);
215
216 enterStateId(recordNo, stateId);
217 enterTime(recordNo, time);
218 } else if (lastTime != time) {
219 leaveTime(lastRecordNo, lastTime);
220 enterTime(recordNo, time);
221 }
222 }
223 count++;
224 Bool result = visitRecord(recordNo, observationId, feedId, fieldId, dataDescId,
225 scanNo, stateId, time);
226
227 { // epilogue
228 lastRecordNo = recordNo;
229
230 lastObservationId = observationId;
231 lastFeedId = feedId;
232 lastFieldId = fieldId;
233 lastDataDescId = dataDescId;
234 lastScanNo = scanNo;
235 lastStateId = stateId;
236 lastTime = time;
237 }
238 return result ;
239 }
240
241 virtual void finish() {
242 if (count > 0) {
243 leaveTime(lastRecordNo, lastTime);
244 leaveStateId(lastRecordNo, lastStateId);
245 leaveScanNo(lastRecordNo, lastScanNo);
246 leaveDataDescId(lastRecordNo, lastDataDescId);
247 leaveFieldId(lastRecordNo, lastFieldId);
248 leaveFeedId(lastRecordNo, lastFeedId);
249 leaveObservationId(lastRecordNo, lastObservationId);
250 }
251 }
252};
253
254class MSFillerVisitor: public BaseMSFillerVisitor, public MSFillerUtils {
255public:
256 MSFillerVisitor(const Table &from, Scantable &to)
257 : BaseMSFillerVisitor(from),
258 scantable( to )
259 {
260 antennaId = 0 ;
261 rowidx = 0 ;
262 tablerow = TableRow( scantable.table() ) ;
263 feedEntry = Vector<Int>( 64, -1 ) ;
264 nbeam = 0 ;
265 ifmap.clear() ;
266 const TableDesc &desc = table.tableDesc() ;
267 if ( desc.isColumn( "DATA" ) )
268 dataColumnName = "DATA" ;
269 else if ( desc.isColumn( "FLOAT_DATA" ) )
270 dataColumnName = "FLOAT_DATA" ;
271 getpt = False ;
272 isWeather = False ;
273 isSysCal = False ;
274 isTcal = False ;
275 cycleNo = 0 ;
276 numSysCalRow = 0 ;
277 header = scantable.getHeader() ;
278 fluxUnit( header.fluxunit ) ;
279
280 // MS subtables
281 const TableRecord &hdr = table.keywordSet();
282 obstab = hdr.asTable( "OBSERVATION" ) ;
283 spwtab = hdr.asTable( "SPECTRAL_WINDOW" ) ;
284 statetab = hdr.asTable( "STATE" ) ;
285 ddtab = hdr.asTable( "DATA_DESCRIPTION" ) ;
286 poltab = hdr.asTable( "POLARIZATION" ) ;
287 fieldtab = hdr.asTable( "FIELD" ) ;
288 anttab = hdr.asTable( "ANTENNA" ) ;
289 if ( hdr.isDefined( "SYSCAL" ) )
290 sctab = hdr.asTable( "SYSCAL" ) ;
291 if ( hdr.isDefined( "SOURCE" ) )
292 srctab = hdr.asTable( "SOURCE" ) ;
293
294 // attach to columns
295 // MS MAIN
296 intervalCol.attach( table, "INTERVAL" ) ;
297 flagRowCol.attach( table, "FLAG_ROW" ) ;
298 flagCol.attach( table, "FLAG" ) ;
299 if ( dataColumnName.compare( "DATA" ) == 0 )
300 dataCol.attach( table, dataColumnName ) ;
301 else
302 floatDataCol.attach( table, dataColumnName ) ;
303
304 // set dummy epoch
305 mf.set( currentTime ) ;
306
307 // initialize dirType
308 dirType = MDirection::N_Types ;
309
310 //
311 // add rows to scantable
312 //
313 // number of polarization is up to 4
314 uInt addrow = table.nrow() * maxNumPol() ;
315 scantable.table().addRow( addrow ) ;
316
317 // attach to columns
318 // Scantable MAIN
319 TableRecord &r = tablerow.record() ;
320 timeRF.attachToRecord( r, "TIME" ) ;
321 intervalRF.attachToRecord( r, "INTERVAL" ) ;
322 directionRF.attachToRecord( r, "DIRECTION" ) ;
323 azimuthRF.attachToRecord( r, "AZIMUTH" ) ;
324 elevationRF.attachToRecord( r, "ELEVATION" ) ;
325 scanRateRF.attachToRecord( r, "SCANRATE" ) ;
326 weatherIdRF.attachToRecord( r, "WEATHER_ID" ) ;
327 cycleNoRF.attachToRecord( r, "CYCLENO" ) ;
328 flagRowRF.attachToRecord( r, "FLAGROW" ) ;
329 polNoRF.attachToRecord( r, "POLNO" ) ;
330 tcalIdRF.attachToRecord( r, "TCAL_ID" ) ;
331 spectraRF.attachToRecord( r, "SPECTRA" ) ;
332 flagtraRF.attachToRecord( r, "FLAGTRA" ) ;
333 tsysRF.attachToRecord( r, "TSYS" ) ;
334 beamNoRF.attachToRecord( r, "BEAMNO" ) ;
335 ifNoRF.attachToRecord( r, "IFNO" ) ;
336 freqIdRF.attachToRecord( r, "FREQ_ID" ) ;
337 moleculeIdRF.attachToRecord( r, "MOLECULE_ID" ) ;
338 sourceNameRF.attachToRecord( r, "SRCNAME" ) ;
339 sourceProperMotionRF.attachToRecord( r, "SRCPROPERMOTION" ) ;
340 sourceDirectionRF.attachToRecord( r, "SRCDIRECTION" ) ;
341 sourceVelocityRF.attachToRecord( r, "SRCVELOCITY" ) ;
342 focusIdRF.attachToRecord( r, "FOCUS_ID" ) ;
343 fieldNameRF.attachToRecord( r, "FIELDNAME" ) ;
344 sourceTypeRF.attachToRecord( r, "SRCTYPE" ) ;
345 scanNoRF.attachToRecord( r, "SCANNO" ) ;
346
347 // put values
348 RecordFieldPtr<Int> refBeamNoRF( r, "REFBEAMNO" ) ;
349 *refBeamNoRF = -1 ;
350 RecordFieldPtr<Int> fitIdRF( r, "FIT_ID" ) ;
351 *fitIdRF = -1 ;
352 RecordFieldPtr<Float> opacityRF( r, "OPACITY" ) ;
353 *opacityRF = 0.0 ;
354 }
355
356 virtual void enterObservationId(const uInt recordNo, Int columnValue) {
357 //printf("%u: ObservationId: %d\n", recordNo, columnValue);
358 // update header
359 if ( header.observer.empty() )
360 getScalar( String("OBSERVER"), (uInt)columnValue, obstab, header.observer ) ;
361 if ( header.project.empty() )
362 getScalar( "PROJECT", (uInt)columnValue, obstab, header.project ) ;
363 if ( header.utc == 0.0 ) {
364 Vector<MEpoch> amp ;
365 getArrayMeas( "TIME_RANGE", (uInt)columnValue, obstab, amp ) ;
366 header.utc = amp[0].get( "d" ).getValue() ;
367 }
368 if ( header.antennaname.empty() )
369 getScalar( "TELESCOPE_NAME", (uInt)columnValue, obstab, header.antennaname ) ;
370 }
371 virtual void leaveObservationId(const uInt recordNo, Int columnValue) {
372 // update header
373 header.nbeam = max( header.nbeam, (Int)nbeam ) ;
374
375 nbeam = 0 ;
376 feedEntry = -1 ;
377 }
378 virtual void enterFeedId(const uInt recordNo, Int columnValue) {
379 //printf("%u: FeedId: %d\n", recordNo, columnValue);
380
381 // update feed entry
382 if ( allNE( feedEntry, columnValue ) ) {
383 feedEntry[nbeam] = columnValue ;
384 nbeam++ ;
385 }
386
387 // put values
388 *beamNoRF = (uInt)columnValue ;
389 *focusIdRF = (uInt)0 ;
390 }
391 virtual void leaveFeedId(const uInt recordNo, Int columnValue) {
392 Int nelem = (Int)feedEntry.nelements() ;
393 if ( nbeam > nelem ) {
394 feedEntry.resize( nelem+64, True ) ;
395 Slicer slice( IPosition( 1, nelem ), IPosition( 1, feedEntry.nelements()-1 ) ) ;
396 feedEntry( slice ) = -1 ;
397 }
398 }
399 virtual void enterFieldId(const uInt recordNo, Int columnValue) {
400 //printf("%u: FieldId: %d\n", recordNo, columnValue);
401 // update sourceId and fieldName
402 getScalar( "SOURCE_ID", (uInt)columnValue, fieldtab, sourceId ) ;
403 String fieldName ;
404 getScalar( "NAME", (uInt)columnValue, fieldtab, fieldName ) ;
405 fieldName += "__" + String::toString( columnValue ) ;
406
407 // put values
408 *fieldNameRF = fieldName ;
409 }
410 virtual void leaveFieldId(const uInt recordNo, Int columnValue) {
411 sourceId = -1 ;
412 }
413 virtual void enterDataDescId(const uInt recordNo, Int columnValue) {
414 //printf("%u: DataDescId: %d\n", recordNo, columnValue);
415 // update polarization and spectral window ids
416 getScalar( "POLARIZATION_ID", (uInt)columnValue, ddtab, polId ) ;
417 getScalar( "SPECTRAL_WINDOW_ID", (uInt)columnValue, ddtab, spwId ) ;
418
419 // polarization setup
420 getScalar( "NUM_CORR", (uInt)polId, poltab, npol ) ;
421 Vector<Int> corrtype ;
422 getArray( "CORR_TYPE", (uInt)polId, poltab, corrtype ) ;
423 polnos = getPolNos( corrtype ) ;
424
425 // process SOURCE table
426 String sourceName ;
427 Vector<Double> sourcePM, restFreqs, sysVels ;
428 Vector<String> transition ;
429 processSource( sourceId, spwId, sourceName, sourceDir, sourcePM,
430 restFreqs, transition, sysVels ) ;
431
432 // spectral setup
433 uInt freqId ;
434 Double reffreq, bandwidth ;
435 String freqref ;
436 getScalar( "NUM_CHAN", (uInt)spwId, spwtab, nchan ) ;
437 Bool iswvr = (Bool)(nchan == 4) ;
438 map<Int,uInt>::iterator iter = ifmap.find( spwId ) ;
439 if ( iter == ifmap.end() ) {
440 MEpoch me ;
441 getScalarMeas( "TIME", recordNo, table, me ) ;
442 spectralSetup( spwId, me, antpos, sourceDir,
443 freqId, nchan,
444 freqref, reffreq, bandwidth ) ;
445 ifmap.insert( pair<Int,uInt>(spwId,freqId) ) ;
446 }
447 else {
448 freqId = iter->second ;
449 }
450 sp.resize( npol, nchan ) ;
451 fl.resize( npol, nchan ) ;
452
453
454 // molecular setup
455 STMolecules mtab = scantable.molecules() ;
456 uInt molId = mtab.addEntry( restFreqs, transition, transition ) ;
457
458 // process SYSCAL table
459 if ( isSysCal )
460 processSysCal( spwId ) ;
461
462 // update header
463 if ( !iswvr ) {
464 header.nchan = max( header.nchan, nchan ) ;
465 header.bandwidth = max( header.bandwidth, bandwidth ) ;
466 if ( header.reffreq == -1.0 )
467 header.reffreq = reffreq ;
468 header.npol = max( header.npol, npol ) ;
469 if ( header.poltype.empty() )
470 header.poltype = getPolType( corrtype[0] ) ;
471 if ( header.freqref.empty() )
472 header.freqref = freqref ;
473 }
474
475 // put values
476 *ifNoRF = (uInt)spwId ;
477 *freqIdRF = freqId ;
478 *moleculeIdRF = molId ;
479 *sourceNameRF = sourceName ;
480 sourceProperMotionRF.define( sourcePM ) ;
481 Vector<Double> srcD = sourceDir.getAngle().getValue( "rad" ) ;
482 sourceDirectionRF.define( srcD ) ;
483 if ( !sysVels.empty() )
484 *sourceVelocityRF = sysVels[0] ;
485 else {
486 *sourceVelocityRF = (Double)0.0 ;
487 }
488 }
489 virtual void leaveDataDescId(const uInt recordNo, Int columnValue) {
490 npol = 0 ;
491 nchan = 0 ;
492 numSysCalRow = 0 ;
493 }
494 virtual void enterScanNo(const uInt recordNo, Int columnValue) {
495 //printf("%u: ScanNo: %d\n", recordNo, columnValue);
496 // put value
497 // scan number is 1-based in MS while 0-based in Scantable
498 *scanNoRF = (uInt)columnValue - 1 ;
499 }
500 virtual void leaveScanNo(const uInt recordNo, Int columnValue) {
501 cycleNo = 0 ;
502 }
503 virtual void enterStateId(const uInt recordNo, Int columnValue) {
504 //printf("%u: StateId: %d\n", recordNo, columnValue);
505 // SRCTYPE
506 Int srcType = getSrcType( columnValue ) ;
507
508 // update header
509 if ( header.obstype.empty() )
510 getScalar( "OBS_MODE", (uInt)columnValue, statetab, header.obstype ) ;
511
512 // put value
513 *sourceTypeRF = srcType ;
514 }
515 virtual void leaveStateId(const uInt recordNo, Int columnValue) { }
516 virtual void enterTime(const uInt recordNo, Double columnValue) {
517 //printf("%u: Time: %f\n", recordNo, columnValue);
518 currentTime = MEpoch( Quantity( columnValue, "s" ), MEpoch::UTC ) ;
519
520 // DIRECTION, AZEL, and SCANRATE
521 Vector<Double> direction, azel ;
522 Vector<Double> scanrate( 2, 0.0 ) ;
523 if ( getpt )
524 getDirection( direction, azel, scanrate ) ;
525 else
526 getSourceDirection( direction, azel, scanrate ) ;
527
528 // INTERVAL
529 Double interval = intervalCol.asdouble( recordNo ) ;
530
531 // WEATHER_ID
532 uInt wid = 0 ;
533 if ( isWeather )
534 wid = getWeatherId() ;
535
536 // put value
537 Double t = currentTime.get( "d" ).getValue() ;
538 *timeRF = t ;
539 *intervalRF = interval ;
540 directionRF.define( direction ) ;
541 *azimuthRF = (Float)azel[0] ;
542 *elevationRF = (Float)azel[1] ;
543 scanRateRF.define( scanrate ) ;
544 *weatherIdRF = wid ;
545 }
546 virtual void leaveTime(const uInt recordNo, Double columnValue) { }
547 virtual Bool visitRecord(const uInt recordNo,
548 const Int observationId,
549 const Int feedId,
550 const Int fieldId,
551 const Int dataDescId,
552 const Int scanNo,
553 const Int stateId,
554 const Double time)
555 {
556 //printf("%u: %d, %d, %d, %d, %d, %d, %f\n", recordNo,
557 //observationId, feedId, fieldId, dataDescId, scanNo, stateId, time);
558
559 // SPECTRA and FLAGTRA
560 //Matrix<Float> sp;
561 //Matrix<uChar> fl;
562 spectraAndFlagtra( recordNo, sp, fl ) ;
563
564 // FLAGROW
565 Bool flr = flagRowCol.asBool( recordNo ) ;
566
567 // TSYS
568 Matrix<Float> tsys ;
569 uInt scIdx = getSysCalIndex() ;
570 if ( numSysCalRow > 0 ) {
571 tsys = sysCalTsysCol( syscalRow[scIdx] ) ;
572 }
573 else {
574 tsys.resize( npol, 1 ) ;
575 tsys = 1.0 ;
576 }
577
578 // TCAL_ID
579 Block<uInt> tcalids( npol, 0 ) ;
580 if ( numSysCalRow > 0 ) {
581 tcalids = getTcalId( syscalTime[scIdx] ) ;
582 }
583
584 // put value
585 *cycleNoRF = cycleNo ;
586 *flagRowRF = (uInt)flr ;
587
588 // for each polarization component
589 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
590 // put value depending on polarization component
591 *polNoRF = polnos[ipol] ;
592 *tcalIdRF = tcalids[ipol] ;
593 spectraRF.define( sp.row( ipol ) ) ;
594 flagtraRF.define( fl.row( ipol ) ) ;
595 tsysRF.define( tsys.row( ipol ) ) ;
596
597 // commit row
598 tablerow.put( rowidx ) ;
599 rowidx++ ;
600 }
601
602 // increment CYCLENO
603 cycleNo++ ;
604
605 return True ;
606 }
607 virtual void finish()
608 {
609 BaseMSFillerVisitor::finish();
610 //printf("Total: %u\n", count);
611 // remove redundant rows
612 //cout << "filled " << rowidx << " rows out of " << scantable.nrow() << " rows" << endl ;
613 if ( scantable.nrow() > rowidx ) {
614 uInt numRemove = scantable.nrow() - rowidx ;
615 //cout << "numRemove = " << numRemove << endl ;
616 Vector<uInt> rows( numRemove ) ;
617 indgen( rows, rowidx ) ;
618 scantable.table().removeRow( rows ) ;
619 }
620
621 // antenna name and station name
622 String antennaName ;
623 getScalar( "NAME", (uInt)antennaId, anttab, antennaName ) ;
624 String stationName ;
625 getScalar( "STATION", (uInt)antennaId, anttab, stationName ) ;
626
627 // update header
628 header.nif = ifmap.size() ;
629 header.antennaposition = antpos.get( "m" ).getValue() ;
630 if ( header.antennaname.empty() || header.antennaname == antennaName )
631 header.antennaname = antennaName ;
632 else
633 header.antennaname += "//" + antennaName ;
634 if ( !stationName.empty() && stationName != antennaName )
635 header.antennaname += "@" + stationName ;
636 if ( header.fluxunit.empty() || header.fluxunit == "CNTS" )
637 header.fluxunit = "K" ;
638 header.epoch = "UTC" ;
639 header.equinox = 2000.0 ;
640 if (header.freqref == "TOPO") {
641 header.freqref = "TOPOCENT";
642 } else if (header.freqref == "GEO") {
643 header.freqref = "GEOCENTR";
644 } else if (header.freqref == "BARY") {
645 header.freqref = "BARYCENT";
646 } else if (header.freqref == "GALACTO") {
647 header.freqref = "GALACTOC";
648 } else if (header.freqref == "LGROUP") {
649 header.freqref = "LOCALGRP";
650 } else if (header.freqref == "CMB") {
651 header.freqref = "CMBDIPOL";
652 } else if (header.freqref == "REST") {
653 header.freqref = "SOURCE";
654 }
655 scantable.setHeader( header ) ;
656 }
657 void setAntenna( Int id )
658 {
659 antennaId = id ;
660
661 Vector< Quantum<Double> > pos ;
662 getArrayQuant( "POSITION", (uInt)antennaId, anttab, pos ) ;
663 antpos = MPosition( MVPosition( pos ), MPosition::ITRF ) ;
664 mf.set( antpos ) ;
665 }
666 void setPointingTable( const Table &tab, String columnToUse="DIRECTION" )
667 {
668 // input POINTING table must be
669 // 1) selected by antenna
670 // 2) sorted by TIME
671 ROScalarColumn<Double> tcol( tab, "TIME" ) ;
672 ROArrayColumn<Double> dcol( tab, columnToUse ) ;
673 tcol.getColumn( pointingTime ) ;
674 dcol.getColumn( pointingDirection ) ;
675 const TableRecord &rec = dcol.keywordSet() ;
676 String pointingRef = rec.asRecord( "MEASINFO" ).asString( "Ref" ) ;
677 MDirection::getType( dirType, pointingRef ) ;
678 getpt = True ;
679 }
680 void setWeatherTime( const Vector<Double> &t, const Vector<Double> &it )
681 {
682 isWeather = True ;
683 weatherTime = t ;
684 weatherInterval = it ;
685 }
686 void setSysCalRecord( const Record &r )
687 //void setSysCalRecord( const map< String,Vector<uInt> > &r )
688 {
689 isSysCal = True ;
690 isTcal = True ;
691 syscalRecord = r ;
692 if ( syscalRecord.nfields() == 0 )
693 isTcal = False ;
694
695 const TableDesc &desc = sctab.tableDesc() ;
696 uInt nrow = sctab.nrow() ;
697 syscalRow.resize( nrow ) ;
698 syscalTime.resize( nrow ) ;
699 syscalInterval.resize( nrow ) ;
700 String tsysCol = "NONE" ;
701 Vector<String> tsysCols = stringToVector( "TSYS_SPECTRUM,TSYS" ) ;
702 for ( uInt i = 0 ; i < tsysCols.nelements() ; i++ ) {
703 if ( tsysCol == "NONE" && desc.isColumn( tsysCols[i] ) )
704 tsysCol = tsysCols[i] ;
705 }
706 sysCalTsysCol.attach( sctab, tsysCol ) ;
707 }
708 STHeader getHeader() { return header ; }
709 uInt getNumBeam() { return nbeam ; }
710 uInt getFilledRowNum() { return rowidx ; }
711private:
712 void fluxUnit( String &u )
713 {
714 ROTableColumn col( table, dataColumnName ) ;
715 const TableRecord &rec = col.keywordSet() ;
716 if ( rec.isDefined( "UNIT" ) )
717 u = rec.asString( "UNIT" ) ;
718 else if ( rec.isDefined( "QuantumUnits" ) )
719 u = rec.asString( "QuantumUnits" ) ;
720 if ( u.empty() )
721 u = "K" ;
722 }
723 void processSource( Int sourceId, Int spwId,
724 String &name, MDirection &dir, Vector<Double> &pm,
725 Vector<Double> &rf, Vector<String> &trans, Vector<Double> &vel )
726 {
727 // find row
728 uInt nrow = srctab.nrow() ;
729 Int idx = -1 ;
730 ROTableRow row( srctab ) ;
731 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
732 const TableRecord &r = row.get( irow ) ;
733 if ( r.asInt( "SOURCE_ID" ) == sourceId ) {
734 Int tmpSpwId = r.asInt( "SPECTRAL_WINDOW_ID" ) ;
735 if ( tmpSpwId == spwId || tmpSpwId == -1 ) {
736 idx = (Int)irow ;
737 break ;
738 }
739 }
740 }
741
742 // fill
743 Int numLines = 0 ;
744 if ( idx != -1 ) {
745 const TableRecord &r = row.get( idx ) ;
746 name = r.asString( "NAME" ) ;
747 getScalarMeas( "DIRECTION", idx, srctab, dir ) ;
748 pm = r.toArrayDouble( "PROPER_MOTION" ) ;
749 numLines = r.asInt( "NUM_LINES" ) ;
750 }
751 else {
752 name = "" ;
753 pm = Vector<Double>( 2, 0.0 ) ;
754 dir = MDirection( Quantum<Double>(0.0,Unit("rad")), Quantum<Double>(0.0,Unit("rad")) ) ;
755 }
756 if ( !getpt ) {
757 String ref = dir.getRefString() ;
758 MDirection::getType( dirType, ref ) ;
759 }
760
761 rf.resize( numLines ) ;
762 trans.resize( numLines ) ;
763 vel.resize( numLines ) ;
764 if ( numLines > 0 ) {
765 Block<Bool> isDefined = row.getDefined() ;
766 Vector<String> colNames = row.columnNames() ;
767 Vector<Int> indexes( 3, -1 ) ;
768 Vector<String> cols = stringToVector( "REST_FREQUENCY,TRANSITION,SYSVEL" ) ;
769 for ( uInt icol = 0 ; icol < colNames.nelements() ; icol++ ) {
770 if ( anyEQ( indexes, -1 ) ) {
771 for ( uInt jcol = 0 ; jcol < cols.nelements() ; jcol++ ) {
772 if ( colNames[icol] == cols[jcol] )
773 indexes[jcol] = icol ;
774 }
775 }
776 }
777 if ( indexes[0] != -1 && isDefined[indexes[0]] == True ) {
778 Vector< Quantum<Double> > qrf ;
779 getArrayQuant( "REST_FREQUENCY", idx, srctab, qrf ) ;
780 for ( int i = 0 ; i < numLines ; i++ )
781 rf[i] = qrf[i].getValue( "Hz" ) ;
782 }
783 if ( indexes[1] != -1 && isDefined[indexes[1]] == True ) {
784 getArray( "TRANSITION", idx, srctab, trans ) ;
785 }
786 if ( indexes[2] != -1 && isDefined[indexes[2]] == True ) {
787 Vector< Quantum<Double> > qsv ;
788 getArrayQuant( "SYSVEL", idx, srctab, qsv ) ;
789 for ( int i = 0 ; i < numLines ; i++ )
790 vel[i] = qsv[i].getValue( "m/s" ) ;
791 }
792 }
793 }
794 void spectralSetup( Int &spwId, MEpoch &me, MPosition &mp, MDirection &md,
795 uInt &freqId, Int &nchan,
796 String &freqref, Double &reffreq, Double &bandwidth )
797 {
798 // fill
799 Int measFreqRef ;
800 getScalar( "MEAS_FREQ_REF", spwId, spwtab, measFreqRef ) ;
801 MFrequency::Types freqRef = MFrequency::castType( measFreqRef ) ;
802 //freqref = MFrequency::showType( freqRef ) ;
803 freqref = "LSRK" ;
804 Quantum<Double> q ;
805 getScalarQuant( "TOTAL_BANDWIDTH", spwId, spwtab, q ) ;
806 bandwidth = q.getValue( "Hz" ) ;
807 getScalarQuant( "REF_FREQUENCY", spwId, spwtab, q ) ;
808 reffreq = q.getValue( "Hz" ) ;
809 Double refpix = 0.5 * ( (Double)nchan-1.0 ) ;
810 Int refchan = ( nchan - 1 ) / 2 ;
811 Bool even = (Bool)( nchan % 2 == 0 ) ;
812 Vector< Quantum<Double> > qa ;
813 getArrayQuant( "CHAN_WIDTH", spwId, spwtab, qa ) ;
814// Double increment = qa[refchan].getValue( "Hz" ) ;
815 Double increment = abs(qa[refchan].getValue( "Hz" )) ;
816 getArrayQuant( "CHAN_FREQ", spwId, spwtab, qa ) ;
817 if ( nchan == 1 ) {
818 Int netSideband ;
819 getScalar( "NET_SIDEBAND", spwId, spwtab, netSideband ) ;
820 if ( netSideband == 1 ) increment *= -1.0 ;
821 }
822 else {
823 if ( qa[0].getValue( "Hz" ) > qa[1].getValue( "Hz" ) )
824 increment *= -1.0 ;
825 }
826 Double refval = qa[refchan].getValue( "Hz" ) ;
827 if ( even )
828 refval = 0.5 * ( refval + qa[refchan+1].getValue( "Hz" ) ) ;
829 if ( freqRef != MFrequency::LSRK ) {
830 MeasFrame mframe( me, mp, md ) ;
831 MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mframe ) ) ;
832 refval = tolsr( Quantum<Double>( refval, "Hz" ) ).get( "Hz" ).getValue() ;
833 }
834
835 // add new row to FREQUENCIES
836 Table ftab = scantable.frequencies().table() ;
837 freqId = ftab.nrow() ;
838 ftab.addRow() ;
839 TableRow row( ftab ) ;
840 TableRecord &r = row.record() ;
841 RecordFieldPtr<uInt> idRF( r, "ID" ) ;
842 *idRF = freqId ;
843 RecordFieldPtr<Double> refpixRF( r, "REFPIX" ) ;
844 RecordFieldPtr<Double> refvalRF( r, "REFVAL" ) ;
845 RecordFieldPtr<Double> incrRF( r, "INCREMENT" ) ;
846 *refpixRF = refpix ;
847 *refvalRF = refval ;
848 *incrRF = increment ;
849 row.put( freqId ) ;
850 }
851 void spectraAndFlagtra( uInt recordNo, Matrix<Float> &sp, Matrix<uChar> &fl )
852 {
853 Matrix<Bool> b = flagCol( recordNo ) ;
854 if ( dataColumnName.compare( "FLOAT_DATA" ) == 0 ) {
855 sp = floatDataCol( recordNo ) ;
856 convertArray( fl, b ) ;
857 }
858 else {
859 Bool notyet = True ;
860 Matrix<Complex> c = dataCol( recordNo ) ;
861 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
862 if ( ( header.poltype == "linear" || header.poltype == "circular" )
863 && ( polnos[ipol] == 2 || polnos[ipol] == 3 ) ) {
864 if ( notyet ) {
865 Vector<Float> tmp = ComplexToReal( c.row( ipol ) ) ;
866 IPosition start( 1, 0 ) ;
867 IPosition end( 1, 2*nchan-1 ) ;
868 IPosition inc( 1, 2 ) ;
869 if ( polnos[ipol] == 2 ) {
870 sp.row( ipol ) = tmp( start, end, inc ) ;
871 Vector<Bool> br = b.row( ipol ) ;
872 Vector<uChar> flr = fl.row( ipol ) ;
873 convertArray( flr, br ) ;
874 start = IPosition( 1, 1 ) ;
875 Int jpol = ipol+1 ;
876 while( polnos[jpol] != 3 && jpol < npol )
877 jpol++ ;
878 sp.row( jpol ) = tmp( start, end, inc ) ;
879 flr.reference( fl.row( jpol ) ) ;
880 convertArray( flr, br ) ;
881 }
882 else if ( polnos[ipol] == 3 ) {
883 sp.row( ipol ) = sp.row( ipol ) * (Float)(-1.0) ;
884 Int jpol = ipol+1 ;
885 while( polnos[jpol] != 2 && jpol < npol )
886 jpol++ ;
887 Vector<Bool> br = b.row( ipol ) ;
888 Vector<uChar> flr = fl.row( jpol ) ;
889 sp.row( jpol ) = tmp( start, end, inc ) ;
890 convertArray( flr, br ) ;
891 start = IPosition( 1, 1 ) ;
892 sp.row( ipol ) = tmp( start, end, inc ) * (Float)(-1.0) ;
893 flr.reference( fl.row( ipol ) ) ;
894 convertArray( flr, br ) ;
895 }
896 notyet = False ;
897 }
898 }
899 else {
900 Vector<Float> tmp = ComplexToReal( c.row( ipol ) ) ;
901 IPosition start( 1, 0 ) ;
902 IPosition end( 1, 2*nchan-1 ) ;
903 IPosition inc( 1, 2 ) ;
904 sp.row( ipol ) = tmp( start, end, inc ) ;
905 Vector<Bool> br = b.row( ipol ) ;
906 Vector<uChar> flr = fl.row( ipol ) ;
907 convertArray( flr, br ) ;
908 }
909 }
910 }
911 }
912 uInt binarySearch( Vector<Double> &timeList, Double target )
913 {
914 Int low = 0 ;
915 Int high = timeList.nelements() ;
916 uInt idx = 0 ;
917
918 while ( low <= high ) {
919 idx = (Int)( 0.5 * ( low + high ) ) ;
920 Double t = timeList[idx] ;
921 if ( t < target )
922 low = idx + 1 ;
923 else if ( t > target )
924 high = idx - 1 ;
925 else {
926 return idx ;
927 }
928 }
929
930 idx = max( 0, min( low, high ) ) ;
931 return idx ;
932 }
933 void getDirection( Vector<Double> &dir, Vector<Double> &azel, Vector<Double> &srate )
934 {
935 // @todo At the moment, do binary search every time
936 // if this is bottleneck, frequency of binary search must be reduced
937 Double t = currentTime.get( "s" ).getValue() ;
938 uInt idx = min( binarySearch( pointingTime, t ), pointingTime.nelements()-1 ) ;
939 Matrix<Double> d ;
940 if ( pointingTime[idx] == t )
941 d = pointingDirection.xyPlane( idx ) ;
942 else if ( pointingTime[idx] < t ) {
943 if ( idx == pointingTime.nelements()-1 )
944 d = pointingDirection.xyPlane( idx ) ;
945 else
946 d = interp( pointingTime[idx], pointingTime[idx+1], t,
947 pointingDirection.xyPlane( idx ), pointingDirection.xyPlane( idx+1 ) ) ;
948 }
949 else {
950 if ( idx == 0 )
951 d = pointingDirection.xyPlane( idx ) ;
952 else
953 d = interp( pointingTime[idx-1], pointingTime[idx], t,
954 pointingDirection.xyPlane( idx-1 ), pointingDirection.xyPlane( idx ) ) ;
955 }
956 mf.resetEpoch( currentTime ) ;
957 Quantum< Vector<Double> > tmp( d.column( 0 ), Unit( "rad" ) ) ;
958 if ( dirType != MDirection::J2000 ) {
959 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
960 dir = toj2000( tmp ).getAngle( "rad" ).getValue() ;
961 }
962 else {
963 dir = d.column( 0 ) ;
964 }
965 if ( dirType != MDirection::AZELGEO ) {
966 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ;
967 //MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
968 azel = toazel( tmp ).getAngle( "rad" ).getValue() ;
969 }
970 else {
971 azel = d.column( 0 ) ;
972 }
973 if ( d.ncolumn() > 1 )
974 srate = d.column( 1 ) ;
975 }
976 void getSourceDirection( Vector<Double> &dir, Vector<Double> &azel, Vector<Double> &srate )
977 {
978 dir = sourceDir.getAngle( "rad" ).getValue() ;
979 mf.resetEpoch( currentTime ) ;
980 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ;
981 azel = toazel( Quantum< Vector<Double> >( dir, Unit("rad") ) ).getAngle( "rad" ).getValue() ;
982 if ( dirType != MDirection::J2000 ) {
983 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
984 dir = toj2000( Quantum< Vector<Double> >( dir, Unit("rad") ) ).getAngle( "rad" ).getValue() ;
985 }
986 }
987 String detectSeparator( String &s )
988 {
989 String tmp = s.substr( 0, s.find_first_of( "," ) ) ;
990 Char *separators[] = { ":", "#", ".", "_" } ;
991 uInt nsep = 4 ;
992 for ( uInt i = 0 ; i < nsep ; i++ ) {
993 if ( tmp.find( separators[i] ) != String::npos )
994 return separators[i] ;
995 }
996 return "" ;
997 }
998 Int getSrcType( Int stateId )
999 {
1000 // get values
1001 Bool sig ;
1002 getScalar( "SIG", stateId, statetab, sig ) ;
1003 Bool ref ;
1004 getScalar( "REF", stateId, statetab, ref ) ;
1005 Double cal ;
1006 getScalar( "CAL", stateId, statetab, cal ) ;
1007 String obsmode ;
1008 getScalar( "OBS_MODE", stateId, statetab, obsmode ) ;
1009 String sep = detectSeparator( obsmode ) ;
1010
1011 Int srcType = SrcType::NOTYPE ;
1012 if ( sep == ":" )
1013 srcTypeGBT( srcType, sep, obsmode, sig, ref, cal ) ;
1014 else if ( sep == "." || sep == "#" )
1015 srcTypeALMA( srcType, sep, obsmode ) ;
1016 else if ( sep == "_" )
1017 srcTypeOldALMA( srcType, sep, obsmode, sig, ref ) ;
1018 else
1019 srcTypeDefault( srcType, sig, ref ) ;
1020
1021 return srcType ;
1022 }
1023 void srcTypeDefault( Int &st, Bool &sig, Bool &ref )
1024 {
1025 if ( sig ) st = SrcType::SIG ;
1026 else if ( ref ) st = SrcType::REF ;
1027 }
1028 void srcTypeGBT( Int &st, String &sep, String &mode, Bool &sig, Bool &ref, Double &cal )
1029 {
1030 Int epos = mode.find_first_of( sep ) ;
1031 Int nextpos = mode.find_first_of( sep, epos+1 ) ;
1032 String m1 = mode.substr( 0, epos ) ;
1033 String m2 = mode.substr( epos+1, nextpos-epos-1 ) ;
1034 if ( m1 == "Nod" ) {
1035 st = SrcType::NOD ;
1036 }
1037 else if ( m1 == "OffOn" ) {
1038 if ( m2 == "PSWITCHON" ) st = SrcType::PSON ;
1039 if ( m2 == "PSWITCHOFF" ) st = SrcType::PSOFF ;
1040 }
1041 else {
1042 if ( m2 == "FSWITCH" ) {
1043 if ( sig ) st = SrcType::FSON ;
1044 else if ( ref ) st = SrcType::FSOFF ;
1045 }
1046 }
1047 if ( cal > 0.0 ) {
1048 if ( st == SrcType::NOD )
1049 st = SrcType::NODCAL ;
1050 else if ( st == SrcType::PSON )
1051 st = SrcType::PONCAL ;
1052 else if ( st == SrcType::PSOFF )
1053 st = SrcType::POFFCAL ;
1054 else if ( st == SrcType::FSON )
1055 st = SrcType::FONCAL ;
1056 else if ( st == SrcType::FSOFF )
1057 st = SrcType::FOFFCAL ;
1058 else
1059 st = SrcType::CAL ;
1060 }
1061 }
1062 void srcTypeALMA( Int &st, String &sep, String &mode )
1063 {
1064 Int epos = mode.find_first_of( "," ) ;
1065 String first = mode.substr( 0, epos ) ;
1066 epos = first.find_first_of( sep ) ;
1067 Int nextpos = first.find_first_of( sep, epos+1 ) ;
1068 String m1 = first.substr( 0, epos ) ;
1069 String m2 = first.substr( epos+1, nextpos-epos-1 ) ;
1070 if ( m1.find( "CALIBRATE_" ) == 0 ) {
1071 if ( m2.find( "ON_SOURCE" ) == 0 )
1072 st = SrcType::PONCAL ;
1073 else if ( m2.find( "OFF_SOURCE" ) == 0 )
1074 st = SrcType::POFFCAL ;
1075 }
1076 else if ( m1.find( "OBSERVE_TARGET" ) == 0 ) {
1077 if ( m2.find( "ON_SOURCE" ) == 0 )
1078 st = SrcType::PSON ;
1079 else if ( m2.find( "OFF_SOURCE" ) == 0 )
1080 st = SrcType::PSOFF ;
1081 }
1082 }
1083 void srcTypeOldALMA( Int &st, String &sep, String &mode, Bool &sig, Bool &ref )
1084 {
1085 Int epos = mode.find_first_of( "," ) ;
1086 String first = mode.substr( 0, epos ) ;
1087 string substr[4] ;
1088 int numSubstr = split( first, substr, 4, sep ) ;
1089 String m1( substr[0] ) ;
1090 String m2( substr[2] ) ;
1091 if ( numSubstr == 4 ) {
1092 if ( m1.find( "CALIBRATE" ) == 0 ) {
1093 if ( m2.find( "ON" ) == 0 )
1094 st = SrcType::PONCAL ;
1095 else if ( m2.find( "OFF" ) == 0 )
1096 st = SrcType::POFFCAL ;
1097 }
1098 else if ( m1.find( "OBSERVE" ) == 0 ) {
1099 if ( m2.find( "ON" ) == 0 )
1100 st = SrcType::PSON ;
1101 else if ( m2.find( "OFF" ) == 0 )
1102 st = SrcType::PSOFF ;
1103 }
1104 }
1105 else {
1106 if ( sig ) st = SrcType::SIG ;
1107 else if ( ref ) st = SrcType::REF ;
1108 }
1109 }
1110 Block<uInt> getPolNos( Vector<Int> &corr )
1111 {
1112 Block<uInt> polnos( npol ) ;
1113 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
1114 if ( corr[ipol] == Stokes::I || corr[ipol] == Stokes::RR || corr[ipol] == Stokes::XX )
1115 polnos[ipol] = 0 ;
1116 else if ( corr[ipol] == Stokes::Q || corr[ipol] == Stokes::LL || corr[ipol] == Stokes::YY )
1117 polnos[ipol] = 1 ;
1118 else if ( corr[ipol] == Stokes::U || corr[ipol] == Stokes::RL || corr[ipol] == Stokes::XY )
1119 polnos[ipol] = 2 ;
1120 else if ( corr[ipol] == Stokes::V || corr[ipol] == Stokes::LR || corr[ipol] == Stokes::YX )
1121 polnos[ipol] = 3 ;
1122 }
1123 return polnos ;
1124 }
1125 String getPolType( Int &corr )
1126 {
1127 String poltype = "" ;
1128 if ( corr == Stokes::I || corr == Stokes::Q || corr == Stokes::U || corr == Stokes::V )
1129 poltype = "stokes" ;
1130 else if ( corr == Stokes::XX || corr == Stokes::YY || corr == Stokes::XY || corr == Stokes::YX )
1131 poltype = "linear" ;
1132 else if ( corr == Stokes::RR || corr == Stokes::LL || corr == Stokes::RL || corr == Stokes::LR )
1133 poltype = "circular" ;
1134 else if ( corr == Stokes::Plinear || corr == Stokes::Pangle )
1135 poltype = "linpol" ;
1136 return poltype ;
1137 }
1138 uInt getWeatherId()
1139 {
1140 // if only one row, return 0
1141 if ( weatherTime.nelements() == 1 )
1142 return 0 ;
1143
1144 // @todo At the moment, do binary search every time
1145 // if this is bottleneck, frequency of binary search must be reduced
1146 Double t = currentTime.get( "s" ).getValue() ;
1147 uInt idx = min( binarySearch( weatherTime, t ), weatherTime.nelements()-1 ) ;
1148 if ( weatherTime[idx] < t ) {
1149 if ( idx != weatherTime.nelements()-1 ) {
1150 if ( weatherTime[idx+1] - t < 0.5 * weatherInterval[idx+1] )
1151 idx++ ;
1152 }
1153 }
1154 else if ( weatherTime[idx] > t ) {
1155 if ( idx != 0 ) {
1156 if ( weatherTime[idx] - t > 0.5 * weatherInterval[idx] )
1157 idx-- ;
1158 }
1159 }
1160 return idx ;
1161 }
1162 void processSysCal( Int &spwId )
1163 {
1164 // get feedId from row
1165 Int feedId = (Int)tablerow.record().asuInt( "BEAMNO" ) ;
1166
1167 uInt nrow = sctab.nrow() ;
1168 ROScalarColumn<Int> col( sctab, "ANTENNA_ID" ) ;
1169 Vector<Int> aids = col.getColumn() ;
1170 col.attach( sctab, "FEED_ID" ) ;
1171 Vector<Int> fids = col.getColumn() ;
1172 col.attach( sctab, "SPECTRAL_WINDOW_ID" ) ;
1173 Vector<Int> sids = col.getColumn() ;
1174 ROScalarColumn<Double> timeCol( sctab, "TIME" ) ;
1175 ROScalarColumn<Double> intCol( sctab, "INTERVAL" ) ;
1176 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
1177 if ( aids[irow] == antennaId
1178 && fids[irow] == feedId
1179 && sids[irow] == spwId ) {
1180 syscalRow[numSysCalRow] = irow ;
1181 syscalTime[numSysCalRow] = timeCol( irow ) ;
1182 syscalInterval[numSysCalRow] = intCol( irow ) ;
1183 numSysCalRow++ ;
1184 }
1185 }
1186 }
1187 uInt getSysCalIndex()
1188 {
1189 // if only one row, return 0
1190 if ( numSysCalRow == 1 || !isSysCal )
1191 return 0 ;
1192
1193 // @todo At the moment, do binary search every time
1194 // if this is bottleneck, frequency of binary search must be reduced
1195 Double t = currentTime.get( "s" ).getValue() ;
1196 Vector<Double> tslice = syscalTime( Slice(0, numSysCalRow) ) ;
1197 uInt idx = min( binarySearch( tslice, t ), numSysCalRow-1 ) ;
1198 if ( syscalTime[idx] < t ) {
1199 if ( idx != numSysCalRow-1 ) {
1200 if ( syscalTime[idx+1] - t < 0.5 * syscalInterval[idx+1] )
1201 idx++ ;
1202 }
1203 }
1204 else if ( syscalTime[idx] > t ) {
1205 if ( idx != 0 ) {
1206 if ( syscalTime[idx] - t > 0.5 * syscalInterval[idx] )
1207 idx-- ;
1208 }
1209 }
1210 return idx ;
1211 }
1212 Block<uInt> getTcalId( Double &t )
1213 {
1214 // return 0 if no SysCal table
1215 if ( !isSysCal or !isTcal ) {
1216 return Block<uInt>( 4, 0 ) ;
1217 }
1218
1219 // get feedId from row
1220 Int feedId = (Int)tablerow.record().asuInt( "BEAMNO" ) ;
1221
1222 // key
1223 String key = keyTcal( feedId, spwId, t ) ;
1224
1225 // retrieve ids
1226 Vector<uInt> ids = syscalRecord.asArrayuInt( key ) ;
1227 //Vector<uInt> ids = syscalRecord[key] ;
1228 uInt np = ids[1] - ids[0] + 1 ;
1229 Block<uInt> tcalids( np ) ;
1230 if ( np > 0 ) {
1231 tcalids[0] = ids[0] ;
1232 if ( np > 1 ) {
1233 tcalids[1] = ids[1] ;
1234 for ( uInt ip = 2 ; ip < np ; ip++ )
1235 tcalids[ip] = ids[0] + ip - 1 ;
1236 }
1237 }
1238 return tcalids ;
1239 }
1240 uInt maxNumPol()
1241 {
1242 ROScalarColumn<Int> numCorrCol( poltab, "NUM_CORR" ) ;
1243 return max( numCorrCol.getColumn() ) ;
1244 }
1245
1246 Scantable &scantable;
1247 Int antennaId;
1248 uInt rowidx;
1249 String dataColumnName;
1250 TableRow tablerow;
1251 STHeader header;
1252 Vector<Int> feedEntry;
1253 uInt nbeam;
1254 Int npol;
1255 Int nchan;
1256 Int sourceId;
1257 Int polId;
1258 Int spwId;
1259 uInt cycleNo;
1260 MDirection sourceDir;
1261 MPosition antpos;
1262 MEpoch currentTime;
1263 map<Int,uInt> ifmap;
1264 Block<uInt> polnos;
1265 Bool getpt;
1266 Vector<Double> pointingTime;
1267 Cube<Double> pointingDirection;
1268 MDirection::Types dirType;
1269 Bool isWeather;
1270 Vector<Double> weatherTime;
1271 Vector<Double> weatherInterval;
1272 Bool isSysCal;
1273 Bool isTcal;
1274 Record syscalRecord;
1275 //map< String,Vector<uInt> > syscalRecord;
1276 uInt numSysCalRow ;
1277 Vector<uInt> syscalRow;
1278 Vector<Double> syscalTime;
1279 Vector<Double> syscalInterval;
1280 //String tsysCol;
1281 //String tcalCol;
1282
1283 // MS subtables
1284 Table obstab;
1285 Table sctab;
1286 Table spwtab;
1287 Table statetab;
1288 Table ddtab;
1289 Table poltab;
1290 Table fieldtab;
1291 Table anttab;
1292 Table srctab;
1293 Matrix<Float> sp;
1294 Matrix<uChar> fl;
1295 MeasFrame mf;
1296
1297 // MS MAIN columns
1298 ROTableColumn intervalCol;
1299 ROTableColumn flagRowCol;
1300 ROArrayColumn<Float> floatDataCol;
1301 ROArrayColumn<Complex> dataCol;
1302 ROArrayColumn<Bool> flagCol;
1303
1304 // MS SYSCAL columns
1305 ROArrayColumn<Float> sysCalTsysCol;
1306
1307 // Scantable MAIN columns
1308 RecordFieldPtr<Double> timeRF,intervalRF,sourceVelocityRF;
1309 RecordFieldPtr< Vector<Double> > directionRF,scanRateRF,
1310 sourceProperMotionRF,sourceDirectionRF;
1311 RecordFieldPtr<Float> azimuthRF,elevationRF;
1312 RecordFieldPtr<uInt> weatherIdRF,cycleNoRF,flagRowRF,polNoRF,tcalIdRF,
1313 ifNoRF,freqIdRF,moleculeIdRF,beamNoRF,focusIdRF,scanNoRF;
1314 RecordFieldPtr< Vector<Float> > spectraRF,tsysRF;
1315 RecordFieldPtr< Vector<uChar> > flagtraRF;
1316 RecordFieldPtr<String> sourceNameRF,fieldNameRF;
1317 RecordFieldPtr<Int> sourceTypeRF;
1318};
1319
1320class BaseTcalVisitor: public TableVisitor {
1321 uInt lastRecordNo ;
1322 Int lastAntennaId ;
1323 Int lastFeedId ;
1324 Int lastSpwId ;
1325 Double lastTime ;
1326protected:
1327 const Table &table;
1328 uInt count;
1329public:
1330 BaseTcalVisitor(const Table &table)
1331 : table(table)
1332 {
1333 count = 0;
1334 }
1335
1336 virtual void enterAntennaId(const uInt recordNo, Int columnValue) { }
1337 virtual void leaveAntennaId(const uInt recordNo, Int columnValue) { }
1338 virtual void enterFeedId(const uInt recordNo, Int columnValue) { }
1339 virtual void leaveFeedId(const uInt recordNo, Int columnValue) { }
1340 virtual void enterSpwId(const uInt recordNo, Int columnValue) { }
1341 virtual void leaveSpwId(const uInt recordNo, Int columnValue) { }
1342 virtual void enterTime(const uInt recordNo, Double columnValue) { }
1343 virtual void leaveTime(const uInt recordNo, Double columnValue) { }
1344
1345 virtual Bool visitRecord(const uInt recordNo,
1346 const Int antennaId,
1347 const Int feedId,
1348 const Int spwId,
1349 const Double time) { return True ; }
1350
1351 virtual Bool visit(Bool isFirst, const uInt recordNo,
1352 const uInt nCols, void const *const colValues[]) {
1353 Int antennaId, feedId, spwId;
1354 Double time;
1355 { // prologue
1356 uInt i = 0;
1357 {
1358 const Int *col = (const Int *)colValues[i++];
1359 antennaId = col[recordNo];
1360 }
1361 {
1362 const Int *col = (const Int *)colValues[i++];
1363 feedId = col[recordNo];
1364 }
1365 {
1366 const Int *col = (const Int *)colValues[i++];
1367 spwId = col[recordNo];
1368 }
1369 {
1370 const Double *col = (const Double *)colValues[i++];
1371 time = col[recordNo];
1372 }
1373 assert(nCols == i);
1374 }
1375
1376 if (isFirst) {
1377 enterAntennaId(recordNo, antennaId);
1378 enterFeedId(recordNo, feedId);
1379 enterSpwId(recordNo, spwId);
1380 enterTime(recordNo, time);
1381 } else {
1382 if ( lastAntennaId != antennaId ) {
1383 leaveTime(lastRecordNo, lastTime);
1384 leaveSpwId(lastRecordNo, lastSpwId);
1385 leaveFeedId(lastRecordNo, lastFeedId);
1386 leaveAntennaId(lastRecordNo, lastAntennaId);
1387
1388 enterAntennaId(recordNo, antennaId);
1389 enterFeedId(recordNo, feedId);
1390 enterSpwId(recordNo, spwId);
1391 enterTime(recordNo, time);
1392 }
1393 else if (lastFeedId != feedId) {
1394 leaveTime(lastRecordNo, lastTime);
1395 leaveSpwId(lastRecordNo, lastSpwId);
1396 leaveFeedId(lastRecordNo, lastFeedId);
1397
1398 enterFeedId(recordNo, feedId);
1399 enterSpwId(recordNo, spwId);
1400 enterTime(recordNo, time);
1401 } else if (lastSpwId != spwId) {
1402 leaveTime(lastRecordNo, lastTime);
1403 leaveSpwId(lastRecordNo, lastSpwId);
1404
1405 enterSpwId(recordNo, spwId);
1406 enterTime(recordNo, time);
1407 } else if (lastTime != time) {
1408 leaveTime(lastRecordNo, lastTime);
1409 enterTime(recordNo, time);
1410 }
1411 }
1412 count++;
1413 Bool result = visitRecord(recordNo, antennaId, feedId, spwId, time);
1414
1415 { // epilogue
1416 lastRecordNo = recordNo;
1417
1418 lastAntennaId = antennaId;
1419 lastFeedId = feedId;
1420 lastSpwId = spwId;
1421 lastTime = time;
1422 }
1423 return result ;
1424 }
1425
1426 virtual void finish() {
1427 if (count > 0) {
1428 leaveTime(lastRecordNo, lastTime);
1429 leaveSpwId(lastRecordNo, lastSpwId);
1430 leaveFeedId(lastRecordNo, lastFeedId);
1431 leaveAntennaId(lastRecordNo, lastAntennaId);
1432 }
1433 }
1434};
1435
1436class TcalVisitor: public BaseTcalVisitor, public MSFillerUtils {
1437public:
1438 TcalVisitor(const Table &table, Table &tcaltab, Record &r, Int aid )
1439 //TcalVisitor(const Table &table, Table &tcaltab, map< String,Vector<uInt> > &r, Int aid )
1440 : BaseTcalVisitor( table ),
1441 tcal(tcaltab),
1442 rec(r),
1443 antenna(aid)
1444 {
1445 process = False ;
1446 rowidx = 0 ;
1447
1448 // attach to SYSCAL columns
1449 timeCol.attach( table, "TIME" ) ;
1450
1451 // add rows
1452 uInt addrow = table.nrow() * 4 ;
1453 tcal.addRow( addrow ) ;
1454
1455 // attach to TCAL columns
1456 row = TableRow( tcal ) ;
1457 TableRecord &trec = row.record() ;
1458 idRF.attachToRecord( trec, "ID" ) ;
1459 timeRF.attachToRecord( trec, "TIME" ) ;
1460 tcalRF.attachToRecord( trec, "TCAL" ) ;
1461 }
1462
1463 virtual void enterAntennaId(const uInt recordNo, Int columnValue) {
1464 if ( columnValue == antenna )
1465 process = True ;
1466 }
1467 virtual void leaveAntennaId(const uInt recordNo, Int columnValue) {
1468 process = False ;
1469 }
1470 virtual void enterFeedId(const uInt recordNo, Int columnValue) { }
1471 virtual void leaveFeedId(const uInt recordNo, Int columnValue) { }
1472 virtual void enterSpwId(const uInt recordNo, Int columnValue) { }
1473 virtual void leaveSpwId(const uInt recordNo, Int columnValue) { }
1474 virtual void enterTime(const uInt recordNo, Double columnValue) {
1475 qtime = timeCol( recordNo ) ;
1476 }
1477 virtual void leaveTime(const uInt recordNo, Double columnValue) { }
1478 virtual Bool visitRecord(const uInt recordNo,
1479 const Int antennaId,
1480 const Int feedId,
1481 const Int spwId,
1482 const Double time)
1483 {
1484 //cout << "(" << recordNo << "," << antennaId << "," << feedId << "," << spwId << ")" << endl ;
1485 if ( process ) {
1486 String sTime = MVTime( qtime ).string( MVTime::YMD ) ;
1487 *timeRF = sTime ;
1488 uInt oldidx = rowidx ;
1489 Matrix<Float> subtcal = tcalCol( recordNo ) ;
1490 Vector<uInt> idminmax( 2 ) ;
1491 for ( uInt ipol = 0 ; ipol < subtcal.nrow() ; ipol++ ) {
1492 *idRF = rowidx ;
1493 tcalRF.define( subtcal.row( ipol ) ) ;
1494
1495 // commit row
1496 row.put( rowidx ) ;
1497 rowidx++ ;
1498 }
1499
1500 idminmax[0] = oldidx ;
1501 idminmax[1] = rowidx - 1 ;
1502
1503 String key = keyTcal( feedId, spwId, sTime ) ;
1504 rec.define( key, idminmax ) ;
1505 //rec[key] = idminmax ;
1506 }
1507 return True ;
1508 }
1509 virtual void finish()
1510 {
1511 BaseTcalVisitor::finish() ;
1512
1513 if ( tcal.nrow() > rowidx ) {
1514 uInt numRemove = tcal.nrow() - rowidx ;
1515 //cout << "numRemove = " << numRemove << endl ;
1516 Vector<uInt> rows( numRemove ) ;
1517 indgen( rows, rowidx ) ;
1518 tcal.removeRow( rows ) ;
1519 }
1520
1521 }
1522 void setTcalColumn( String &col )
1523 {
1524 //colName = col ;
1525 tcalCol.attach( table, col ) ;
1526 }
1527private:
1528 Table &tcal;
1529 Record &rec;
1530 //map< String,Vector<uInt> > &rec;
1531 Int antenna;
1532 uInt rowidx;
1533 Bool process;
1534 Quantum<Double> qtime;
1535 TableRow row;
1536 String colName;
1537
1538 // MS SYSCAL columns
1539 ROScalarQuantColumn<Double> timeCol;
1540 ROArrayColumn<Float> tcalCol;
1541
1542 // TCAL columns
1543 RecordFieldPtr<uInt> idRF;
1544 RecordFieldPtr<String> timeRF;
1545 RecordFieldPtr< Vector<Float> > tcalRF;
1546};
1547
1548MSFiller::MSFiller( casa::CountedPtr<Scantable> stable )
1549 : table_( stable ),
1550 tablename_( "" ),
1551 antenna_( -1 ),
1552 antennaStr_(""),
1553 getPt_( True ),
1554 isFloatData_( False ),
1555 isData_( False ),
1556 isDoppler_( False ),
1557 isFlagCmd_( False ),
1558 isFreqOffset_( False ),
1559 isHistory_( False ),
1560 isProcessor_( False ),
1561 isSysCal_( False ),
1562 isWeather_( False ),
1563 colTsys_( "TSYS_SPECTRUM" ),
1564 colTcal_( "TCAL_SPECTRUM" )
1565{
1566 os_ = LogIO() ;
1567 os_.origin( LogOrigin( "MSFiller", "MSFiller()", WHERE ) ) ;
1568}
1569
1570MSFiller::~MSFiller()
1571{
1572 os_.origin( LogOrigin( "MSFiller", "~MSFiller()", WHERE ) ) ;
1573}
1574
1575bool MSFiller::open( const std::string &filename, const casa::Record &rec )
1576{
1577 os_.origin( LogOrigin( "MSFiller", "open()", WHERE ) ) ;
1578 //double startSec = mathutil::gettimeofday_sec() ;
1579 //os_ << "start MSFiller::open() startsec=" << startSec << LogIO::POST ;
1580 //os_ << " filename = " << filename << endl ;
1581
1582 // parsing MS options
1583 if ( rec.isDefined( "ms" ) ) {
1584 Record msrec = rec.asRecord( "ms" ) ;
1585 if ( msrec.isDefined( "getpt" ) ) {
1586 getPt_ = msrec.asBool( "getpt" ) ;
1587 }
1588 if ( msrec.isDefined( "antenna" ) ) {
1589 if ( msrec.type( msrec.fieldNumber( "antenna" ) ) == TpInt ) {
1590 antenna_ = msrec.asInt( "antenna" ) ;
1591 }
1592 else {
1593 //antenna_ = atoi( msrec.asString( "antenna" ).c_str() ) ;
1594 antennaStr_ = msrec.asString( "antenna" ) ;
1595 }
1596 }
1597 else {
1598 antenna_ = 0 ;
1599 }
1600 }
1601
1602 MeasurementSet *tmpMS = new MeasurementSet( filename, Table::Old ) ;
1603 tablename_ = tmpMS->tableName() ;
1604 if ( antenna_ == -1 && antennaStr_.size() > 0 ) {
1605 MSAntennaIndex msAntIdx( tmpMS->antenna() ) ;
1606 Vector<Int> id = msAntIdx.matchAntennaName( antennaStr_ ) ;
1607 if ( id.size() > 0 )
1608 antenna_ = id[0] ;
1609 else {
1610 delete tmpMS ;
1611 //throw( AipsError( "Antenna " + antennaStr_ + " doesn't exist." ) ) ;
1612 os_ << LogIO::SEVERE << "Antenna " << antennaStr_ << " doesn't exist." << LogIO::POST ;
1613 return False ;
1614 }
1615 }
1616
1617 os_ << "Parsing MS options" << endl ;
1618 os_ << " getPt = " << getPt_ << endl ;
1619 os_ << " antenna = " << antenna_ << endl ;
1620 os_ << " antennaStr = " << antennaStr_ << LogIO::POST ;
1621
1622 mstable_ = MeasurementSet( (*tmpMS)( tmpMS->col("ANTENNA1") == antenna_
1623 && tmpMS->col("ANTENNA1") == tmpMS->col("ANTENNA2") ) ) ;
1624
1625 delete tmpMS ;
1626
1627 // check which data column exists
1628 isFloatData_ = mstable_.tableDesc().isColumn( "FLOAT_DATA" ) ;
1629 isData_ = mstable_.tableDesc().isColumn( "DATA" ) ;
1630
1631 //double endSec = mathutil::gettimeofday_sec() ;
1632 //os_ << "end MSFiller::open() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1633 return true ;
1634}
1635
1636void MSFiller::fill()
1637{
1638 //double startSec = mathutil::gettimeofday_sec() ;
1639 //os_ << "start MSFiller::fill() startSec=" << startSec << LogIO::POST ;
1640
1641 os_.origin( LogOrigin( "MSFiller", "fill()", WHERE ) ) ;
1642
1643 // Initialize header
1644 STHeader sdh ;
1645 initHeader( sdh ) ;
1646 table_->setHeader( sdh ) ;
1647
1648 // check if optional table exists
1649 const TableRecord &msrec = mstable_.keywordSet() ;
1650 isDoppler_ = msrec.isDefined( "DOPPLER" ) ;
1651 if ( isDoppler_ )
1652 if ( mstable_.doppler().nrow() == 0 )
1653 isDoppler_ = False ;
1654 isFlagCmd_ = msrec.isDefined( "FLAG_CMD" ) ;
1655 if ( isFlagCmd_ )
1656 if ( mstable_.flagCmd().nrow() == 0 )
1657 isFlagCmd_ = False ;
1658 isFreqOffset_ = msrec.isDefined( "FREQ_OFFSET" ) ;
1659 if ( isFreqOffset_ )
1660 if ( mstable_.freqOffset().nrow() == 0 )
1661 isFreqOffset_ = False ;
1662 isHistory_ = msrec.isDefined( "HISTORY" ) ;
1663 if ( isHistory_ )
1664 if ( mstable_.history().nrow() == 0 )
1665 isHistory_ = False ;
1666 isProcessor_ = msrec.isDefined( "PROCESSOR" ) ;
1667 if ( isProcessor_ )
1668 if ( mstable_.processor().nrow() == 0 )
1669 isProcessor_ = False ;
1670 isSysCal_ = msrec.isDefined( "SYSCAL" ) ;
1671 if ( isSysCal_ )
1672 if ( mstable_.sysCal().nrow() == 0 )
1673 isSysCal_ = False ;
1674 isWeather_ = msrec.isDefined( "WEATHER" ) ;
1675 if ( isWeather_ )
1676 if ( mstable_.weather().nrow() == 0 )
1677 isWeather_ = False ;
1678
1679 // column name for Tsys and Tcal
1680 if ( isSysCal_ ) {
1681 const MSSysCal &caltab = mstable_.sysCal() ;
1682 if ( !caltab.tableDesc().isColumn( colTcal_ ) ) {
1683 colTcal_ = "TCAL" ;
1684 if ( !caltab.tableDesc().isColumn( colTcal_ ) )
1685 colTcal_ = "NONE" ;
1686 }
1687 if ( !caltab.tableDesc().isColumn( colTsys_ ) ) {
1688 colTsys_ = "TSYS" ;
1689 if ( !caltab.tableDesc().isColumn( colTcal_ ) )
1690 colTsys_ = "NONE" ;
1691 }
1692 }
1693 else {
1694 colTcal_ = "NONE" ;
1695 colTsys_ = "NONE" ;
1696 }
1697
1698 // Access to MS subtables
1699 //MSField &fieldtab = mstable_.field() ;
1700 //MSPolarization &poltab = mstable_.polarization() ;
1701 //MSDataDescription &ddtab = mstable_.dataDescription() ;
1702 //MSObservation &obstab = mstable_.observation() ;
1703 //MSSource &srctab = mstable_.source() ;
1704 //MSSpectralWindow &spwtab = mstable_.spectralWindow() ;
1705 //MSSysCal &caltab = mstable_.sysCal() ;
1706 MSPointing &pointtab = mstable_.pointing() ;
1707 //MSState &stattab = mstable_.state() ;
1708 //MSAntenna &anttab = mstable_.antenna() ;
1709
1710 // SUBTABLES: FREQUENCIES
1711 //string freqFrame = getFrame() ;
1712 string freqFrame = "LSRK" ;
1713 table_->frequencies().setFrame( freqFrame ) ;
1714 table_->frequencies().setFrame( freqFrame, True ) ;
1715
1716 // SUBTABLES: WEATHER
1717 fillWeather() ;
1718
1719 // SUBTABLES: FOCUS
1720 fillFocus() ;
1721
1722 // SUBTABLES: TCAL
1723 fillTcal() ;
1724
1725 // SUBTABLES: FIT
1726 //fillFit() ;
1727
1728 // SUBTABLES: HISTORY
1729 //fillHistory() ;
1730
1731 /***
1732 * Start iteration using TableVisitor
1733 ***/
1734 Table stab = table_->table() ;
1735 {
1736 static const char *cols[] = {
1737 "OBSERVATION_ID", "FEED1", "FIELD_ID", "DATA_DESC_ID", "SCAN_NUMBER",
1738 "STATE_ID", "TIME",
1739 NULL
1740 };
1741 static const TypeManagerImpl<Int> tmInt;
1742 static const TypeManagerImpl<Double> tmDouble;
1743 static const TypeManager *const tms[] = {
1744 &tmInt, &tmInt, &tmInt, &tmInt, &tmInt, &tmInt, &tmDouble, NULL
1745 };
1746 //double t0 = mathutil::gettimeofday_sec() ;
1747 MSFillerVisitor myVisitor(mstable_, *table_ );
1748 //double t1 = mathutil::gettimeofday_sec() ;
1749 //cout << "MSFillerVisitor(): elapsed time " << t1-t0 << " sec" << endl ;
1750 myVisitor.setAntenna( antenna_ ) ;
1751 //myVisitor.setHeader( sdh ) ;
1752 if ( getPt_ ) {
1753 Table ptsel = pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort( "TIME" ) ;
1754 myVisitor.setPointingTable( ptsel ) ;
1755 }
1756 if ( isWeather_ )
1757 myVisitor.setWeatherTime( mwTime_, mwInterval_ ) ;
1758 if ( isSysCal_ )
1759 myVisitor.setSysCalRecord( tcalrec_ ) ;
1760
1761 //double t2 = mathutil::gettimeofday_sec() ;
1762 traverseTable(mstable_, cols, tms, &myVisitor);
1763 //double t3 = mathutil::gettimeofday_sec() ;
1764 //cout << "traverseTable(): elapsed time " << t3-t2 << " sec" << endl ;
1765
1766 sdh = myVisitor.getHeader() ;
1767 }
1768 /***
1769 * End iteration using TableVisitor
1770 ***/
1771
1772 // set header
1773 //sdh = myVisitor.getHeader() ;
1774 //table_->setHeader( sdh ) ;
1775
1776 // save path to POINTING table
1777 // 2011/07/06 TN
1778 // Path to POINTING table in original MS will not be written
1779 // if getPt_ is True
1780 Path datapath( tablename_ ) ;
1781 if ( !getPt_ ) {
1782 String pTabName = datapath.absoluteName() + "/POINTING" ;
1783 stab.rwKeywordSet().define( "POINTING", pTabName ) ;
1784 }
1785
1786 // for GBT
1787 if ( sdh.antennaname.contains( "GBT" ) ) {
1788 String goTabName = datapath.absoluteName() + "/GBT_GO" ;
1789 stab.rwKeywordSet().define( "GBT_GO", goTabName ) ;
1790 }
1791
1792 // for MS created from ASDM
1793 const TableRecord &msKeys = mstable_.keywordSet() ;
1794 uInt nfields = msKeys.nfields() ;
1795 for ( uInt ifield = 0 ; ifield < nfields ; ifield++ ) {
1796 String name = msKeys.name( ifield ) ;
1797 //os_ << "name = " << name << LogIO::POST ;
1798 if ( name.find( "ASDM" ) != String::npos ) {
1799 String asdmpath = msKeys.asTable( ifield ).tableName() ;
1800 os_ << "ASDM table: " << asdmpath << LogIO::POST ;
1801 stab.rwKeywordSet().define( name, asdmpath ) ;
1802 }
1803 }
1804
1805 //double endSec = mathutil::gettimeofday_sec() ;
1806 //os_ << "end MSFiller::fill() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1807}
1808
1809void MSFiller::close()
1810{
1811 //tablesel_.closeSubTables() ;
1812 mstable_.closeSubTables() ;
1813 //tablesel_.unlock() ;
1814 mstable_.unlock() ;
1815}
1816
1817void MSFiller::fillWeather()
1818{
1819 //double startSec = mathutil::gettimeofday_sec() ;
1820 //os_ << "start MSFiller::fillWeather() startSec=" << startSec << LogIO::POST ;
1821
1822 if ( !isWeather_ ) {
1823 // add dummy row
1824 table_->weather().table().addRow(1,True) ;
1825 return ;
1826 }
1827
1828 Table mWeather = mstable_.weather() ;
1829 //Table mWeatherSel = mWeather( mWeather.col("ANTENNA_ID") == antenna_ ).sort("TIME") ;
1830 Table mWeatherSel( mWeather( mWeather.col("ANTENNA_ID") == antenna_ ).sort("TIME") ) ;
1831 //os_ << "mWeatherSel.nrow() = " << mWeatherSel.nrow() << LogIO::POST ;
1832 if ( mWeatherSel.nrow() == 0 ) {
1833 os_ << "No rows with ANTENNA_ID = " << antenna_ << " in WEATHER table, Try -1..." << LogIO::POST ;
1834 mWeatherSel = Table( MSWeather( mWeather( mWeather.col("ANTENNA_ID") == -1 ) ) ) ;
1835 if ( mWeatherSel.nrow() == 0 ) {
1836 os_ << "No rows in WEATHER table" << LogIO::POST ;
1837 }
1838 }
1839 uInt wnrow = mWeatherSel.nrow() ;
1840 //os_ << "wnrow = " << wnrow << LogIO::POST ;
1841
1842 if ( wnrow == 0 )
1843 return ;
1844
1845 Table wtab = table_->weather().table() ;
1846 wtab.addRow( wnrow ) ;
1847
1848 Bool stationInfoExists = mWeatherSel.tableDesc().isColumn( "NS_WX_STATION_ID" ) ;
1849 Int stationId = -1 ;
1850 if ( stationInfoExists ) {
1851 // determine which station is closer
1852 ROScalarColumn<Int> stationCol( mWeatherSel, "NS_WX_STATION_ID" ) ;
1853 ROArrayColumn<Double> stationPosCol( mWeatherSel, "NS_WX_STATION_POSITION" ) ;
1854 Vector<Int> stationIds = stationCol.getColumn() ;
1855 Vector<Int> stationIdList( 0 ) ;
1856 Matrix<Double> stationPosList( 0, 3, 0.0 ) ;
1857 uInt numStation = 0 ;
1858 for ( uInt i = 0 ; i < stationIds.size() ; i++ ) {
1859 if ( !anyEQ( stationIdList, stationIds[i] ) ) {
1860 numStation++ ;
1861 stationIdList.resize( numStation, True ) ;
1862 stationIdList[numStation-1] = stationIds[i] ;
1863 stationPosList.resize( numStation, 3, True ) ;
1864 stationPosList.row( numStation-1 ) = stationPosCol( i ) ;
1865 }
1866 }
1867 //os_ << "staionIdList = " << stationIdList << endl ;
1868 Table mAntenna = mstable_.antenna() ;
1869 ROArrayColumn<Double> antposCol( mAntenna, "POSITION" ) ;
1870 Vector<Double> antpos = antposCol( antenna_ ) ;
1871 Double minDiff = -1.0 ;
1872 for ( uInt i = 0 ; i < stationIdList.size() ; i++ ) {
1873 Double diff = sum( square( antpos - stationPosList.row( i ) ) ) ;
1874 if ( minDiff < 0.0 || minDiff > diff ) {
1875 minDiff = diff ;
1876 stationId = stationIdList[i] ;
1877 }
1878 }
1879 }
1880 //os_ << "stationId = " << stationId << endl ;
1881
1882 ScalarColumn<Float> *fCol ;
1883 ROScalarColumn<Float> *sharedFloatCol ;
1884 if ( mWeatherSel.tableDesc().isColumn( "TEMPERATURE" ) ) {
1885 fCol = new ScalarColumn<Float>( wtab, "TEMPERATURE" ) ;
1886 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "TEMPERATURE" ) ;
1887 fCol->putColumn( *sharedFloatCol ) ;
1888 delete sharedFloatCol ;
1889 delete fCol ;
1890 }
1891 if ( mWeatherSel.tableDesc().isColumn( "PRESSURE" ) ) {
1892 fCol = new ScalarColumn<Float>( wtab, "PRESSURE" ) ;
1893 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "PRESSURE" ) ;
1894 fCol->putColumn( *sharedFloatCol ) ;
1895 delete sharedFloatCol ;
1896 delete fCol ;
1897 }
1898 if ( mWeatherSel.tableDesc().isColumn( "REL_HUMIDITY" ) ) {
1899 fCol = new ScalarColumn<Float>( wtab, "HUMIDITY" ) ;
1900 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "REL_HUMIDITY" ) ;
1901 fCol->putColumn( *sharedFloatCol ) ;
1902 delete sharedFloatCol ;
1903 delete fCol ;
1904 }
1905 if ( mWeatherSel.tableDesc().isColumn( "WIND_SPEED" ) ) {
1906 fCol = new ScalarColumn<Float>( wtab, "WINDSPEED" ) ;
1907 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "WIND_SPEED" ) ;
1908 fCol->putColumn( *sharedFloatCol ) ;
1909 delete sharedFloatCol ;
1910 delete fCol ;
1911 }
1912 if ( mWeatherSel.tableDesc().isColumn( "WIND_DIRECTION" ) ) {
1913 fCol = new ScalarColumn<Float>( wtab, "WINDAZ" ) ;
1914 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "WIND_DIRECTION" ) ;
1915 fCol->putColumn( *sharedFloatCol ) ;
1916 delete sharedFloatCol ;
1917 delete fCol ;
1918 }
1919 ScalarColumn<uInt> idCol( wtab, "ID" ) ;
1920 for ( uInt irow = 0 ; irow < wnrow ; irow++ )
1921 idCol.put( irow, irow ) ;
1922
1923 ROScalarQuantColumn<Double> tqCol( mWeatherSel, "TIME" ) ;
1924 ROScalarColumn<Double> tCol( mWeatherSel, "TIME" ) ;
1925 String tUnit = tqCol.getUnits() ;
1926 Vector<Double> mwTime = tCol.getColumn() ;
1927 if ( tUnit == "d" )
1928 mwTime *= 86400.0 ;
1929 tqCol.attach( mWeatherSel, "INTERVAL" ) ;
1930 tCol.attach( mWeatherSel, "INTERVAL" ) ;
1931 String iUnit = tqCol.getUnits() ;
1932 Vector<Double> mwInterval = tCol.getColumn() ;
1933 if ( iUnit == "d" )
1934 mwInterval *= 86400.0 ;
1935
1936 if ( stationId > 0 ) {
1937 ROScalarColumn<Int> stationCol( mWeatherSel, "NS_WX_STATION_ID" ) ;
1938 Vector<Int> stationVec = stationCol.getColumn() ;
1939 uInt wsnrow = ntrue( stationVec == stationId ) ;
1940 mwTime_.resize( wsnrow ) ;
1941 mwInterval_.resize( wsnrow ) ;
1942 mwIndex_.resize( wsnrow ) ;
1943 uInt wsidx = 0 ;
1944 for ( uInt irow = 0 ; irow < wnrow ; irow++ ) {
1945 if ( stationId == stationVec[irow] ) {
1946 mwTime_[wsidx] = mwTime[irow] ;
1947 mwInterval_[wsidx] = mwInterval[irow] ;
1948 mwIndex_[wsidx] = irow ;
1949 wsidx++ ;
1950 }
1951 }
1952 }
1953 else {
1954 mwTime_ = mwTime ;
1955 mwInterval_ = mwInterval ;
1956 mwIndex_.resize( mwTime_.size() ) ;
1957 indgen( mwIndex_ ) ;
1958 }
1959 //os_ << "mwTime[0] = " << mwTime_[0] << " mwInterval[0] = " << mwInterval_[0] << LogIO::POST ;
1960 //double endSec = mathutil::gettimeofday_sec() ;
1961 //os_ << "end MSFiller::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1962}
1963
1964void MSFiller::fillFocus()
1965{
1966 //double startSec = mathutil::gettimeofday_sec() ;
1967 //os_ << "start MSFiller::fillFocus() startSec=" << startSec << LogIO::POST ;
1968 // tentative
1969 table_->focus().addEntry( 0.0, 0.0, 0.0, 0.0 ) ;
1970 //double endSec = mathutil::gettimeofday_sec() ;
1971 //os_ << "end MSFiller::fillFocus() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1972}
1973
1974void MSFiller::fillTcal()
1975{
1976 //double startSec = mathutil::gettimeofday_sec() ;
1977 //os_ << "start MSFiller::fillTcal() startSec=" << startSec << LogIO::POST ;
1978
1979 if ( !isSysCal_ ) {
1980 // add dummy row
1981 os_ << "No SYSCAL rows" << LogIO::POST ;
1982 table_->tcal().table().addRow(1,True) ;
1983 Vector<Float> defaultTcal( 1, 1.0 ) ;
1984 ArrayColumn<Float> tcalCol( table_->tcal().table(), "TCAL" ) ;
1985 tcalCol.put( 0, defaultTcal ) ;
1986 return ;
1987 }
1988
1989 if ( colTcal_ == "NONE" ) {
1990 // add dummy row
1991 os_ << "No TCAL column" << LogIO::POST ;
1992 table_->tcal().table().addRow(1,True) ;
1993 Vector<Float> defaultTcal( 1, 1.0 ) ;
1994 ArrayColumn<Float> tcalCol( table_->tcal().table(), "TCAL" ) ;
1995 tcalCol.put( 0, defaultTcal ) ;
1996 return ;
1997 }
1998
1999 Table &sctab = mstable_.sysCal() ;
2000 if ( sctab.nrow() == 0 ) {
2001 os_ << "No SYSCAL rows" << LogIO::POST ;
2002 return ;
2003 }
2004 ROScalarColumn<Int> antCol( sctab, "ANTENNA_ID" ) ;
2005 Vector<Int> ant = antCol.getColumn() ;
2006 if ( allNE( ant, antenna_ ) ) {
2007 os_ << "No SYSCAL rows" << LogIO::POST ;
2008 return ;
2009 }
2010 ROTableColumn tcalCol( sctab, colTcal_ ) ;
2011 Bool notDefined = False ;
2012 for ( uInt irow = 0 ; irow < sctab.nrow() ; irow++ ) {
2013 if ( ant[irow] == antenna_ && !tcalCol.isDefined( irow ) ) {
2014 notDefined = True ;
2015 break ;
2016 }
2017 }
2018 if ( notDefined ) {
2019 os_ << "No TCAL value" << LogIO::POST ;
2020 table_->tcal().table().addRow(1,True) ;
2021 Vector<Float> defaultTcal( 1, 1.0 ) ;
2022 ArrayColumn<Float> tcalCol( table_->tcal().table(), "TCAL" ) ;
2023 tcalCol.put( 0, defaultTcal ) ;
2024 return ;
2025 }
2026
2027 static const char *cols[] = {
2028 "ANTENNA_ID", "FEED_ID", "SPECTRAL_WINDOW_ID", "TIME",
2029 NULL
2030 };
2031 static const TypeManagerImpl<Int> tmInt;
2032 static const TypeManagerImpl<Double> tmDouble;
2033 static const TypeManager *const tms[] = {
2034 &tmInt, &tmInt, &tmInt, &tmDouble, NULL
2035 };
2036 Table tab = table_->tcal().table() ;
2037 TcalVisitor visitor( sctab, tab, tcalrec_, antenna_ ) ;
2038 visitor.setTcalColumn( colTcal_ ) ;
2039
2040 traverseTable(sctab, cols, tms, &visitor);
2041
2042 //tcalrec_.print( std::cout ) ;
2043 //double endSec = mathutil::gettimeofday_sec() ;
2044 //os_ << "end MSFiller::fillTcal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2045}
2046
2047string MSFiller::getFrame()
2048{
2049 MFrequency::Types frame = MFrequency::DEFAULT ;
2050 ROTableColumn numChanCol( mstable_.spectralWindow(), "NUM_CHAN" ) ;
2051 ROTableColumn measFreqRefCol( mstable_.spectralWindow(), "MEAS_FREQ_REF" ) ;
2052 uInt nrow = numChanCol.nrow() ;
2053 Vector<Int> measFreqRef( nrow, MFrequency::DEFAULT ) ;
2054 uInt nref = 0 ;
2055 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
2056 if ( numChanCol.asInt( irow ) != 4 ) { // exclude WVR
2057 measFreqRef[nref] = measFreqRefCol.asInt( irow ) ;
2058 nref++ ;
2059 }
2060 }
2061 if ( nref > 0 )
2062 frame = (MFrequency::Types)measFreqRef[0] ;
2063
2064 return MFrequency::showType( frame ) ;
2065}
2066
2067void MSFiller::initHeader( STHeader &header )
2068{
2069 header.nchan = 0 ;
2070 header.npol = 0 ;
2071 header.nif = 0 ;
2072 header.nbeam = 0 ;
2073 header.observer = "" ;
2074 header.project = "" ;
2075 header.obstype = "" ;
2076 header.antennaname = "" ;
2077 header.antennaposition.resize( 3 ) ;
2078 header.equinox = 0.0 ;
2079 header.freqref = "" ;
2080 header.reffreq = -1.0 ;
2081 header.bandwidth = 0.0 ;
2082 header.utc = 0.0 ;
2083 header.fluxunit = "" ;
2084 header.epoch = "" ;
2085 header.poltype = "" ;
2086}
2087
2088};
Note: See TracBrowser for help on using the repository browser.