Context Navigation

MSFiller.cpp@ 2031

Visit:

Last change on this file since 2031 was 2028, 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 on getSrcType() that takes place if ALMA data.

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

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

Context Navigation

source: trunk/src/MSFiller.cpp@ 2031

Download in other formats: