Context Navigation

MSFiller.cpp @ 1993

Visit:

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

New Development: No

JIRA Issue: Yes CAS-2718

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

Small improvement on TCAL(_SPECTRUM) and TSYS(_SPECTRUM) handling.

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

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

Context Navigation

source: trunk/src/MSFiller.cpp @ 1993

Download in other formats: