Context Navigation

source: branches/casa-release-4_3/src/MSFiller.cpp@ 3068

Visit:

Last change on this file since 3068 was 2941, checked in by Malte Marquarding, 11 years ago
Fix const char* warning
File size: 68.6 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: