Context Navigation

MSFiller.cpp@ 2023

Visit:

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

Copy ASDM subtable that are optionally created by importasdm task.
For filler, their root names are stored as String in table keyword,
while for writer, existing tables are copied and are stored as Table
in table keyword.

File size: 57.4 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
	121	os_ << "Parsing MS options" << endl ;
	122	os_ << " getPt = " << getPt_ << endl ;
[2021]	123	os_ << " antenna = " << antenna_ << endl ;
	124	os_ << " antennaStr = " << antennaStr_ << LogIO::POST ;
[1974]	125
[1987]	126	MeasurementSet *tmpMS = new MeasurementSet( filename, Table::Old ) ;
	127	//mstable_ = (*tmpMS)( tmpMS->col("ANTENNA1") == antenna_
	128	// && tmpMS->col("ANTENNA1") == tmpMS->col("ANTENNA2") ) ;
	129	tablename_ = tmpMS->tableName() ;
[2021]	130	if ( antenna_ == -1 && antennaStr_.size() > 0 ) {
	131	MSAntennaIndex msAntIdx( tmpMS->antenna() ) ;
	132	Vector<Int> id = msAntIdx.matchAntennaName( antennaStr_ ) ;
	133	if ( id.size() > 0 )
	134	antenna_ = id[0] ;
	135	os_ << "searched antenna_ = " << antenna_ << 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" ) ;
	352	sdh.observer = tcolr->asString( obsId ) ;
	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 ) ;
	569	Double refpix = 0.5*(nchan-1) ;
	570	Double refval = 0.0 ;
[1987]	571	sharedQDArrCol = new ROArrayQuantColumn<Double>( spwtab, "CHAN_WIDTH" ) ;
	572	Double increment = (*sharedQDArrCol)( spwId )( refip ).getValue( "Hz" ) ;
	573	delete sharedQDArrCol ;
[1974]	574	//os_ << "nchan = " << nchan << " refchan = " << refchan << "(even=" << even << ") refpix = " << refpix << LogIO::POST ;
[1987]	575	sharedQDArrCol = new ROArrayQuantColumn<Double>( spwtab, "CHAN_FREQ" ) ;
	576	Vector< Quantum<Double> > chanFreqs = (*sharedQDArrCol)( spwId ) ;
	577	delete sharedQDArrCol ;
[1974]	578	if ( freqRef == MFrequency::LSRK ) {
	579	if ( even ) {
	580	IPosition refip0( 1, refchan-1 ) ;
	581	Double refval0 = chanFreqs(refip0).getValue("Hz") ;
	582	Double refval1 = chanFreqs(refip).getValue("Hz") ;
	583	refval = 0.5 * ( refval0 + refval1 ) ;
	584	}
	585	else {
	586	refval = chanFreqs(refip).getValue("Hz") ;
	587	}
	588	}
	589	else {
	590	MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
	591	if ( even ) {
	592	IPosition refip0( 1, refchan-1 ) ;
	593	Double refval0 = chanFreqs(refip0).getValue("Hz") ;
	594	Double refval1 = chanFreqs(refip).getValue("Hz") ;
	595	refval = 0.5 * ( refval0 + refval1 ) ;
	596	refval = tolsr( refval ).get("Hz").getValue() ;
	597	}
	598	else {
	599	refval = tolsr( chanFreqs(refip) ).get("Hz").getValue() ;
	600	}
	601	}
	602	uInt freqId = table_->frequencies().addEntry( refpix, refval, increment ) ;
	603	if ( ifmap.find( spwId ) == ifmap.end() ) {
	604	ifmap.insert( pair<Int, uInt>(spwId,freqId) ) ;
	605	//os_ << "added to ifmap: (" << spwId << "," << freqId << ")" << LogIO::POST ;
	606	}
[2002]	607
	608	// FREQ_ID
	609	uintRF.attachToRecord( trec, "FREQ_ID" ) ;
	610	*uintRF = freqId ;
	611
[1974]	612	// for TSYS and TCAL
[2004]	613	Vector<MEpoch> scTime ;
	614	Vector<Double> scInterval ;
	615	ROArrayColumn<Float> scTsysCol ;
	616	MSSysCal caltabsel ;
	617	if ( isSysCal_ ) {
	618	caltabsel = caltab( caltab.col("ANTENNA_ID") == antenna_ && caltab.col("FEED_ID") == feedId && caltab.col("SPECTRAL_WINDOW_ID") == spwId ).sort("TIME") ;
	619	ROScalarMeasColumn<MEpoch> scTimeCol( caltabsel, "TIME" ) ;
	620	scTime.resize( caltabsel.nrow() ) ;
	621	for ( uInt irow = 0 ; irow < caltabsel.nrow() ; irow++ )
	622	scTime[irow] = scTimeCol( irow ) ;
	623	ROScalarColumn<Double> *scIntervalCol = new ROScalarColumn<Double>( caltabsel, "INTERVAL" ) ;
	624	scInterval = scIntervalCol->getColumn() ;
	625	delete scIntervalCol ;
	626	scTsysCol.attach( caltabsel, colTsys_ ) ;
	627	}
[2021]	628	// time1 = gettimeofday_sec() ;
	629	// os_ << "end 3rd iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
[1974]	630	//
	631	// ITERATION: SCAN_NUMBER
	632	//
	633	TableIterator iter4( t3, "SCAN_NUMBER" ) ;
	634	while( !iter4.pastEnd() ) {
[2021]	635	// time0 = gettimeofday_sec() ;
	636	// os_ << "start 4th iteration: " << time0 << LogIO::POST ;
[1987]	637	Table t4 = iter4.table() ;
	638	tcolr = tpoolr->construct( t4, "SCAN_NUMBER" ) ;
	639	Int scanNum = tcolr->asInt( 0 ) ;
	640	tpoolr->destroy( tcolr ) ;
[2002]	641
	642	// SCANNO
	643	uintRF.attachToRecord( trec, "SCANNO" ) ;
	644	*uintRF = scanNum - 1 ;
	645
[2021]	646	// time1 = gettimeofday_sec() ;
	647	// os_ << "end 4th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
[1974]	648	//
	649	// ITERATION: STATE_ID
	650	//
	651	TableIterator iter5( t4, "STATE_ID" ) ;
	652	while( !iter5.pastEnd() ) {
[2004]	653	//time0 = gettimeofday_sec() ;
	654	//os_ << "start 5th iteration: " << time0 << LogIO::POST ;
[1987]	655	Table t5 = iter5.table() ;
	656	tcolr = tpoolr->construct( t5, "STATE_ID" ) ;
	657	Int stateId = tcolr->asInt( 0 ) ;
	658	tpoolr->destroy( tcolr ) ;
	659	tcolr = tpoolr->construct( stattab, "OBS_MODE" ) ;
	660	String obstype = tcolr->asString( stateId ) ;
	661	tpoolr->destroy( tcolr ) ;
[1974]	662	if ( sdh.obstype == "" ) sdh.obstype = obstype ;
	663
	664	Int nrow = t5.nrow() ;
[2004]	665	//time1 = gettimeofday_sec() ;
	666	//os_ << "end 5th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
[2002]	667
[2004]	668	uInt cycle = 0 ;
	669
[2002]	670	Cube<Float> spArr ;
	671	Cube<Bool> flArr ;
[1974]	672	if ( isFloatData_ ) {
[2002]	673	ROArrayColumn<Bool> mFlagCol( t5, "FLAG" ) ;
[1974]	674	ROArrayColumn<Float> mFloatDataCol( t5, "FLOAT_DATA" ) ;
[2002]	675	spArr = mFloatDataCol.getColumn() ;
	676	flArr = mFlagCol.getColumn() ;
[1974]	677	if ( sdh.fluxunit == "" ) {
[2002]	678	const TableRecord &dataColKeys = mFloatDataCol.keywordSet() ;
[1974]	679	if ( dataColKeys.isDefined( "UNIT" ) )
	680	sdh.fluxunit = dataColKeys.asString( "UNIT" ) ;
	681	}
	682	}
	683	else if ( isData_ ) {
[2002]	684	spArr.resize( npol, nchan, nrow ) ;
[2004]	685	flArr.resize( npol, nchan, nrow ) ;
[2002]	686	ROArrayColumn<Bool> mFlagCol( t5, "FLAG" ) ;
[1974]	687	ROArrayColumn<Complex> mDataCol( t5, "DATA" ) ;
[1987]	688	for ( Int irow = 0 ; irow < nrow ; irow++ ) {
	689	Bool crossOK = False ;
	690	Matrix<Complex> sp = mDataCol( irow ) ;
[2002]	691	Matrix<Bool> fl = mFlagCol( irow ) ;
[2004]	692	Matrix<Float> spxy = spArr.xyPlane( irow ) ;
	693	Matrix<Bool> flxy = flArr.xyPlane( irow ) ;
[1987]	694	for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
	695	if ( corrtype[ipol] == Stokes::XY \|\| corrtype[ipol] == Stokes::YX
	696	\|\| corrtype[ipol] == Stokes::RL \|\| corrtype[ipol] == Stokes::LR ) {
	697	if ( !crossOK ) {
[2004]	698	spxy.row( ipol ) = real( sp.row( ipol ) ) ;
	699	flxy.row( ipol ) = fl.row( ipol ) ;
[2002]	700	if ( corrtype[ipol] == Stokes::XY \|\| corrtype[ipol] == Stokes::RL ) {
[2004]	701	spxy.row( ipol+1 ) = imag( sp.row( ipol ) ) ;
	702	flxy.row( ipol+1 ) = fl.row( ipol ) ;
[2002]	703	}
	704	else {
[2004]	705	spxy.row( ipol+1 ) = imag( conj( sp.row( ipol ) ) ) ;
	706	flxy.row( ipol+1 ) = fl.row( ipol ) ;
[2002]	707	}
[2004]	708	crossOK = True ;
[1974]	709	}
	710	}
[1987]	711	else {
[2004]	712	spxy.row( ipol ) = real( sp.row( ipol ) ) ;
	713	flxy.row( ipol ) = fl.row( ipol ) ;
[1987]	714	}
[1974]	715	}
[1987]	716	}
[1974]	717	if ( sdh.fluxunit == "" ) {
[2002]	718	const TableRecord &dataColKeys = mDataCol.keywordSet() ;
[1974]	719	if ( dataColKeys.isDefined( "UNIT" ) )
[2002]	720	sdh.fluxunit = dataColKeys.asString( "UNIT" ) ;
	721	}
[1974]	722	}
[1987]	723	ROScalarMeasColumn<MEpoch> *mTimeCol = new ROScalarMeasColumn<MEpoch>( t5, "TIME" ) ;
[2002]	724	Block<MEpoch> mTimeB( nrow ) ;
	725	for ( Int irow = 0 ; irow < nrow ; irow++ )
	726	mTimeB[irow] = (*mTimeCol)( irow ) ;
	727	ROTableColumn *mIntervalCol = tpoolr->construct( t5, "INTERVAL" ) ;
	728	ROTableColumn *mFlagRowCol = tpoolr->construct( t5, "FLAG_ROW" ) ;
	729	Block<Int> sysCalIdx( nrow, -1 ) ;
	730	if ( isSysCal_ ) {
[2017]	731	getSysCalTime( scTime, scInterval, mTimeB, sysCalIdx ) ;
[1974]	732	}
[2002]	733	delete mTimeCol ;
	734	Matrix<Float> defaulttsys( npol, 1, 1.0 ) ;
	735	Int srcType = getSrcType( stateId, tpoolr ) ;
	736	uInt diridx = 0 ;
	737	MDirection::Types dirType ;
	738	uInt wid = 0 ;
	739	Int pidx = 0 ;
	740	Bool crossOK = False ;
	741	Block<uInt> polnos( npol, 99 ) ;
	742	for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
	743	Block<uInt> p = getPolNo( corrtype[ipol] ) ;
	744	if ( p.size() > 1 ) {
	745	if ( crossOK ) continue ;
	746	else {
	747	polnos[pidx] = p[0] ;
	748	pidx++ ;
	749	polnos[pidx] = p[1] ;
	750	pidx++ ;
	751	crossOK = True ;
[1974]	752	}
	753	}
[2002]	754	else {
	755	polnos[pidx] = p[0] ;
	756	pidx++ ;
	757	}
[1974]	758	}
[2002]	759
	760	// SRCTYPE
	761	intRF.attachToRecord( trec, "SRCTYPE" ) ;
	762	*intRF = srcType ;
[1974]	763
[2002]	764	for ( Int irow = 0 ; irow < nrow ; irow++ ) {
	765	// CYCLENO
	766	*cycleRF = cycle ;
[1987]	767
[2002]	768	// FLAGROW
	769	*flrRF = (uInt)mFlagRowCol->asBool( irow ) ;
	770
	771	// SPECTRA, FLAG
	772	Matrix<Float> sp = spArr.xyPlane( irow ) ;
	773	Matrix<Bool> flb = flArr.xyPlane( irow ) ;
	774	Matrix<uChar> fl( flb.shape() ) ;
	775	convertArray( fl, flb ) ;
	776
	777	// TIME
[2004]	778	*timeRF = mTimeB[irow].get("d").getValue() ;
[2002]	779
	780	// INTERVAL
	781	*intervalRF = (Double)(mIntervalCol->asdouble( irow )) ;
	782
	783	// TSYS
	784	Matrix<Float> tsys ;
	785	if ( sysCalIdx[irow] != -1 )
	786	tsys = scTsysCol( irow ) ;
	787	else
	788	tsys = defaulttsys ;
	789
	790	// TCAL_ID
[2004]	791	Block<uInt> tcalids( npol, 0 ) ;
	792	if ( sysCalIdx[irow] != -1 ) {
	793	tcalids = getTcalId( feedId, spwId, scTime[sysCalIdx[irow]] ) ;
[1987]	794	}
	795
[2002]	796	// WEATHER_ID
	797	if ( isWeather_ )
	798	wid = getWeatherId( wid, mTimeB[irow].get("s").getValue() ) ;
	799	*widRF = wid ;
	800
[1987]	801
[2002]	802	// DIRECTION, AZEL, SCANRATE
	803	if ( getPt_ ) {
[1974]	804	Vector<Double> dir ;
	805	Vector<Double> scanrate ;
	806	String refString ;
[2002]	807	diridx = getDirection( diridx, dir, scanrate, refString, pointtab, mTimeB[irow].get("s").getValue() ) ;
[1974]	808	MDirection::getType( dirType, refString ) ;
[2002]	809	mf.resetEpoch( mTimeB[irow] ) ;
[1974]	810	mf.resetDirection( MDirection( MVDirection(dir), dirType ) ) ;
	811	if ( refString == "J2000" ) {
[2002]	812	*dirRF = dir ;
[1974]	813	MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
	814	Vector<Double> azel = toazel( dir ).getAngle("rad").getValue() ;
[2002]	815	*azRF = (Float)azel(0) ;
	816	*elRF = (Float)azel(1) ;
[1974]	817	}
	818	else if ( refString(0,4) == "AZEL" ) {
[2002]	819	*azRF = (Float)dir(0) ;
	820	*elRF = (Float)dir(1) ;
[1974]	821	MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
	822	Vector<Double> newdir = toj2000( dir ).getAngle("rad").getValue() ;
[2002]	823	*dirRF = newdir ;
[1974]	824	}
	825	else {
	826	MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
	827	Vector<Double> azel = toazel( dir ).getAngle("rad").getValue() ;
	828	MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
	829	Vector<Double> newdir = toj2000( dir ).getAngle("rad").getValue() ;
[2002]	830	*dirRF = newdir ;
	831	*azRF = (Float)azel(0) ;
	832	*elRF = (Float)azel(1) ;
[1974]	833	}
	834	if ( scanrate.size() != 0 ) {
[2002]	835	*scrRF = scanrate ;
[1974]	836	}
	837	else {
[2002]	838	*scrRF = defaultScanrate ;
[1974]	839	}
	840	}
[2002]	841	else {
	842	String ref = md.getRefString() ;
[2004]	843	//Vector<Double> defaultDir = srcDir ;
[2002]	844	MDirection::getType( dirType, "J2000" ) ;
	845	if ( ref != "J2000" ) {
	846	ROScalarMeasColumn<MEpoch> tmCol( pointtab, "TIME" ) ;
	847	mf.resetEpoch( tmCol( 0 ) ) ;
	848	MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
	849	defaultDir = toj2000( defaultDir ).getAngle("rad").getValue() ;
[1974]	850	}
[2002]	851	mf.resetEpoch( mTimeB[irow] ) ;
	852	MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
	853	Vector<Double> azel = toazel( defaultDir ).getAngle("rad").getValue() ;
	854	*azRF = (Float)azel(0) ;
	855	*elRF = (Float)azel(1) ;
	856	*dirRF = defaultDir ;
	857	*scrRF = defaultScanrate ;
[1974]	858	}
	859
[2002]	860	// Polarization dependent things
	861	for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
	862	// POLNO
	863	*polnoRF = polnos[ipol] ;
[1987]	864
[2019]	865	//*spRF = sp.row( ipol ) ;
	866	//*ucarrRF = fl.row( ipol ) ;
	867	//*tsysRF = tsys.row( ipol ) ;
	868	spRF.define( sp.row( ipol ) ) ;
	869	ucarrRF.define( fl.row( ipol ) ) ;
	870	tsysRF.define( tsys.row( ipol ) ) ;
[2002]	871	*tcalidRF = tcalids[ipol] ;
[1987]	872
[2002]	873	// Commit row
	874	stab.addRow() ;
	875	row.put( stab.nrow()-1 ) ;
[1987]	876	}
[2004]	877
	878	cycle++ ;
[1987]	879	}
[2002]	880	tpoolr->destroy( mIntervalCol ) ;
	881	tpoolr->destroy( mFlagRowCol ) ;
[1987]	882
[2004]	883	//time1 = gettimeofday_sec() ;
	884	//os_ << "end 5th iteration: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
[1974]	885
	886	iter5.next() ;
	887	}
[1978]	888
	889
[1974]	890	iter4.next() ;
	891	}
[1978]	892
	893
[1974]	894	iter3.next() ;
	895	}
[1978]	896
[1990]	897
[1974]	898	iter2.next() ;
	899	}
[1978]	900
	901
[1974]	902	iter1.next() ;
	903	}
	904	if ( sdh.nbeam < nbeam ) sdh.nbeam = nbeam ;
[1978]	905
[1974]	906	iter0.next() ;
	907	}
	908
[1987]	909
	910	delete tpoolr ;
	911
	912
[1978]	913	// Table Keywords
[2004]	914	// TODO: This must be updated if writer will be replaced from STWriter to MSWriter
	915	// sdh.nif = ifmap.size() ;
	916	sdh.nif++ ;
[1987]	917	if ( ( telescopeName == "" ) \|\| ( antennaName == telescopeName ) ) {
[1974]	918	sdh.antennaname = antennaName ;
	919	}
	920	else {
	921	sdh.antennaname = telescopeName + "//" + antennaName ;
	922	}
	923	if ( stationName != "" ) {
	924	sdh.antennaname += "@" + stationName ;
	925	}
[1987]	926	ROArrayColumn<Double> pdirCol( pointtab, "DIRECTION" ) ;
	927	String dirref = pdirCol.keywordSet().asRecord("MEASINFO").asString("Ref") ;
[1974]	928	if ( dirref == "AZELGEO" \|\| dirref == "AZEL" ) {
	929	dirref = "J2000" ;
	930	}
	931	sscanf( dirref.chars()+1, "%f", &sdh.equinox ) ;
	932	sdh.epoch = "UTC" ;
	933	if (sdh.freqref == "TOPO") {
	934	sdh.freqref = "TOPOCENT";
	935	} else if (sdh.freqref == "GEO") {
	936	sdh.freqref = "GEOCENTR";
	937	} else if (sdh.freqref == "BARY") {
	938	sdh.freqref = "BARYCENT";
	939	} else if (sdh.freqref == "GALACTO") {
	940	sdh.freqref = "GALACTOC";
	941	} else if (sdh.freqref == "LGROUP") {
	942	sdh.freqref = "LOCALGRP";
	943	} else if (sdh.freqref == "CMB") {
	944	sdh.freqref = "CMBDIPOL";
	945	} else if (sdh.freqref == "REST") {
	946	sdh.freqref = "SOURCE";
	947	}
	948	table_->setHeader( sdh ) ;
[1987]	949
	950	// save path to POINTING table
[2022]	951	// 2011/3/2 TN
	952	// So far, path to original POINTING table is always stored
	953	// since sd tasks and regressions don't support getpt control
	954	//if ( !getPt_ ) {
[1987]	955	Path datapath( tablename_ ) ;
	956	String pTabName = datapath.absoluteName() + "/POINTING" ;
[1990]	957	stab.rwKeywordSet().define( "POINTING", pTabName ) ;
[2022]	958	//}
[1987]	959
	960	// for GBT
[2022]	961	if ( antennaName.matches( "GBT" ) ) {
[1987]	962	String goTabName = datapath.absoluteName() + "/GBT_GO" ;
[1990]	963	stab.rwKeywordSet().define( "GBT_GO", goTabName ) ;
[1987]	964	}
[2022]	965
	966	// for MS created from ASDM
	967	//mstable_.keywordSet().print(cout) ;
	968	const TableRecord &msKeys = mstable_.keywordSet() ;
	969	uInt nfields = msKeys.nfields() ;
	970	for ( uInt ifield = 0 ; ifield < nfields ; ifield++ ) {
	971	String name = msKeys.name( ifield ) ;
	972	//os_ << "name = " << name << LogIO::POST ;
	973	if ( name.find( "ASDM" ) != String::npos ) {
	974	String asdmpath = msKeys.asTable( ifield ).tableName() ;
	975	os_ << "ASDM table: " << asdmpath << LogIO::POST ;
	976	stab.rwKeywordSet().define( name, asdmpath ) ;
	977	}
	978	}
	979
[2004]	980	// double endSec = gettimeofday_sec() ;
	981	// os_ << "end MSFiller::fill() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
[1974]	982	}
	983
	984	void MSFiller::close()
	985	{
[1987]	986	//tablesel_.closeSubTables() ;
[1974]	987	mstable_.closeSubTables() ;
[1987]	988	//tablesel_.unlock() ;
[1974]	989	mstable_.unlock() ;
	990	}
	991
[1987]	992	Int MSFiller::getSrcType( Int stateId, boost::object_pool<ROTableColumn> *tpool )
[1974]	993	{
[2004]	994	// double startSec = gettimeofday_sec() ;
	995	// os_ << "start MSFiller::getSrcType() startSec=" << startSec << LogIO::POST ;
[1974]	996
	997	MSState statetab = mstable_.state() ;
[1987]	998	ROTableColumn *sharedCol ;
	999	sharedCol = tpool->construct( statetab, "OBS_MODE" ) ;
	1000	String obsMode = sharedCol->asString( stateId ) ;
	1001	tpool->destroy( sharedCol ) ;
	1002	sharedCol = tpool->construct( statetab, "SIG" ) ;
	1003	Bool sig = sharedCol->asBool( stateId ) ;
	1004	tpool->destroy( sharedCol ) ;
	1005	sharedCol = tpool->construct( statetab, "REF" ) ;
	1006	Bool ref = sharedCol->asBool( stateId ) ;
	1007	tpool->destroy( sharedCol ) ;
	1008	sharedCol = tpool->construct( statetab, "CAL" ) ;
	1009	Double cal = (Double)(sharedCol->asdouble( stateId )) ;
	1010	tpool->destroy( sharedCol ) ;
[1974]	1011	//os_ << "OBS_MODE = " << obsMode << LogIO::POST ;
	1012
	1013	// determine separator
	1014	String sep = "" ;
	1015	if ( obsMode.find( ":" ) != String::npos ) {
	1016	sep = ":" ;
	1017	}
	1018	else if ( obsMode.find( "." ) != String::npos ) {
	1019	sep = "." ;
	1020	}
	1021
	1022	// determine SRCTYPE
	1023	Int srcType = SrcType::NOTYPE ;
	1024	if ( sep == ":" ) {
	1025	// sep == ":"
	1026	//
	1027	// GBT case
	1028	//
	1029	// obsMode1=Nod
	1030	// NOD
	1031	// obsMode1=OffOn
	1032	// obsMode2=PSWITCHON: PSON
	1033	// obsMode2=PSWITCHOFF: PSOFF
	1034	// obsMode1=??
	1035	// obsMode2=FSWITCH:
	1036	// SIG=1: FSON
	1037	// REF=1: FSOFF
[1987]	1038	// Calibration scan if CAL != 0
[1974]	1039	Int epos = obsMode.find_first_of( sep ) ;
	1040	Int nextpos = obsMode.find_first_of( sep, epos+1 ) ;
	1041	String obsMode1 = obsMode.substr( 0, epos ) ;
	1042	String obsMode2 = obsMode.substr( epos+1, nextpos-epos-1 ) ;
	1043	if ( obsMode1 == "Nod" ) {
	1044	srcType = SrcType::NOD ;
	1045	}
	1046	else if ( obsMode1 == "OffOn" ) {
	1047	if ( obsMode2 == "PSWITCHON" ) srcType = SrcType::PSON ;
	1048	if ( obsMode2 == "PSWITCHOFF" ) srcType = SrcType::PSOFF ;
	1049	}
	1050	else {
	1051	if ( obsMode2 == "FSWITCH" ) {
	1052	if ( sig ) srcType = SrcType::FSON ;
	1053	if ( ref ) srcType = SrcType::FSOFF ;
	1054	}
	1055	}
[1987]	1056	if ( cal > 0.0 ) {
	1057	if ( srcType == SrcType::NOD )
	1058	srcType = SrcType::NODCAL ;
	1059	else if ( srcType == SrcType::PSON )
	1060	srcType = SrcType::PONCAL ;
	1061	else if ( srcType == SrcType::PSOFF )
	1062	srcType = SrcType::POFFCAL ;
	1063	else if ( srcType == SrcType::FSON )
	1064	srcType = SrcType::FONCAL ;
	1065	else if ( srcType == SrcType::FSOFF )
	1066	srcType = SrcType::FOFFCAL ;
	1067	else
	1068	srcType = SrcType::CAL ;
	1069	}
[1974]	1070	}
	1071	else if ( sep == "." ) {
	1072	// sep == "."
	1073	//
	1074	// ALMA & EVLA case (MS via ASDM)
	1075	//
	1076	// obsMode1=CALIBRATE_*
	1077	// obsMode2=ON_SOURCE: PONCAL
	1078	// obsMode2=OFF_SOURCE: POFFCAL
	1079	// obsMode1=OBSERVE_TARGET
	1080	// obsMode2=ON_SOURCE: PON
	1081	// obsMode2=OFF_SOURCE: POFF
	1082	Int epos = obsMode.find_first_of( sep ) ;
	1083	Int nextpos = obsMode.find_first_of( sep, epos+1 ) ;
	1084	String obsMode1 = obsMode.substr( 0, epos ) ;
	1085	String obsMode2 = obsMode.substr( epos+1, nextpos-epos-1 ) ;
	1086	if ( obsMode1.find( "CALIBRATE_" ) == 0 ) {
	1087	if ( obsMode2 == "ON_SOURCE" ) srcType = SrcType::PONCAL ;
	1088	if ( obsMode2 == "OFF_SOURCE" ) srcType = SrcType::POFFCAL ;
	1089	}
	1090	else if ( obsMode1 == "OBSERVE_TARGET" ) {
	1091	if ( obsMode2 == "ON_SOURCE" ) srcType = SrcType::PSON ;
	1092	if ( obsMode2 == "OFF_SOURCE" ) srcType = SrcType::PSOFF ;
	1093	}
	1094	}
	1095	else {
	1096	if ( sig ) srcType = SrcType::SIG ;
	1097	if ( ref ) srcType = SrcType::REF ;
	1098	}
	1099
	1100	//os_ << "srcType = " << srcType << LogIO::POST ;
[2004]	1101	// double endSec = gettimeofday_sec() ;
	1102	// os_ << "end MSFiller::getSrcType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
[1974]	1103	return srcType ;
	1104	}
	1105
[2002]	1106	//Vector<uInt> MSFiller::getPolNo( Int corrType )
	1107	Block<uInt> MSFiller::getPolNo( Int corrType )
[1974]	1108	{
[2004]	1109	// double startSec = gettimeofday_sec() ;
	1110	// os_ << "start MSFiller::getPolNo() startSec=" << startSec << LogIO::POST ;
[2002]	1111	Block<uInt> polno( 1 ) ;
[1974]	1112
	1113	if ( corrType == Stokes::I \|\| corrType == Stokes::RR \|\| corrType == Stokes::XX ) {
	1114	polno = 0 ;
	1115	}
	1116	else if ( corrType == Stokes::Q \|\| corrType == Stokes::LL \|\| corrType == Stokes::YY ) {
	1117	polno = 1 ;
	1118	}
	1119	else if ( corrType == Stokes::U ) {
	1120	polno = 2 ;
	1121	}
	1122	else if ( corrType == Stokes::V ) {
	1123	polno = 3 ;
	1124	}
	1125	else if ( corrType == Stokes::RL \|\| corrType == Stokes::XY \|\| corrType == Stokes::LR \|\| corrType == Stokes::RL ) {
	1126	polno.resize( 2 ) ;
	1127	polno[0] = 2 ;
	1128	polno[1] = 3 ;
	1129	}
	1130	else if ( corrType == Stokes::Plinear ) {
	1131	polno[0] = 1 ;
	1132	}
	1133	else if ( corrType == Stokes::Pangle ) {
	1134	polno[0] = 2 ;
	1135	}
	1136	else {
	1137	polno = 99 ;
	1138	}
	1139	//os_ << "polno = " << polno << LogIO::POST ;
[2004]	1140	// double endSec = gettimeofday_sec() ;
	1141	// os_ << "end MSFiller::getPolNo() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
[1974]	1142
	1143	return polno ;
	1144	}
	1145
	1146	String MSFiller::getPolType( Int corrType )
	1147	{
[2004]	1148	// double startSec = gettimeofday_sec() ;
	1149	// os_ << "start MSFiller::getPolType() startSec=" << startSec << LogIO::POST ;
[1974]	1150	String poltype = "" ;
	1151
	1152	if ( corrType == Stokes::I \|\| corrType == Stokes::Q \|\| corrType == Stokes::U \|\| corrType == Stokes::V )
	1153	poltype = "stokes" ;
	1154	else if ( corrType == Stokes::XX \|\| corrType == Stokes::YY \|\| corrType == Stokes::XY \|\| corrType == Stokes::YX )
	1155	poltype = "linear" ;
	1156	else if ( corrType == Stokes::RR \|\| corrType == Stokes::LL \|\| corrType == Stokes::RL \|\| corrType == Stokes::LR )
	1157	poltype = "circular" ;
	1158	else if ( corrType == Stokes::Plinear \|\| corrType == Stokes::Pangle )
	1159	poltype = "linpol" ;
	1160
[2004]	1161	// double endSec = gettimeofday_sec() ;
	1162	// os_ << "end MSFiller::getPolType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
[1974]	1163	return poltype ;
	1164	}
	1165
	1166	void MSFiller::fillWeather()
	1167	{
[2004]	1168	// double startSec = gettimeofday_sec() ;
	1169	// os_ << "start MSFiller::fillWeather() startSec=" << startSec << LogIO::POST ;
	1170
	1171	if ( !isWeather_ ) {
	1172	// add dummy row
	1173	table_->weather().table().addRow(1,True) ;
	1174	return ;
	1175	}
	1176
[1987]	1177	Table mWeather = mstable_.weather() ;
	1178	//Table mWeatherSel = mWeather( mWeather.col("ANTENNA_ID") == antenna_ ).sort("TIME") ;
	1179	Table mWeatherSel( mWeather( mWeather.col("ANTENNA_ID") == antenna_ ).sort("TIME") ) ;
[1974]	1180	//os_ << "mWeatherSel.nrow() = " << mWeatherSel.nrow() << LogIO::POST ;
	1181	if ( mWeatherSel.nrow() == 0 ) {
	1182	os_ << "No rows with ANTENNA_ID = " << antenna_ << ", Try -1..." << LogIO::POST ;
[1987]	1183	mWeatherSel = Table( MSWeather( mWeather( mWeather.col("ANTENNA_ID") == -1 ) ) ) ;
[1974]	1184	if ( mWeatherSel.nrow() == 0 ) {
	1185	os_ << "No rows in WEATHER table" << LogIO::POST ;
	1186	}
	1187	}
[1987]	1188	uInt wnrow = mWeatherSel.nrow() ;
[1974]	1189	//os_ << "wnrow = " << wnrow << LogIO::POST ;
	1190
	1191	if ( wnrow == 0 )
	1192	return ;
	1193
	1194	Table wtab = table_->weather().table() ;
	1195	wtab.addRow( wnrow ) ;
	1196
[1987]	1197	ScalarColumn<Float> *fCol ;
	1198	ROScalarColumn<Float> *sharedFloatCol ;
	1199	if ( mWeatherSel.tableDesc().isColumn( "TEMPERATURE" ) ) {
	1200	fCol = new ScalarColumn<Float>( wtab, "TEMPERATURE" ) ;
	1201	sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "TEMPERATURE" ) ;
	1202	fCol->putColumn( *sharedFloatCol ) ;
	1203	delete sharedFloatCol ;
	1204	delete fCol ;
	1205	}
	1206	if ( mWeatherSel.tableDesc().isColumn( "PRESSURE" ) ) {
	1207	fCol = new ScalarColumn<Float>( wtab, "PRESSURE" ) ;
	1208	sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "PRESSURE" ) ;
	1209	fCol->putColumn( *sharedFloatCol ) ;
	1210	delete sharedFloatCol ;
	1211	delete fCol ;
	1212	}
	1213	if ( mWeatherSel.tableDesc().isColumn( "REL_HUMIDITY" ) ) {
	1214	fCol = new ScalarColumn<Float>( wtab, "HUMIDITY" ) ;
	1215	sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "REL_HUMIDITY" ) ;
	1216	fCol->putColumn( *sharedFloatCol ) ;
	1217	delete sharedFloatCol ;
	1218	delete fCol ;
	1219	}
	1220	if ( mWeatherSel.tableDesc().isColumn( "WIND_SPEED" ) ) {
	1221	fCol = new ScalarColumn<Float>( wtab, "WINDSPEED" ) ;
	1222	sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "WIND_SPEED" ) ;
	1223	fCol->putColumn( *sharedFloatCol ) ;
	1224	delete sharedFloatCol ;
	1225	delete fCol ;
	1226	}
	1227	if ( mWeatherSel.tableDesc().isColumn( "WIND_DIRECTION" ) ) {
	1228	fCol = new ScalarColumn<Float>( wtab, "WINDAZ" ) ;
	1229	sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "WIND_DIRECTION" ) ;
	1230	fCol->putColumn( *sharedFloatCol ) ;
	1231	delete sharedFloatCol ;
	1232	delete fCol ;
	1233	}
[1974]	1234	ScalarColumn<uInt> idCol( wtab, "ID" ) ;
[1987]	1235	for ( uInt irow = 0 ; irow < wnrow ; irow++ )
	1236	idCol.put( irow, irow ) ;
[1974]	1237
[1987]	1238	ROScalarQuantColumn<Double> tqCol( mWeatherSel, "TIME" ) ;
	1239	ROScalarColumn<Double> tCol( mWeatherSel, "TIME" ) ;
	1240	String tUnit = tqCol.getUnits() ;
	1241	mwTime_ = tCol.getColumn() ;
[1974]	1242	if ( tUnit == "d" )
	1243	mwTime_ *= 86400.0 ;
[1987]	1244	tqCol.attach( mWeatherSel, "INTERVAL" ) ;
	1245	tCol.attach( mWeatherSel, "INTERVAL" ) ;
	1246	String iUnit = tqCol.getUnits() ;
	1247	mwInterval_ = tCol.getColumn() ;
[1974]	1248	if ( iUnit == "d" )
	1249	mwInterval_ *= 86400.0 ;
	1250	//os_ << "mwTime[0] = " << mwTime_[0] << " mwInterval[0] = " << mwInterval_[0] << LogIO::POST ;
[2004]	1251	// double endSec = gettimeofday_sec() ;
	1252	// os_ << "end MSFiller::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
[1974]	1253	}
	1254
	1255	void MSFiller::fillFocus()
	1256	{
[2004]	1257	// double startSec = gettimeofday_sec() ;
	1258	// os_ << "start MSFiller::fillFocus() startSec=" << startSec << LogIO::POST ;
[1974]	1259	// tentative
	1260	Table tab = table_->focus().table() ;
	1261	tab.addRow( 1 ) ;
	1262	ScalarColumn<uInt> idCol( tab, "ID" ) ;
	1263	idCol.put( 0, 0 ) ;
[2004]	1264	// double endSec = gettimeofday_sec() ;
	1265	// os_ << "end MSFiller::fillFocus() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
[1974]	1266	}
	1267
[2002]	1268	void MSFiller::fillTcal( boost::object_pool<ROTableColumn> *tpoolr )
[1974]	1269	{
[2004]	1270	// double startSec = gettimeofday_sec() ;
	1271	// os_ << "start MSFiller::fillTcal() startSec=" << startSec << LogIO::POST ;
[1987]	1272
[2004]	1273	if ( !isSysCal_ ) {
	1274	// add dummy row
	1275	os_ << "No SysCal rows" << LogIO::POST ;
	1276	table_->tcal().table().addRow(1,True) ;
[2021]	1277	Vector<Float> defaultTcal( 1, 1.0 ) ;
	1278	ArrayColumn<Float> tcalCol( table_->tcal().table(), "TCAL" ) ;
	1279	tcalCol.put( 0, defaultTcal ) ;
[2004]	1280	return ;
	1281	}
	1282
[1987]	1283	Table sctab = mstable_.sysCal() ;
[1974]	1284	if ( sctab.nrow() == 0 ) {
	1285	os_ << "No SysCal rows" << LogIO::POST ;
	1286	return ;
	1287	}
[1987]	1288	Table sctabsel( sctab( sctab.col("ANTENNA_ID") == antenna_ ) ) ;
[1974]	1289	if ( sctabsel.nrow() == 0 ) {
	1290	os_ << "No SysCal rows" << LogIO::POST ;
	1291	return ;
	1292	}
[1987]	1293	ROArrayColumn<Float> *tmpTcalCol = new ROArrayColumn<Float>( sctabsel, "TCAL" ) ;
	1294	uInt npol = tmpTcalCol->shape( 0 )(0) ;
	1295	delete tmpTcalCol ;
[1974]	1296	//os_ << "fillTcal(): npol = " << npol << LogIO::POST ;
	1297	Table tab = table_->tcal().table() ;
	1298	ArrayColumn<Float> tcalCol( tab, "TCAL" ) ;
[1987]	1299	ROTableColumn *sharedCol ;
[1974]	1300	uInt oldnr = 0 ;
	1301	uInt newnr = 0 ;
[2002]	1302	TableRow row( tab ) ;
	1303	TableRecord &trec = row.record() ;
	1304	RecordFieldPtr<uInt> idRF( trec, "ID" ) ;
	1305	RecordFieldPtr<String> timeRF( trec, "TIME" ) ;
	1306	RecordFieldPtr< Array<Float> > tcalRF( trec, "TCAL" ) ;
[1974]	1307	TableIterator iter0( sctabsel, "FEED_ID" ) ;
	1308	while( !iter0.pastEnd() ) {
[1987]	1309	Table t0 = iter0.table() ;
[1988]	1310	sharedCol = tpoolr->construct( t0, "FEED_ID" ) ;
[1987]	1311	Int feedId = sharedCol->asInt( 0 ) ;
[1988]	1312	tpoolr->destroy( sharedCol ) ;
[1974]	1313	TableIterator iter1( t0, "SPECTRAL_WINDOW_ID" ) ;
	1314	while( !iter1.pastEnd() ) {
[1988]	1315	Table t1 = iter1.table() ;
	1316	sharedCol = tpoolr->construct( t1, "SPECTRAL_WINDOW_ID" ) ;
[1987]	1317	Int spwId = sharedCol->asInt( 0 ) ;
[1988]	1318	tpoolr->destroy( sharedCol ) ;
[2002]	1319	tmpTcalCol = new ROArrayColumn<Float>( t1, colTcal_ ) ;
	1320	ROScalarQuantColumn<Double> scTimeCol( t1, "TIME" ) ;
	1321	Vector<uInt> idminmax( 2, oldnr ) ;
	1322	for ( uInt irow = 0 ; irow < t1.nrow() ; irow++ ) {
	1323	String sTime = MVTime( scTimeCol(irow) ).string( MVTime::YMD ) ;
	1324	*timeRF = sTime ;
[1988]	1325	uInt idx = oldnr ;
[2002]	1326	Matrix<Float> subtcal = (*tmpTcalCol)( irow ) ;
[1988]	1327	for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
[2002]	1328	*idRF = idx++ ;
[2019]	1329	//*tcalRF = subtcal.row( ipol ) ;
	1330	tcalRF.define( subtcal.row( ipol ) ) ;
[2002]	1331
	1332	// commit row
	1333	tab.addRow() ;
	1334	row.put( tab.nrow()-1 ) ;
	1335
	1336	newnr++ ;
[1974]	1337	}
[2002]	1338	idminmax[0] = oldnr ;
[1988]	1339	idminmax[1] = newnr - 1 ;
	1340	oldnr = newnr ;
	1341
[1990]	1342	String key = keyTcal( feedId, spwId, sTime ) ;
[1988]	1343	tcalrec_.define( key, idminmax ) ;
[1974]	1344	}
[2002]	1345	delete tmpTcalCol ;
[1974]	1346	iter1++ ;
	1347	}
	1348	iter0++ ;
	1349	}
	1350
	1351	//tcalrec_.print( std::cout ) ;
[2004]	1352	// double endSec = gettimeofday_sec() ;
	1353	// os_ << "end MSFiller::fillTcal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
[1974]	1354	}
	1355
	1356	uInt MSFiller::getWeatherId( uInt idx, Double wtime )
	1357	{
[2004]	1358	// double startSec = gettimeofday_sec() ;
	1359	// os_ << "start MSFiller::getWeatherId() startSec=" << startSec << LogIO::POST ;
[1974]	1360	uInt nrow = mwTime_.size() ;
[2019]	1361	if ( nrow < 2 )
[1974]	1362	return 0 ;
	1363	uInt wid = nrow ;
	1364	for ( uInt i = idx ; i < nrow-1 ; i++ ) {
	1365	Double tStart = mwTime_[i]-0.5*mwInterval_[i] ;
	1366	// use of INTERVAL column is problematic
	1367	// since there are "blank" time of weather monitoring
	1368	//Double tEnd = tStart + mwInterval_[i] ;
	1369	Double tEnd = mwTime_[i+1]-0.5*mwInterval_[i+1] ;
	1370	//os_ << "tStart = " << tStart << " dtEnd = " << tEnd-tStart << " dwtime = " << wtime-tStart << LogIO::POST ;
	1371	if ( wtime >= tStart && wtime <= tEnd ) {
	1372	wid = i ;
	1373	break ;
	1374	}
	1375	}
	1376	if ( wid == nrow ) {
	1377	uInt i = nrow - 1 ;
	1378	Double tStart = mwTime_[i-1]+0.5*mwInterval_[i-1] ;
	1379	Double tEnd = mwTime_[i]+0.5*mwInterval_[i] ;
	1380	//os_ << "tStart = " << tStart << " dtEnd = " << tEnd-tStart << " dwtime = " << wtime-tStart << LogIO::POST ;
	1381	if ( wtime >= tStart && wtime <= tEnd )
[2021]	1382	wid = i-1 ;
	1383	else
[1974]	1384	wid = i ;
	1385	}
	1386
	1387	//if ( wid == nrow )
	1388	//os_ << LogIO::WARN << "Couldn't find correct WEATHER_ID for time " << wtime << LogIO::POST ;
	1389
[2004]	1390	// double endSec = gettimeofday_sec() ;
	1391	// os_ << "end MSFiller::getWeatherId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
[1974]	1392	return wid ;
	1393	}
	1394
[2017]	1395	void MSFiller::getSysCalTime( Vector<MEpoch> &scTime, Vector<Double> &scInterval, Block<MEpoch> &tcol, Block<Int> &tidx )
[1974]	1396	{
[2004]	1397	// double startSec = gettimeofday_sec() ;
	1398	// os_ << "start MSFiller::getSysCalTime() startSec=" << startSec << LogIO::POST ;
	1399
	1400	if ( !isSysCal_ )
	1401	return ;
	1402
[2017]	1403	uInt nrow = tidx.nelements() ;
[2002]	1404	if ( scTime.nelements() == 0 )
[1993]	1405	return ;
[2019]	1406	else if ( scTime.nelements() == 1 ) {
	1407	tidx[0] = 0 ;
	1408	return ;
	1409	}
[2002]	1410	uInt scnrow = scTime.nelements() ;
[1974]	1411	uInt idx = 0 ;
	1412	const Double half = 0.5e0 ;
[2003]	1413	// execute binary search
	1414	idx = binarySearch( scTime, tcol[0].get( "s" ).getValue() ) ;
	1415	if ( idx != 0 )
	1416	idx -= 1 ;
[1974]	1417	for ( uInt i = 0 ; i < nrow ; i++ ) {
[2002]	1418	Double t = tcol[i].get( "s" ).getValue() ;
[2017]	1419	Double tsc = scTime[0].get( "s" ).getValue() ;
	1420	if ( t < tsc ) {
	1421	tidx[i] = 0 ;
	1422	continue ;
	1423	}
[1974]	1424	for ( uInt j = idx ; j < scnrow-1 ; j++ ) {
[2002]	1425	Double tsc1 = scTime[j].get( "s" ).getValue() ;
	1426	Double dt1 = scInterval[j] ;
	1427	Double tsc2 = scTime[j+1].get( "s" ).getValue() ;
	1428	Double dt2 = scInterval[j+1] ;
[1974]	1429	if ( t > tsc1-halfdt1 && t <= tsc2-halfdt2 ) {
[2002]	1430	tidx[i] = j ;
[1974]	1431	idx = j ;
	1432	break ;
	1433	}
	1434	}
[2017]	1435	if ( tidx[i] == -1 ) {
[2022]	1436	// Double tsc = scTime[scnrow-1].get( "s" ).getValue() ;
	1437	// Double dt = scInterval[scnrow-1] ;
[2017]	1438	// if ( t <= tsc+0.5*dt ) {
	1439	// tidx[i] = scnrow-1 ;
	1440	// }
	1441	tidx[i] = scnrow-1 ;
[1974]	1442	}
	1443	}
[2004]	1444	// double endSec = gettimeofday_sec() ;
	1445	// os_ << "end MSFiller::getSysCalTime() endSec=" << endSec << " (" << endSec-startSec << "sec) scnrow = " << scnrow << " tcol.nelements = " << tcol.nelements() << LogIO::POST ;
[1993]	1446	return ;
[1974]	1447	}
	1448
[2004]	1449	Block<uInt> MSFiller::getTcalId( Int fid, Int spwid, MEpoch &t )
	1450	{
[2002]	1451	// double startSec = gettimeofday_sec() ;
[2004]	1452	// os_ << "start MSFiller::getTcalId() startSec=" << startSec << LogIO::POST ;
	1453	//if ( table_->tcal().table().nrow() == 0 ) {
	1454	if ( !isSysCal_ ) {
[1988]	1455	os_ << "No TCAL rows" << LogIO::POST ;
[2002]	1456	Block<uInt> tcalids( 0 ) ;
[1988]	1457	return tcalids ;
	1458	}
[2004]	1459	//String sctime = MVTime( Quantum<Double>(t,"s") ).string(MVTime::YMD) ;
	1460	String sctime = MVTime( t.getValue() ).string(MVTime::YMD) ;
[1990]	1461	String key = keyTcal( fid, spwid, sctime ) ;
[1988]	1462	if ( !tcalrec_.isDefined( key ) ) {
[1974]	1463	os_ << "No TCAL rows" << LogIO::POST ;
[2002]	1464	Block<uInt> tcalids( 0 ) ;
[1974]	1465	return tcalids ;
	1466	}
[1988]	1467	Vector<uInt> ids = tcalrec_.asArrayuInt( key ) ;
	1468	uInt npol = ids[1] - ids[0] + 1 ;
[2002]	1469	Block<uInt> tcalids( npol ) ;
[1988]	1470	tcalids[0] = ids[0] ;
	1471	tcalids[1] = ids[1] ;
	1472	for ( uInt ipol = 2 ; ipol < npol ; ipol++ )
	1473	tcalids[ipol] = ids[0] + ipol - 1 ;
	1474
[2004]	1475	// double endSec = gettimeofday_sec() ;
	1476	// os_ << "end MSFiller::getTcalId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
[1974]	1477	return tcalids ;
	1478	}
	1479
[1987]	1480	uInt MSFiller::getDirection( uInt idx, Vector<Double> &dir, Vector<Double> &srate, String &ref, MSPointing &tab, Double t )
[1974]	1481	{
[2004]	1482	// double startSec = gettimeofday_sec() ;
	1483	// os_ << "start MSFiller::getDirection() startSec=" << startSec << LogIO::POST ;
[1974]	1484	// assume that cols is sorted by TIME
	1485	Bool doInterp = False ;
[1987]	1486	//uInt nrow = cols.nrow() ;
	1487	uInt nrow = tab.nrow() ;
[1974]	1488	if ( nrow == 0 )
	1489	return 0 ;
[1987]	1490	ROScalarMeasColumn<MEpoch> tcol( tab, "TIME" ) ;
	1491	ROArrayMeasColumn<MDirection> dmcol( tab, "DIRECTION" ) ;
	1492	ROArrayColumn<Double> dcol( tab, "DIRECTION" ) ;
[2003]	1493	// binary search if idx == 0
	1494	if ( idx == 0 ) {
	1495	uInt nrowb = 75000 ;
[2004]	1496	if ( nrow > nrowb ) {
[2003]	1497	uInt nblock = nrow / nrowb + 1 ;
	1498	for ( uInt iblock = 0 ; iblock < nblock ; iblock++ ) {
[2004]	1499	uInt high = min( nrowb, nrow-iblock*nrowb ) ;
[2003]	1500
	1501	if ( tcol( high-1 ).get( "s" ).getValue() < t ) {
[2004]	1502	idx = iblock * nrowb ;
[2003]	1503	continue ;
	1504	}
	1505
	1506	Vector<MEpoch> tarr( high ) ;
	1507	for ( uInt irow = 0 ; irow < high ; irow++ ) {
[2004]	1508	tarr[irow] = tcol( iblock*nrowb+irow ) ;
[2003]	1509	}
	1510
	1511	uInt bidx = binarySearch( tarr, t ) ;
	1512
[2004]	1513	idx = iblock * nrowb + bidx ;
[2003]	1514	break ;
	1515	}
	1516	}
	1517	else {
	1518	Vector<MEpoch> tarr( nrow ) ;
	1519	for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
	1520	tarr[irow] = tcol( irow ) ;
	1521	}
	1522	idx = binarySearch( tarr, t ) ;
	1523	}
	1524	}
[1974]	1525	// ensure that tcol(idx) < t
	1526	//os_ << "tcol(idx) = " << tcol(idx).get("s").getValue() << " t = " << t << " diff = " << tcol(idx).get("s").getValue()-t << endl ;
	1527	while ( tcol(idx).get("s").getValue() > t && idx > 0 )
	1528	idx-- ;
	1529	//os_ << "idx = " << idx << LogIO::POST ;
	1530
	1531	// index search
	1532	for ( uInt i = idx ; i < nrow ; i++ ) {
	1533	Double tref = tcol( i ).get( "s" ).getValue() ;
	1534	if ( tref == t ) {
	1535	idx = i ;
	1536	break ;
	1537	}
	1538	else if ( tref > t ) {
	1539	if ( i == 0 ) {
	1540	idx = i ;
	1541	}
	1542	else {
	1543	idx = i-1 ;
	1544	doInterp = True ;
	1545	}
	1546	break ;
	1547	}
	1548	else {
	1549	idx = nrow - 1 ;
	1550	}
	1551	}
	1552	//os_ << "searched idx = " << idx << LogIO::POST ;
	1553
	1554	//os_ << "dmcol(idx).shape() = " << dmcol(idx).shape() << LogIO::POST ;
	1555	IPosition ip( dmcol(idx).shape().nelements(), 0 ) ;
	1556	//os_ << "ip = " << ip << LogIO::POST ;
	1557	ref = dmcol(idx)(ip).getRefString() ;
	1558	//os_ << "ref = " << ref << LogIO::POST ;
	1559	if ( doInterp ) {
	1560	//os_ << "do interpolation" << LogIO::POST ;
	1561	//os_ << "dcol(idx).shape() = " << dcol(idx).shape() << LogIO::POST ;
	1562	Double tref0 = tcol(idx).get("s").getValue() ;
	1563	Double tref1 = tcol(idx+1).get("s").getValue() ;
[1987]	1564	Matrix<Double> mdir0 = dcol( idx ) ;
	1565	Matrix<Double> mdir1 = dcol( idx+1 ) ;
	1566	Vector<Double> dir0 = mdir0.column( 0 ) ;
[1974]	1567	//os_ << "dir0 = " << dir0 << LogIO::POST ;
[1987]	1568	Vector<Double> dir1 = mdir1.column( 0 ) ;
[1974]	1569	//os_ << "dir1 = " << dir1 << LogIO::POST ;
	1570	Double dt0 = t - tref0 ;
	1571	Double dt1 = tref1 - t ;
	1572	dir.reference( (dt0dir1+dt1dir0)/(dt0+dt1) ) ;
[1987]	1573	if ( mdir0.ncolumn() > 1 ) {
	1574	if ( dt0 >= dt1 )
	1575	srate.reference( mdir0.column( 1 ) ) ;
	1576	else
	1577	srate.reference( mdir1.column( 1 ) ) ;
[1974]	1578	}
	1579	//os_ << "dir = " << dir << LogIO::POST ;
	1580	}
	1581	else {
	1582	//os_ << "no interpolation" << LogIO::POST ;
[1987]	1583	Matrix<Double> mdir0 = dcol( idx ) ;
	1584	dir.reference( mdir0.column( 0 ) ) ;
	1585	if ( mdir0.ncolumn() > 1 )
	1586	srate.reference( mdir0.column( 1 ) ) ;
[1974]	1587	}
	1588
[2004]	1589	// double endSec = gettimeofday_sec() ;
	1590	// os_ << "end MSFiller::getDirection() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
[1974]	1591	return idx ;
	1592	}
	1593
[1990]	1594	String MSFiller::keyTcal( Int feedid, Int spwid, String stime )
	1595	{
	1596	String sfeed = "FEED" + String::toString( feedid ) ;
	1597	String sspw = "SPW" + String::toString( spwid ) ;
	1598	return sfeed+":"+sspw+":"+stime ;
	1599	}
	1600
[2003]	1601	uInt MSFiller::binarySearch( Vector<MEpoch> &timeList, Double target )
	1602	{
	1603	Int low = 0 ;
	1604	Int high = timeList.nelements() ;
	1605	uInt idx = 0 ;
	1606
	1607	while ( low <= high ) {
	1608	idx = (Int)( 0.5 * ( low + high ) ) ;
	1609	Double t = timeList[idx].get( "s" ).getValue() ;
	1610	if ( t < target )
	1611	low = idx + 1 ;
	1612	else if ( t > target )
	1613	high = idx - 1 ;
	1614	else
	1615	return idx ;
	1616	}
	1617
	1618	idx = max( 0, min( low, high ) ) ;
	1619
	1620	return idx ;
	1621
	1622	}
	1623
[1974]	1624	} ;
	1625

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

Context Navigation

source: trunk/src/MSFiller.cpp@ 2023

Download in other formats: