Context Navigation

source: branches/hpc34/src/MSFiller.cpp@ 2806

Visit:

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

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

Download in other formats: