Context Navigation

MSFiller.cpp@ 2020

Visit:

Last change on this file since 2020 was 2019, checked in by Takeshi Nakazato, 14 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...

Bug fix

1) supporting variable shaped array input to some column (SPECTRAL DATA, TSYS, and TCAL)

2) put [1.0] instead of empty array as default TCAL in NROFiller

3) Check if some subtables have any effective rows in MSWriter

File size: 56.1 KB

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

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

Context Navigation

source: trunk/src/MSFiller.cpp@ 2020

Download in other formats: