Context Navigation

MSFiller.cpp@ 2283

Visit:

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

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

INCREMENT for channel averaged data depends on NET_SIDEBAND value
of MS SPW table.

File size: 68.5 KB

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

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

Context Navigation

source: trunk/src/MSFiller.cpp@ 2283

Download in other formats: