Context Navigation

MSFiller.cpp @ 1988

Visit:

Last change on this file since 1988 was 1988, 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...

More speed up filler.

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

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

Context Navigation

source: trunk/src/MSFiller.cpp @ 1988

Download in other formats: