Context Navigation

PKSMS2reader.cc @ 1880

Visit:

Last change on this file since 1880 was 1880, checked in by Takeshi Nakazato, 14 years ago

New Development: No

JIRA Issue: No

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

An implementation of asapmath.splitant() is changed since
the new filler doesn't support antenna parameter at the
moment.

The default behavior of MS filler also changed.

File size: 42.4 KB

Line
1	//#---------------------------------------------------------------------------
2	//# PKSMS2reader.cc: Class to read Parkes Multibeam data from a v2 MS.
3	//#---------------------------------------------------------------------------
4	//# livedata - processing pipeline for single-dish, multibeam spectral data.
5	//# Copyright (C) 2000-2009, Australia Telescope National Facility, CSIRO
6	//#
7	//# This file is part of livedata.
8	//#
9	//# livedata is free software: you can redistribute it and/or modify it under
10	//# the terms of the GNU General Public License as published by the Free
11	//# Software Foundation, either version 3 of the License, or (at your option)
12	//# any later version.
13	//#
14	//# livedata is distributed in the hope that it will be useful, but WITHOUT
15	//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16	//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
17	//# more details.
18	//#
19	//# You should have received a copy of the GNU General Public License along
20	//# with livedata. If not, see <http://www.gnu.org/licenses/>.
21	//#
22	//# Correspondence concerning livedata may be directed to:
23	//# Internet email: mcalabre@atnf.csiro.au
24	//# Postal address: Dr. Mark Calabretta
25	//# Australia Telescope National Facility, CSIRO
26	//# PO Box 76
27	//# Epping NSW 1710
28	//# AUSTRALIA
29	//#
30	//# http://www.atnf.csiro.au/computing/software/livedata.html
31	//# $Id: PKSMS2reader.cc,v 19.23 2009-09-29 07:33:38 cal103 Exp $
32	//#---------------------------------------------------------------------------
33	//# Original: 2000/08/03, Mark Calabretta, ATNF
34	//#---------------------------------------------------------------------------
35
36	// AIPS++ includes.
37	#include <casa/stdio.h>
38	#include <casa/Arrays/ArrayMath.h>
39	#include <casa/Arrays/Slice.h>
40	#include <ms/MeasurementSets/MSColumns.h>
41	#include <tables/Tables.h>
42	#include <casa/Quanta/MVTime.h>
43	#include <casa/Quanta/MVAngle.h>
44	#include <casa/BasicMath/Math.h>
45	#include <casa/Logging/LogIO.h>
46	#include <casa/Utilities/Sort.h>
47	#include <measures/Measures/MeasConvert.h>
48	#include <measures/Measures/MEpoch.h>
49	#include <measures/Measures/MeasRef.h>
50
51
52	// Parkes includes.
53	#include <atnf/pks/pks_maths.h>
54	#include <atnf/PKSIO/PKSrecord.h>
55	#include <atnf/PKSIO/PKSMS2reader.h>
56
57
58	//------------------------------------------------- PKSMS2reader::PKSMS2reader
59
60	// Default constructor.
61
62	PKSMS2reader::PKSMS2reader()
63	{
64	cMSopen = False;
65	}
66
67	//------------------------------------------------ PKSMS2reader::~PKSMS2reader
68
69	PKSMS2reader::~PKSMS2reader()
70	{
71	close();
72	}
73
74	//--------------------------------------------------------- PKSMS2reader::open
75
76	// Open the MS for reading.
77
78	Int PKSMS2reader::open(
79	const String msName,
80	const String antenna,
81	Vector<Bool> &beams,
82	Vector<Bool> &IFs,
83	Vector<uInt> &nChan,
84	Vector<uInt> &nPol,
85	Vector<Bool> &haveXPol,
86	Bool &haveBase,
87	Bool &haveSpectra)
88	{
89	// Check that MS is readable.
90	if (!MS::isReadable(msName)) {
91	return 1;
92	}
93
94	if (cMSopen) {
95	close();
96	}
97
98	cPKSMS = MeasurementSet(msName);
99
100	// data selection by antenna
101	if ( antenna.length() == 0 ) {
102	cAntId.resize( 1 ) ;
103	//cAntId[0] = 0 ;
104	ROScalarColumn<Int> ant1Col( cPKSMS, "ANTENNA1" ) ;
105	cAntId[0] = ant1Col(0) ;
106	}
107	else {
108	setupAntennaList( antenna ) ;
109	if ( cAntId.size() > 1 ) {
110	LogIO os( LogOrigin( "PKSMS2reader", "open()", WHERE ) ) ;
111	os << LogIO::WARN << "PKSMS2reader is not ready for multiple antenna selection. Use first antenna id " << cAntId[0] << "."<< LogIO::POST ;
112	Int tmp = cAntId[0] ;
113	cAntId.resize( 1 ) ;
114	cAntId[0] = tmp ;
115	}
116	stringstream ss ;
117	ss << "SELECT FROM $1 WHERE ANTENNA1 == ANTENNA2 && ANTENNA1 IN [" ;
118	for ( uInt i = 0 ; i < cAntId.size() ; i++ ) {
119	ss << cAntId[i] ;
120	if ( i == cAntId.size()-1 ) {
121	ss << "]" ;
122	}
123	else {
124	ss << "," ;
125	}
126	}
127	string taql = ss.str() ;
128	//cerr << "taql = " << taql << endl ;
129	cPKSMS = MeasurementSet( tableCommand( taql, cPKSMS ) ) ;
130	}
131
132	// taql access to the syscal table
133	cHaveSysCal = False;
134	if (cHaveSysCal=Table::isReadable(cPKSMS.sysCalTableName())) {
135	cSysCalTab = Table(cPKSMS.sysCalTableName());
136	}
137
138	// Lock the table for read access.
139	cPKSMS.lock(False);
140
141	cIdx = 0;
142	lastmjd = 0.0;
143	cNRow = cPKSMS.nrow();
144	cMSopen = True;
145
146	// Main MS table and subtable column access.
147	ROMSMainColumns msCols(cPKSMS);
148	ROMSDataDescColumns dataDescCols(cPKSMS.dataDescription());
149	ROMSFeedColumns feedCols(cPKSMS.feed());
150	ROMSFieldColumns fieldCols(cPKSMS.field());
151	ROMSPointingColumns pointingCols(cPKSMS.pointing());
152	ROMSPolarizationColumns polarizationCols(cPKSMS.polarization());
153	ROMSSourceColumns sourceCols(cPKSMS.source());
154	ROMSSpWindowColumns spWinCols(cPKSMS.spectralWindow());
155	ROMSStateColumns stateCols(cPKSMS.state());
156	ROMSSysCalColumns sysCalCols(cPKSMS.sysCal());
157	ROMSWeatherColumns weatherCols(cPKSMS.weather());
158	ROMSAntennaColumns antennaCols(cPKSMS.antenna());
159
160	// Column accessors for required columns.
161	cScanNoCol.reference(msCols.scanNumber());
162	cTimeCol.reference(msCols.time());
163	cIntervalCol.reference(msCols.interval());
164
165	cFieldIdCol.reference(msCols.fieldId());
166	cFieldNameCol.reference(fieldCols.name());
167	cFieldDelayDirCol.reference(fieldCols.delayDir());
168
169	cSrcIdCol.reference(fieldCols.sourceId());
170	cSrcId2Col.reference(sourceCols.sourceId());
171	cSrcNameCol.reference(sourceCols.name());
172	cSrcDirCol.reference(sourceCols.direction());
173	cSrcPMCol.reference(sourceCols.properMotion());
174	cSrcRestFrqCol.reference(sourceCols.restFrequency());
175
176	cStateIdCol.reference(msCols.stateId());
177	cObsModeCol.reference(stateCols.obsMode());
178	cCalCol.reference(stateCols.cal());
179	cSigStateCol.reference(stateCols.sig());
180	cRefStateCol.reference(stateCols.ref());
181
182	cDataDescIdCol.reference(msCols.dataDescId());
183	cSpWinIdCol.reference(dataDescCols.spectralWindowId());
184	cChanFreqCol.reference(spWinCols.chanFreq());
185	cTotBWCol.reference(spWinCols.totalBandwidth());
186
187	cWeatherTimeCol.reference(weatherCols.time());
188	cTemperatureCol.reference(weatherCols.temperature());
189	cPressureCol.reference(weatherCols.pressure());
190	cHumidityCol.reference(weatherCols.relHumidity());
191
192	cBeamNoCol.reference(msCols.feed1());
193	cPointingCol.reference(pointingCols.direction());
194	cPointingTimeCol.reference(pointingCols.time());
195	cSigmaCol.reference(msCols.sigma());
196	cNumReceptorCol.reference(feedCols.numReceptors());
197
198	// Optional columns.
199	cHaveTsys = False;
200	cHaveTcal = False;
201	if ((cHaveSrcVel = cPKSMS.source().tableDesc().isColumn("SYSVEL"))) {
202	cSrcVelCol.attach(cPKSMS.source(), "SYSVEL");
203	}
204
205	if (cHaveSysCal && (cHaveTsys = cPKSMS.sysCal().tableDesc().isColumn("TSYS"))) {
206	cTsysCol.attach(cPKSMS.sysCal(), "TSYS");
207	}
208
209	if (cHaveSysCal && (cHaveTcal = cPKSMS.sysCal().tableDesc().isColumn("TCAL"))) {
210	cTcalCol.attach(cPKSMS.sysCal(), "TCAL");
211	}
212
213	if ((cHaveCalFctr = cPKSMS.tableDesc().isColumn("CALFCTR"))) {
214	cCalFctrCol.attach(cPKSMS, "CALFCTR");
215	}
216
217	if ((cHaveBaseLin = cPKSMS.tableDesc().isColumn("BASELIN"))) {
218	cBaseLinCol.attach(cPKSMS, "BASELIN");
219	cBaseSubCol.attach(cPKSMS, "BASESUB");
220	}
221
222	// Spectral data should always be present.
223	haveSpectra = True;
224	cHaveDataCol = False;
225	cHaveCorrectedDataCol = False;
226	ROMSObservationColumns observationCols(cPKSMS.observation());
227	//String telName = observationCols.telescopeName()(0);
228	cTelName = observationCols.telescopeName()(0);
229	//cATF = cTelName.contains("ATF");
230	//cOSF = cTelName.contains("OSF");
231	//cALMA = cTelName.contains("ALMA");
232	cALMA = cTelName.contains("ATF")\|\|cTelName.contains("OSF")\|\|
233	cTelName.contains("ALMA");
234
235	if (cHaveDataCol = cPKSMS.isColumn(MSMainEnums::DATA)) {
236	if (cALMA) {
237	//try to read a single baseline interferometeric data
238	//and treat it as single dish data
239	//maybe extended for ALMA commissioning later
240	cDataCol.reference(msCols.data());
241	if (cHaveCorrectedDataCol = cPKSMS.isColumn(MSMainEnums::CORRECTED_DATA)) {
242	//cerr<<"Do have CORRECTED_DATA column"<<endl;
243	cCorrectedDataCol.reference(msCols.correctedData());
244	}
245	}
246	}
247	else {
248	cFloatDataCol.reference(msCols.floatData());
249	}
250	cFlagCol.reference(msCols.flag());
251	cFlagRowCol.reference(msCols.flagRow());
252
253	if (cGetXPol = (cPKSMS.isColumn(MSMainEnums::DATA) && (!cALMA))) {
254	if ((cHaveXCalFctr = cPKSMS.tableDesc().isColumn("XCALFCTR"))) {
255	cXCalFctrCol.attach(cPKSMS, "XCALFCTR");
256	}
257
258	cDataCol.reference(msCols.data());
259	}
260
261	// Find which beams are present in the data.
262	Vector<Int> beamNos = cBeamNoCol.getColumn();
263	Int maxBeamNo = max(beamNos) + 1;
264	beams.resize(maxBeamNo);
265
266	beams = False;
267	for (uInt irow = 0; irow < beamNos.nelements(); irow++) {
268	beams(beamNos(irow)) = True;
269	}
270
271
272	// Number of IFs.
273	//uInt nIF = dataDescCols.nrow();
274	uInt nIF =spWinCols.nrow();
275	Vector<Int> spWinIds = cSpWinIdCol.getColumn() ;
276	IFs.resize(nIF);
277	IFs = True;
278	for ( Int ispw = 0 ; ispw < nIF ; ispw++ ) {
279	if ( allNE( ispw, spWinIds ) ) {
280	IFs(ispw) = False ;
281	}
282	}
283
284	// Number of polarizations and channels in each IF.
285	ROScalarColumn<Int> numChanCol(spWinCols.numChan());
286
287	ROScalarColumn<Int> polIdCol(dataDescCols.polarizationId());
288	ROScalarColumn<Int> numPolCol(polarizationCols.numCorr());
289
290	nChan.resize(nIF);
291	nPol.resize(nIF);
292	for (uInt iIF = 0; iIF < nIF; iIF++) {
293	if ( IFs(iIF) ) {
294	nChan(iIF) = numChanCol(cSpWinIdCol(iIF)) ;
295	nPol(iIF) = numPolCol(polIdCol(iIF)) ;
296	}
297	else {
298	nChan(iIF) = 0 ;
299	nPol(iIF) = 0 ;
300	}
301	}
302
303	// Cross-polarization data present?
304	haveXPol.resize(nIF);
305	haveXPol = False;
306
307	if (cGetXPol && !(cALMA)) {
308	for (Int irow = 0; irow < cNRow; irow++) {
309	if (cDataCol.isDefined(irow)) {
310	Int iIF = cDataDescIdCol(irow);
311	haveXPol(iIF) = True;
312	}
313	}
314	}
315
316
317	// Initialize member data.
318	cBeams.assign(beams);
319	cIFs.assign(IFs);
320	cNChan.assign(nChan);
321	cNPol.assign(nPol);
322	cHaveXPol.assign(haveXPol);
323
324
325	// Default channel range selection.
326	cStartChan.resize(nIF);
327	cEndChan.resize(nIF);
328	cRefChan.resize(nIF);
329
330	for (uInt iIF = 0; iIF < nIF; iIF++) {
331	cStartChan(iIF) = 1;
332	cEndChan(iIF) = cNChan(iIF);
333	cRefChan(iIF) = cNChan(iIF)/2 + 1;
334	}
335
336	Slice all;
337	cDataSel.resize(nIF);
338	cDataSel = Slicer(all, all);
339
340	cScanNo = 0;
341	cCycleNo = 1;
342	cTime = cTimeCol(0);
343
344	return 0;
345	}
346
347	//---------------------------------------------------- PKSMS2reader::getHeader
348
349	// Get parameters describing the data.
350
351	Int PKSMS2reader::getHeader(
352	String &observer,
353	String &project,
354	String &antName,
355	Vector<Double> &antPosition,
356	// before merge...
357	//String &obsMode,
358	String &obsType,
359	String &bunit,
360	Float &equinox,
361	String &dopplerFrame,
362	Double &mjd,
363	Double &refFreq,
364	Double &bandwidth)
365	{
366	if (!cMSopen) {
367	return 1;
368	}
369
370	// Observer and project.
371	ROMSObservationColumns observationCols(cPKSMS.observation());
372	observer = observationCols.observer()(0);
373	project = observationCols.project()(0);
374
375	// Antenna name and ITRF coordinates.
376	ROMSAntennaColumns antennaCols(cPKSMS.antenna());
377	//antName = antennaCols.name()(0);
378	antName = antennaCols.name()(cAntId[0]);
379	if (cALMA) {
380	antName = cTelName + "-" + antName;
381	}
382	//antPosition = antennaCols.position()(0);
383	antPosition = antennaCols.position()(cAntId[0]);
384
385	// Observation type.
386	if (cObsModeCol.nrow()) {
387	obsType = cObsModeCol(0);
388	if (obsType == "\0") obsType = "RF";
389	} else {
390	obsType = "RF";
391	}
392
393	bunit = "";
394	if (cHaveDataCol) {
395	const TableRecord& keywordSet2
396	= cDataCol.columnDesc().keywordSet();
397	if(keywordSet2.isDefined("UNIT")) {
398	bunit = keywordSet2.asString("UNIT");
399	}
400	} else {
401	const TableRecord& keywordSet
402	= cFloatDataCol.columnDesc().keywordSet();
403	if(keywordSet.isDefined("UNIT")) {
404	bunit = keywordSet.asString("UNIT");
405	}
406	}
407
408	/***
409	const TableRecord& keywordSet
410	= cFloatDataCol.columnDesc().keywordSet();
411	if(keywordSet.isDefined("UNIT")) {
412	fluxunit = keywordSet.asString("UNIT");
413	}
414	***/
415	// Coordinate equinox.
416	ROMSPointingColumns pointingCols(cPKSMS.pointing());
417	String dirref = pointingCols.direction().keywordSet().asRecord("MEASINFO").
418	asString("Ref");
419	cDirRef = dirref;
420	if (dirref =="AZELGEO" \|\| dirref == "AZEL") {
421	dirref = "J2000";
422	}
423	sscanf(dirref.chars()+1, "%f", &equinox);
424
425	// Frequency/velocity reference frame.
426	ROMSSpWindowColumns spWinCols(cPKSMS.spectralWindow());
427	dopplerFrame = MFrequency::showType(spWinCols.measFreqRef()(0));
428
429	// Translate to FITS standard names.
430	if (dopplerFrame == "TOPO") {
431	dopplerFrame = "TOPOCENT";
432	} else if (dopplerFrame == "GEO") {
433	dopplerFrame = "GEOCENTR";
434	} else if (dopplerFrame == "BARY") {
435	dopplerFrame = "BARYCENT";
436	} else if (dopplerFrame == "GALACTO") {
437	dopplerFrame = "GALACTOC";
438	} else if (dopplerFrame == "LGROUP") {
439	dopplerFrame = "LOCALGRP";
440	} else if (dopplerFrame == "CMB") {
441	dopplerFrame = "CMBDIPOL";
442	} else if (dopplerFrame == "REST") {
443	dopplerFrame = "SOURCE";
444	}
445
446	// MJD at start of observation.
447	mjd = cTimeCol(0)/86400.0;
448
449	// Reference frequency and bandwidth.
450	refFreq = spWinCols.refFrequency()(0);
451	bandwidth = spWinCols.totalBandwidth()(0);
452
453	return 0;
454	}
455
456	//-------------------------------------------------- PKSMS2reader::getFreqInfo
457
458	// Get frequency parameters for each IF.
459
460	Int PKSMS2reader::getFreqInfo(
461	Vector<Double> &startFreq,
462	Vector<Double> &endFreq)
463	{
464	uInt nIF = cIFs.nelements();
465	startFreq.resize(nIF);
466	endFreq.resize(nIF);
467
468	for (uInt iIF = 0; iIF < nIF; iIF++) {
469	Vector<Double> chanFreq = cChanFreqCol(iIF);
470
471	Int nChan = chanFreq.nelements();
472	startFreq(iIF) = chanFreq(0);
473	endFreq(iIF) = chanFreq(nChan-1);
474	}
475
476	return 0;
477	}
478
479	//------------------------------------------------------- PKSMS2reader::select
480
481	// Set data selection by beam number and channel.
482
483	uInt PKSMS2reader::select(
484	const Vector<Bool> beamSel,
485	const Vector<Bool> IFsel,
486	const Vector<Int> startChan,
487	const Vector<Int> endChan,
488	const Vector<Int> refChan,
489	const Bool getSpectra,
490	const Bool getXPol,
491	const Bool getFeedPos,
492	const Bool getPointing,
493	const Int coordSys)
494	{
495	if (!cMSopen) {
496	return 1;
497	}
498
499	// Beam selection.
500	uInt nBeam = cBeams.nelements();
501	uInt nBeamSel = beamSel.nelements();
502	for (uInt ibeam = 0; ibeam < nBeam; ibeam++) {
503	if (ibeam < nBeamSel) {
504	cBeams(ibeam) = beamSel(ibeam);
505	} else {
506	cBeams(ibeam) = False;
507	}
508	}
509
510	uInt nIF = cIFs.nelements();
511	uInt maxNChan = 0;
512	for (uInt iIF = 0; iIF < nIF; iIF++) {
513	// IF selection.
514	if (iIF < IFsel.nelements()) {
515	cIFs(iIF) = IFsel(iIF);
516	} else {
517	cIFs(iIF) = False;
518	}
519
520	if (!cIFs(iIF)) continue;
521
522
523	// Channel selection.
524	if (iIF < startChan.nelements()) {
525	cStartChan(iIF) = startChan(iIF);
526
527	if (cStartChan(iIF) <= 0) {
528	cStartChan(iIF) += cNChan(iIF);
529	} else if (cStartChan(iIF) > Int(cNChan(iIF))) {
530	cStartChan(iIF) = cNChan(iIF);
531	}
532	}
533
534	if (iIF < endChan.nelements()) {
535	cEndChan(iIF) = endChan(iIF);
536
537	if (cEndChan(iIF) <= 0) {
538	cEndChan(iIF) += cNChan(iIF);
539	} else if (cEndChan(iIF) > Int(cNChan(iIF))) {
540	cEndChan(iIF) = cNChan(iIF);
541	}
542	}
543
544	if (iIF < refChan.nelements()) {
545	cRefChan(iIF) = refChan(iIF);
546	} else {
547	cRefChan(iIF) = cStartChan(iIF);
548	if (cStartChan(iIF) <= cEndChan(iIF)) {
549	cRefChan(iIF) += (cEndChan(iIF) - cStartChan(iIF) + 1)/2;
550	} else {
551	cRefChan(iIF) -= (cStartChan(iIF) - cEndChan(iIF) + 1)/2;
552	}
553	}
554
555	uInt nChan = abs(cEndChan(iIF) - cStartChan(iIF)) + 1;
556	if (maxNChan < nChan) {
557	maxNChan = nChan;
558	}
559
560	// Inverted Slices are not allowed.
561	Slice outPols;
562	Slice outChans(min(cStartChan(iIF),cEndChan(iIF))-1, nChan);
563	cDataSel(iIF) = Slicer(outPols, outChans);
564	}
565
566	// Get spectral data?
567	cGetSpectra = getSpectra;
568
569	// Get cross-polarization data?
570	cGetXPol = cGetXPol && getXPol;
571
572	// Get feed positions? (Not available.)
573	cGetFeedPos = False;
574
575	// Get Pointing data (for MS)
576	cGetPointing = getPointing;
577
578	// Coordinate system? (Only equatorial available.)
579	cCoordSys = 0;
580
581	return maxNChan;
582	}
583
584	//---------------------------------------------------- PKSMS2reader::findRange
585
586	// Find the range of the data in time and position.
587
588	Int PKSMS2reader::findRange(
589	Int &nRow,
590	Int &nSel,
591	Vector<Double> &timeSpan,
592	Matrix<Double> &positions)
593	{
594	if (!cMSopen) {
595	return 1;
596	}
597
598	nRow = cNRow;
599
600	// Find the number of rows selected.
601	nSel = 0;
602	Vector<Bool> sel(nRow);
603	for (Int irow = 0; irow < nRow; irow++) {
604	if ((sel(irow) = cBeams(cBeamNoCol(irow)) &&
605	cIFs(cDataDescIdCol(irow)))) {
606	nSel++;
607	}
608	}
609
610	// Find the time range (s).
611	timeSpan.resize(2);
612	timeSpan(0) = cTimeCol(0);
613	timeSpan(1) = cTimeCol(nRow-1);
614
615	// Retrieve positions for selected data.
616	Int isel = 0;
617	positions.resize(2,nSel);
618	for (Int irow = 0; irow < nRow; irow++) {
619	if (sel(irow)) {
620	Matrix<Double> pointingDir = cPointingCol(cFieldIdCol(irow));
621	positions.column(isel++) = pointingDir.column(0);
622	}
623	}
624
625	return 0;
626	}
627
628	//--------------------------------------------------------- PKSMS2reader::read
629
630	// Read the next data record.
631
632	/**
633	Int PKSMS2reader::read(
634	Int &scanNo,
635	Int &cycleNo,
636	Double &mjd,
637	Double &interval,
638	String &fieldName,
639	String &srcName,
640	Vector<Double> &srcDir,
641	Vector<Double> &srcPM,
642	Double &srcVel,
643	String &obsMode,
644	Int &IFno,
645	Double &refFreq,
646	Double &bandwidth,
647	Double &freqInc,
648	Vector<Double> &restFreq,
649	Vector<Float> &tcal,
650	String &tcalTime,
651	Float &azimuth,
652	Float &elevation,
653	Float &parAngle,
654	Float &focusAxi,
655	Float &focusTan,
656	Float &focusRot,
657	Float &temperature,
658	Float &pressure,
659	Float &humidity,
660	Float &windSpeed,
661	Float &windAz,
662	Int &refBeam,
663	Int &beamNo,
664	Vector<Double> &direction,
665	Vector<Double> &scanRate,
666	Vector<Float> &tsys,
667	Vector<Float> &sigma,
668	Vector<Float> &calFctr,
669	Matrix<Float> &baseLin,
670	Matrix<Float> &baseSub,
671	Matrix<Float> &spectra,
672	Matrix<uChar> &flagged,
673	uInt &flagrow,
674	Complex &xCalFctr,
675	Vector<Complex> &xPol)
676	**/
677	Int PKSMS2reader::read(PKSrecord &pksrec)
678	{
679	LogIO os( LogOrigin( "PKSMS2reader", "read()", WHERE ) ) ;
680
681	if (!cMSopen) {
682	return 1;
683	}
684
685	// Check for EOF.
686	if (cIdx >= cNRow) {
687	return -1;
688	}
689
690	// Find the next selected beam and IF.
691	Int ibeam;
692	Int iIF;
693	Int iDataDesc;
694
695	while (True) {
696	ibeam = cBeamNoCol(cIdx);
697	iDataDesc = cDataDescIdCol(cIdx);
698	iIF =cSpWinIdCol(iDataDesc);
699	if (cBeams(ibeam) && cIFs(iIF)) {
700	break;
701	}
702
703	// Check for EOF.
704	if (++cIdx >= cNRow) {
705	return -1;
706	}
707	}
708	// Renumerate scan no. Here still is 1-based
709	//scanNo = cScanNoCol(cIdx) - cScanNoCol(0) + 1;
710	//scanNo = cScanNoCol(cIdx);
711	pksrec.scanNo = cScanNoCol(cIdx);
712
713	if (pksrec.scanNo != cScanNo) {
714	// Start of new scan.
715	cScanNo = pksrec.scanNo;
716	cCycleNo = 1;
717	cTime = cTimeCol(cIdx);
718	}
719
720	Double time = cTimeCol(cIdx);
721	pksrec.mjd = time/86400.0;
722	pksrec.interval = cIntervalCol(cIdx);
723
724	// Reconstruct the integration cycle number; due to small latencies the
725	// integration time is usually slightly less than the time between cycles,
726	// resetting cTime will prevent the difference from accumulating.
727	cCycleNo += nint((time - cTime)/pksrec.interval);
728	pksrec.cycleNo = cCycleNo;
729	cTime = time;
730
731	Int fieldId = cFieldIdCol(cIdx);
732	pksrec.fieldName = cFieldNameCol(fieldId);
733
734	Int srcId = cSrcIdCol(fieldId);
735	//For source with multiple spectral window setting, this is not
736	// correct. Source name of srcId may not be at 'srcId'th row of SrcNameCol
737	//srcName = cSrcNameCol(srcId);
738	for (uInt irow = 0; irow < cSrcId2Col.nrow(); irow++) {
739	if (cSrcId2Col(irow) == srcId) {
740	//srcName = cSrcNameCol(irow);
741	pksrec.srcName = cSrcNameCol(irow);
742	}
743	}
744
745	pksrec.srcDir = cSrcDirCol(srcId);
746	pksrec.srcPM = cSrcPMCol(srcId);
747
748	// Systemic velocity.
749	if (!cHaveSrcVel \|\| cALMA) {
750	pksrec.srcVel = 0.0f;
751	} else {
752	pksrec.srcVel = cSrcVelCol(srcId)(IPosition(1,0));
753	}
754
755	ROMSAntennaColumns antennaCols(cPKSMS.antenna());
756	//String telescope = antennaCols.name()(0);
757	String telescope = antennaCols.name()(cAntId[0]);
758	Bool cGBT = telescope.contains("GBT");
759	//Bool cPM = telescope.contains("PM"); // ACA TP antenna
760	//Bool cDV = telescope.contains("DV"); // VERTEX
761	//Bool cCM = telescope.contains("CM"); // ACA 7m antenna
762	//Bool cALMA = cPM \|\| cDV \|\| cCM ;
763	// Observation type.
764	// check if State Table exist
765	//Bool cHaveStateTab=Table::isReadable(cPKSMS.stateTableName());
766	Int stateId = 0;
767	Int StateNRow = 0;
768	StateNRow=cObsModeCol.nrow();
769	if (Table::isReadable(cPKSMS.stateTableName())) {
770	pksrec.obsType = " ";
771	if (StateNRow > 0) {
772	stateId = cStateIdCol(cIdx);
773	if (stateId == -1) {
774	//pksrec.obsType = " ";
775	} else {
776	pksrec.obsType = cObsModeCol(stateId);
777	Bool sigState =cSigStateCol(stateId);
778	Bool refState =cRefStateCol(stateId);
779	//DEBUG
780	//cerr <<"stateid="<<stateId<<" obsmode="<<pksrec.obsType<<endl;
781	if (cGBT) {
782	// split the obsType string and append a proper label
783	// (these are GBT specific)
784	int epos = pksrec.obsType.find_first_of(':');
785	int nextpos = pksrec.obsType.find_first_of(':',epos+1);
786	string obsMode1 = pksrec.obsType.substr(0,epos);
787	string obsMode2 = pksrec.obsType.substr(epos+1,nextpos-epos-1);
788
789	//cerr <<"obsMode2= "<<obsMode2<<endl;
790	if (!pksrec.srcName.contains("_ps")
791	&&!pksrec.srcName.contains("_psr")
792	&&!pksrec.srcName.contains("_nod")
793	&&!pksrec.srcName.contains("_fs")
794	&&!pksrec.srcName.contains("_fsr")) {
795	// if Nod mode observation , append '_nod'
796	if (obsMode1 == "Nod") {
797	//pksrec.srcName.append("_nod");
798	pksrec.srcType = SrcType::NOD ;
799	} else if (obsMode1 == "OffOn") {
800	// for GBT position switch observations (OffOn or OnOff)
801	//if (obsMode2 == "PSWITCHON") pksrec.srcName.append("_ps");
802	//if (obsMode2 == "PSWITCHOFF") pksrec.srcName.append("_psr");
803	if (obsMode2 == "PSWITCHON") pksrec.srcType = SrcType::PSON ;
804	if (obsMode2 == "PSWITCHOFF") pksrec.srcType = SrcType::PSOFF ;
805	} else {
806	if (obsMode2 == "FSWITCH") {
807	// for GBT frequency switch mode
808	//if (sigState) pksrec.srcName.append("_fs");
809	//if (refState) pksrec.srcName.append("_fsr");
810	if (sigState) pksrec.srcType = SrcType::FSON ;
811	if (refState) pksrec.srcType = SrcType::FSOFF ;
812	}
813	}
814	}
815	}
816	else if (cALMA) {
817	// ALMA tag
818	// split the obsType string and append a proper label
819	string substr[1] ;
820	int numSubstr = split( pksrec.obsType, substr, 1, "," );
821	String obsType = String( substr[0] );
822	int epos = obsType.find_first_of('.');
823	int nextpos = obsType.find_first_of('.',epos+1);
824	string obsMode1 = obsType.substr(0,epos);
825	string obsMode2 = obsType.substr(epos+1,nextpos-epos-1);
826
827	//cerr <<"obsMode2= "<<obsMode2<<endl;
828	// Current OBS_MODE format:
829	//
830	// ON: OBSERVE_TARGET.ON_SOURCE
831	// OFF: OBSERVE_TARGET.OFF_SOURCE
832	//
833	if (obsMode1 == "OBSERVE_TARGET") {
834	//if (obsMode2 == "ON_SOURCE") pksrec.srcName.append("_pson");
835	//if (obsMode2 == "OFF_SOURCE") pksrec.srcName.append("_psoff");
836	if (obsMode2 == "ON_SOURCE") pksrec.srcType = SrcType::PSON ;
837	if (obsMode2 == "OFF_SOURCE") pksrec.srcType = SrcType::PSOFF ;
838	}
839	}
840	}
841	}
842	}
843	// CAL state
844	// this should be apply just for GBT data?
845	Double Cal;
846	if (stateId==-1 \|\| StateNRow==0) {
847	Cal = 0;
848	} else {
849	Cal = cCalCol(stateId);
850	}
851	if (cGBT) {
852	if (Cal > 0 && !pksrec.srcName.contains("_calon")) {
853	//pksrec.srcName.append("_calon");
854	if ( pksrec.srcType == SrcType::NOD )
855	pksrec.srcType = SrcType::NODCAL ;
856	else if ( pksrec.srcType == SrcType::PSON )
857	pksrec.srcType = SrcType::PONCAL ;
858	else if ( pksrec.srcType == SrcType::PSOFF )
859	pksrec.srcType = SrcType::POFFCAL ;
860	else if ( pksrec.srcType == SrcType::FSON )
861	pksrec.srcType = SrcType::FONCAL ;
862	else if ( pksrec.srcType == SrcType::FSOFF )
863	pksrec.srcType = SrcType::FOFFCAL ;
864	else
865	pksrec.srcName.append("_calon");
866	}
867	}
868
869	pksrec.IFno = iIF + 1;
870	Int nChan = abs(cEndChan(iIF) - cStartChan(iIF)) + 1;
871
872	// Minimal handling on continuum data.
873	Vector<Double> chanFreq = cChanFreqCol(iIF);
874	if (nChan == 1) {
875	//pksrec.freqInc = chanFreq(0);
876	pksrec.freqInc = cTotBWCol(iIF);
877	pksrec.refFreq = chanFreq(0);
878	pksrec.restFreq.resize(1);
879	pksrec.restFreq[0] = 0.0f;
880	} else {
881
882	if (cStartChan(iIF) <= cEndChan(iIF)) {
883	pksrec.freqInc = chanFreq(1) - chanFreq(0);
884	} else {
885	pksrec.freqInc = chanFreq(0) - chanFreq(1);
886	}
887
888	pksrec.refFreq = chanFreq(cRefChan(iIF)-1);
889
890	Bool HaveSrcRestFreq= cSrcRestFrqCol.isDefined(srcId);
891	if (HaveSrcRestFreq) {
892	//restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0));
893	//restFreq = cSrcRestFrqCol(srcId);
894	pksrec.restFreq = cSrcRestFrqCol(srcId);
895	} else {
896	pksrec.restFreq.resize(1);
897	pksrec.restFreq[0] = 0.0f;
898	}
899	}
900	//pksrec.bandwidth = abs(pksrec.freqInc * nChan);
901	pksrec.bandwidth = abs(cTotBWCol(0));
902
903	pksrec.tcal.resize(cNPol(iIF));
904	pksrec.tcal = 0.0f;
905	pksrec.tcalTime = "";
906	// pksrec.azimuth = 0.0f;
907	// pksrec.elevation = 0.0f;
908	pksrec.parAngle = 0.0f;
909
910	pksrec.focusAxi = 0.0f;
911	pksrec.focusTan = 0.0f;
912	pksrec.focusRot = 0.0f;
913
914	// Find the appropriate entry in the WEATHER subtable.
915	//Bool cHaveStateTab=Table::isReadable(cPKSMS.stateTableName());
916	Bool cHaveWeatherTab = Table::isReadable(cPKSMS.weatherTableName());
917	Int weatherIdx=-1;
918	if (cHaveWeatherTab) {
919	Vector<Double> wTimes = cWeatherTimeCol.getColumn();
920	for (weatherIdx = wTimes.nelements()-1; weatherIdx >= 0; weatherIdx--) {
921	if (cWeatherTimeCol(weatherIdx) <= time) {
922	break;
923	}
924	}
925	}
926
927	if (weatherIdx < 0 \|\| !cHaveWeatherTab) {
928	// No appropriate WEATHER entry.
929	pksrec.temperature = 0.0f;
930	pksrec.pressure = 0.0f;
931	pksrec.humidity = 0.0f;
932	} else {
933	pksrec.temperature = cTemperatureCol(weatherIdx);
934	pksrec.pressure = cPressureCol(weatherIdx);
935	pksrec.humidity = cHumidityCol(weatherIdx);
936	}
937
938	pksrec.windSpeed = 0.0f;
939	pksrec.windAz = 0.0f;
940
941	pksrec.refBeam = 0;
942	pksrec.beamNo = ibeam + 1;
943
944	//pointing/azel
945	MVPosition mvpos(antennaCols.position()(cAntId[0]));
946	MPosition mp(mvpos);
947	Quantum<Double> qt(time,"s");
948	MVEpoch mvt(qt);
949	MEpoch me(mvt);
950	MeasFrame frame(mp, me);
951	MDirection md;
952	pksrec.pCode = 0;
953	pksrec.rateAge = 0.0f;
954	pksrec.paRate = 0.0f;
955	if (cGetPointing) {
956	//cerr << "get pointing data ...." << endl;
957	ROScalarColumn<Int> pAntIdCol ;
958	ROScalarColumn<Double> psTimeCol ;
959	Table ptTable = cPKSMS.pointing() ;
960	MSPointing selPtTab( ptTable( ptTable.col("ANTENNA_ID") == cAntId[0] ) ) ;
961	pAntIdCol.attach( selPtTab, "ANTENNA_ID" ) ;
962	Vector<Int> antIds = pAntIdCol.getColumn() ;
963	psTimeCol.attach( selPtTab, "TIME" ) ;
964	Vector<Double> pTimes = psTimeCol.getColumn();
965	Bool doInterp = False ;
966	Int PtIdx=-1;
967	for (PtIdx = pTimes.nelements()-1; PtIdx >= 0; PtIdx--) {
968	if ( pTimes[PtIdx] == time ) {
969	break ;
970	}
971	else if ( pTimes[PtIdx] < time ) {
972	if ( PtIdx != pTimes.nelements()-1 ) {
973	doInterp = True ;
974	}
975	break ;
976	}
977	}
978	if ( PtIdx == -1 ) {
979	PtIdx = 0 ;
980	}
981	//cerr << "got index=" << PtIdx << endl;
982	Matrix<Double> pointingDir = cPointingCol(PtIdx);
983	ROMSPointingColumns PtCols( selPtTab ) ;
984	Vector<Double> pointingDirVec ;
985	if ( doInterp ) {
986	Double dt1 = time - pTimes[PtIdx] ;
987	Double dt2 = pTimes[PtIdx+1] - time ;
988	Vector<Double> dirVec1 = pointingDir.column(0) ;
989	Matrix<Double> pointingDir2 = cPointingCol(PtIdx+1) ;
990	Vector<Double> dirVec2 = pointingDir2.column(0) ;
991	pointingDirVec = (dt1dirVec2+dt2dirVec1)/(dt1+dt2) ;
992	Vector<MDirection> vmd1(1) ;
993	Vector<MDirection> vmd2(1) ;
994	PtCols.directionMeasCol().get(PtIdx,vmd1) ;
995	Vector<Double> angle1 = vmd1(0).getAngle().getValue("rad") ;
996	PtCols.directionMeasCol().get(PtIdx+1,vmd2) ;
997	Vector<Double> angle2 = vmd2(0).getAngle().getValue("rad") ;
998	Vector<Double> angle = (dt1angle2+dt2angle1)/(dt1+dt2) ;
999	Quantum< Vector<Double> > qangle( angle, "rad" ) ;
1000	String typeStr = vmd1(0).getRefString() ;
1001	//cerr << "vmd1.getRefString()=" << typeStr << endl ;
1002	MDirection::Types mdType ;
1003	MDirection::getType( mdType, typeStr ) ;
1004	//cerr << "mdType=" << mdType << endl ;
1005	md = MDirection( qangle, mdType ) ;
1006	//cerr << "md=" << md.getAngle().getValue("rad") << endl ;
1007	}
1008	else {
1009	pointingDirVec = pointingDir.column(0) ;
1010	Vector<MDirection> vmd(1);
1011	PtCols.directionMeasCol().get(PtIdx,vmd);
1012	md = vmd[0];
1013	}
1014	// put J2000 coordinates in "direction"
1015	if (cDirRef =="J2000") {
1016	pksrec.direction = pointingDirVec ;
1017	}
1018	else {
1019	pksrec.direction =
1020	MDirection::Convert(md, MDirection::Ref(MDirection::J2000,
1021	frame)
1022	)().getAngle("rad").getValue();
1023
1024	}
1025	uInt ncols = pointingDir.ncolumn();
1026	pksrec.scanRate.resize(2);
1027	if (ncols == 1) {
1028	pksrec.scanRate = 0.0f;
1029	} else {
1030	pksrec.scanRate(0) = pointingDir.column(1)(0);
1031	pksrec.scanRate(1) = pointingDir.column(1)(1);
1032	}
1033	}
1034	else {
1035	// Get direction from FIELD table
1036	// here, assume direction to be the field direction not pointing
1037	Matrix<Double> delayDir = cFieldDelayDirCol(fieldId);
1038	pksrec.direction = delayDir.column(0);
1039	uInt ncols = delayDir.ncolumn();
1040	pksrec.scanRate.resize(2);
1041	if (ncols == 1) {
1042	pksrec.scanRate = 0.0f;
1043	} else {
1044	pksrec.scanRate(0) = delayDir.column(1)(0);
1045	pksrec.scanRate(1) = delayDir.column(1)(1);
1046	}
1047	}
1048	// caluculate azimuth and elevation
1049	// first, get the reference frame
1050	/**
1051	MVPosition mvpos(antennaCols.position()(0));
1052	MPosition mp(mvpos);
1053	Quantum<Double> qt(time,"s");
1054	MVEpoch mvt(qt);
1055	MEpoch me(mvt);
1056	MeasFrame frame(mp, me);
1057	**/
1058	//
1059	ROMSFieldColumns fldCols(cPKSMS.field());
1060	Vector<MDirection> vmd(1);
1061	//MDirection md;
1062	fldCols.delayDirMeasCol().get(fieldId,vmd);
1063	md = vmd[0];
1064	//Vector<Double> dircheck = md.getAngle("rad").getValue();
1065	//cerr<<"dircheck="<<dircheck<<endl;
1066
1067	Vector<Double> azel =
1068	MDirection::Convert(md, MDirection::Ref(MDirection::AZEL,
1069	frame)
1070	)().getAngle("rad").getValue();
1071	//cerr<<"azel="<<azel<<endl;
1072	pksrec.azimuth = azel[0];
1073	pksrec.elevation = azel[1];
1074
1075	// Get Tsys assuming that entries in the SYSCAL table match the main table.
1076	if (cHaveTsys) {
1077	Int nTsysColRow = cTsysCol.nrow();
1078	if (nTsysColRow != cNRow) {
1079	cHaveTsys=0;
1080	}
1081	}
1082	if (cHaveTsys) {
1083	cTsysCol.get(cIdx, pksrec.tsys, True);
1084	} else {
1085	Int numReceptor;
1086	cNumReceptorCol.get(0, numReceptor);
1087	pksrec.tsys.resize(numReceptor);
1088	pksrec.tsys = 1.0f;
1089	}
1090	cSigmaCol.get(cIdx, pksrec.sigma, True);
1091
1092	//get Tcal if available
1093	if (cHaveTcal) {
1094	Int nTcalColRow = cTcalCol.nrow();
1095	uInt nBeam = cBeams.nelements();
1096	uInt nIF = cIFs.nelements();
1097	uInt nrws = nBeam * nIF;
1098	if (nTcalColRow > 0) {
1099	// find tcal match with the data with the data time stamp
1100	Double mjds = pksrec.mjd(243600);
1101	Double dtcalTime;
1102	if ( pksrec.mjd > lastmjd \|\| cIdx==0 ) {
1103	//Table tmptab = cSysCalTab(near(cSysCalTab.col("TIME"),mjds));
1104	tmptab = cSysCalTab(near(cSysCalTab.col("TIME"),mjds), nrws);
1105	//DEBUG
1106	//if (cIdx == 0) {
1107	// cerr<<"inital table retrieved"<<endl;
1108	//}
1109
1110	}
1111
1112	if (nBeam == 1) {
1113	tmptab2 = tmptab( tmptab.col("SPECTRAL_WINDOW_ID") == iIF, 1);
1114	} else {
1115	tmptab2 = tmptab( tmptab.col("SPECTRAL_WINDOW_ID") == iIF &&
1116	tmptab.col("FEED_ID") == ibeam , 1);
1117	}
1118	//cerr<<"first subtab rows="<<tmptab.nrow()<<endl;
1119	int syscalrow = tmptab2.nrow();
1120	ROArrayColumn<Float> tcalCol(tmptab2, "TCAL");
1121	ROScalarColumn<Double> tcalTimeCol(tmptab2, "TIME");
1122	if (syscalrow==0) {
1123	os << LogIO::NORMAL
1124	<<"Cannot find any matching Tcal at/near the data timestamp."
1125	<< " Set Tcal=0.0" << LogIO::POST ;
1126	} else {
1127	tcalCol.get(0, pksrec.tcal);
1128	tcalTimeCol.get(0,dtcalTime);
1129	pksrec.tcalTime = MVTime(dtcalTime/(24*3600)).string(MVTime::YMD);
1130	//DEBUG
1131	//cerr<<"cIdx:"<<cIdx<<" tcal="<<tcal<<" tcalTime="<<tcalTime<<endl;
1132	tmptab.markForDelete();
1133	tmptab2.markForDelete();
1134	}
1135	}
1136	lastmjd = pksrec.mjd;
1137	}
1138
1139	// Calibration factors (if available).
1140	pksrec.calFctr.resize(cNPol(iIF));
1141	if (cHaveCalFctr) {
1142	cCalFctrCol.get(cIdx, pksrec.calFctr);
1143	} else {
1144	pksrec.calFctr = 0.0f;
1145	}
1146
1147	// Baseline parameters (if available).
1148	if (cHaveBaseLin) {
1149	pksrec.baseLin.resize(2,cNPol(iIF));
1150	cBaseLinCol.get(cIdx, pksrec.baseLin);
1151
1152	pksrec.baseSub.resize(24,cNPol(iIF));
1153	cBaseSubCol.get(cIdx, pksrec.baseSub);
1154
1155	} else {
1156	pksrec.baseLin.resize(0,0);
1157	pksrec.baseSub.resize(0,0);
1158	}
1159
1160
1161	// Get spectral data.
1162	if (cGetSpectra) {
1163	Matrix<Float> tmpData;
1164	Matrix<Bool> tmpFlag;
1165	if (cHaveDataCol) {
1166	Matrix<Complex> tmpCmplxData;
1167	Matrix<Float> tmpReData;
1168	Matrix<Float> tmpImData;
1169	//cerr<<"reading spectra..."<<endl;
1170	//# TODO - should have a flag to user to select DATA or CORRECTED_DATA
1171	//# currently just automatically determined, --- read CORRECTED one
1172	//# if the column exist.
1173	if (cHaveCorrectedDataCol) {
1174	cCorrectedDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True);
1175	} else {
1176	cDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True);
1177	}
1178	tmpReData = real(tmpCmplxData);
1179	tmpImData = imag(tmpCmplxData);
1180	tmpData = sqrt(tmpReDatatmpReData + tmpImDatatmpImData);
1181	} else {
1182	cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
1183	}
1184	cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
1185
1186	// Transpose spectra.
1187	Int nPol = tmpData.nrow();
1188	pksrec.spectra.resize(nChan, nPol);
1189	pksrec.flagged.resize(nChan, nPol);
1190	if (cEndChan(iIF) >= cStartChan(iIF)) {
1191	// Simple transposition.
1192	for (Int ipol = 0; ipol < nPol; ipol++) {
1193	for (Int ichan = 0; ichan < nChan; ichan++) {
1194	pksrec.spectra(ichan,ipol) = tmpData(ipol,ichan);
1195	pksrec.flagged(ichan,ipol) = tmpFlag(ipol,ichan);
1196	}
1197	}
1198
1199	} else {
1200	// Transpose with inversion.
1201	Int jchan = nChan - 1;
1202	for (Int ipol = 0; ipol < nPol; ipol++) {
1203	for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
1204	pksrec.spectra(ichan,ipol) = tmpData(ipol,jchan);
1205	pksrec.flagged(ichan,ipol) = tmpFlag(ipol,jchan);
1206	}
1207	}
1208	}
1209
1210	// Row-based flagging info. (True:1, False:0)
1211	pksrec.flagrow = (cFlagRowCol(cIdx) ? 1 : 0);
1212	}
1213
1214	// Get cross-polarization data.
1215	if (cGetXPol) {
1216	//cerr<<"cGetXPol="<<cGetXPol<<endl;
1217	//cerr<<"cHaveXCalFctr="<<cHaveXCalFctr<<endl;
1218
1219	if (cHaveXCalFctr) {
1220	cXCalFctrCol.get(cIdx, pksrec.xCalFctr);
1221	} else {
1222	pksrec.xCalFctr = Complex(0.0f, 0.0f);
1223	}
1224
1225	if(!cALMA) {
1226	cDataCol.get(cIdx, pksrec.xPol, True);
1227
1228	if (cEndChan(iIF) < cStartChan(iIF)) {
1229	Complex ctmp;
1230	Int jchan = nChan - 1;
1231	for (Int ichan = 0; ichan < nChan/2; ichan++, jchan--) {
1232	ctmp = pksrec.xPol(ichan);
1233	pksrec.xPol(ichan) = pksrec.xPol(jchan);
1234	pksrec.xPol(jchan) = ctmp;
1235	}
1236	}
1237	}
1238	}
1239	/**
1240	cerr<<"scanNo="<<scanNo<<endl;
1241	cerr<<"cycleNo="<<cycleNo<<endl;
1242	cerr<<"mjd="<<mjd<<endl;
1243	cerr<<"interval="<<interval<<endl;
1244	cerr<<"fieldName="<<fieldName<<endl;
1245	cerr<<"srcNmae="<<srcName<<endl;
1246	cerr<<"srcDir="<<srcDir<<endl;
1247	cerr<<"srcPM="<<srcPM<<endl;
1248	cerr<<"srcVel="<<srcVel<<endl;
1249	cerr<<"obsMode="<<obsMode<<endl;
1250	cerr<<"IFno="<<IFno<<endl;
1251	cerr<<"refFreq="<<refFreq<<endl;
1252	cerr<<"tcal="<<tcal<<endl;
1253	cerr<<"direction="<<direction<<endl;
1254	cerr<<"scanRate="<<scanRate<<endl;
1255	cerr<<"tsys="<<tsys<<endl;
1256	cerr<<"sigma="<<sigma<<endl;
1257	cerr<<"calFctr="<<calFctr<<endl;
1258	cerr<<"baseLin="<<baseLin<<endl;
1259	cerr<<"baseSub="<<baseSub<<endl;
1260	cerr<<"spectra="<<spectra.shape()<<endl;
1261	cerr<<"flagged="<<flagged.shape()<<endl;
1262	cerr<<"xCalFctr="<<xCalFctr<<endl;
1263	cerr<<"xPol="<<xPol<<endl;
1264	**/
1265	cIdx++;
1266
1267	return 0;
1268	}
1269
1270	//--------------------------------------------------------- PKSMS2reader::read
1271
1272	// Read the next data record, just the basics.
1273
1274	Int PKSMS2reader::read(
1275	Int &IFno,
1276	Vector<Float> &tsys,
1277	Vector<Float> &calFctr,
1278	Matrix<Float> &baseLin,
1279	Matrix<Float> &baseSub,
1280	Matrix<Float> &spectra,
1281	Matrix<uChar> &flagged)
1282	{
1283	if (!cMSopen) {
1284	return 1;
1285	}
1286
1287	// Check for EOF.
1288	if (cIdx >= cNRow) {
1289	return -1;
1290	}
1291
1292	// Find the next selected beam and IF.
1293	Int ibeam;
1294	Int iIF;
1295	Int iDataDesc;
1296	while (True) {
1297	ibeam = cBeamNoCol(cIdx);
1298	//iIF = cDataDescIdCol(cIdx);
1299	iDataDesc = cDataDescIdCol(cIdx);
1300	iIF = cSpWinIdCol(iDataDesc);
1301	if (cBeams(ibeam) && cIFs(iIF)) {
1302	break;
1303	}
1304
1305	// Check for EOF.
1306	if (++cIdx >= cNRow) {
1307	return -1;
1308	}
1309	}
1310
1311	IFno = iIF + 1;
1312	// Get Tsys assuming that entries in the SYSCAL table match the main table.
1313	cTsysCol.get(cIdx, tsys, True);
1314
1315	// Calibration factors (if available).
1316	if (cHaveCalFctr) {
1317	cCalFctrCol.get(cIdx, calFctr, True);
1318	} else {
1319	calFctr.resize(cNPol(iIF));
1320	calFctr = 0.0f;
1321	}
1322
1323	// Baseline parameters (if available).
1324	if (cHaveBaseLin) {
1325	baseLin.resize(2,cNPol(iIF));
1326	cBaseLinCol.get(cIdx, baseLin);
1327
1328	baseSub.resize(24,cNPol(iIF));
1329	cBaseSubCol.get(cIdx, baseSub);
1330
1331	} else {
1332	baseLin.resize(0,0);
1333	baseSub.resize(0,0);
1334	}
1335
1336	if (cGetSpectra) {
1337	// Get spectral data.
1338	Matrix<Float> tmpData;
1339	Matrix<Bool> tmpFlag;
1340	if (cHaveDataCol) {
1341	Matrix<Complex> tmpCmplxData;
1342	cDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True);
1343	tmpData = real(tmpCmplxData);
1344	} else {
1345	cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
1346	}
1347	cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
1348
1349	// Transpose spectra.
1350	Int nChan = tmpData.ncolumn();
1351	Int nPol = tmpData.nrow();
1352	spectra.resize(nChan, nPol);
1353	flagged.resize(nChan, nPol);
1354	if (cEndChan(iIF) >= cStartChan(iIF)) {
1355	// Simple transposition.
1356	for (Int ipol = 0; ipol < nPol; ipol++) {
1357	for (Int ichan = 0; ichan < nChan; ichan++) {
1358	spectra(ichan,ipol) = tmpData(ipol,ichan);
1359	flagged(ichan,ipol) = tmpFlag(ipol,ichan);
1360	}
1361	}
1362
1363	} else {
1364	// Transpose with inversion.
1365	Int jchan = nChan - 1;
1366	for (Int ipol = 0; ipol < nPol; ipol++) {
1367	for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
1368	spectra(ichan,ipol) = tmpData(ipol,jchan);
1369	flagged(ichan,ipol) = tmpFlag(ipol,jchan);
1370	}
1371	}
1372	}
1373	}
1374
1375	cIdx++;
1376
1377	return 0;
1378	}
1379
1380	//-------------------------------------------------------- PKSMS2reader::close
1381
1382	// Close the MS.
1383
1384	void PKSMS2reader::close()
1385	{
1386	cPKSMS = MeasurementSet();
1387	cMSopen = False;
1388	}
1389
1390	//-------------------------------------------------------- PKSMS2reader::splitAntenanSelectionString
1391
1392	// split antenna selection string
1393	// delimiter is ','
1394
1395	Vector<String> PKSMS2reader::splitAntennaSelectionString( const String s )
1396	{
1397	Char delim = ',' ;
1398	Int n = s.freq( delim ) + 1 ;
1399	Vector<String> antlist ;
1400	string sl[n] ;
1401	Int numSubstr = split( s, sl, n, "," );
1402	antlist.resize( numSubstr ) ;
1403	for ( Int i = 0 ; i < numSubstr ; i++ ) {
1404	antlist[i] = String( sl[i] ) ;
1405	antlist[i].trim() ;
1406	}
1407	//cerr << "antlist = " << antlist << endl ;
1408	return antlist ;
1409	}
1410
1411	//-------------------------------------------------------- PKSMS2reader::setupAntennaList
1412
1413	// Fill cAntenna and cAntId
1414
1415	void PKSMS2reader::setupAntennaList( const String s )
1416	{
1417	LogIO os( LogOrigin( "PKSMS2reader", "setupAntennaList()", WHERE ) ) ;
1418	//cerr << "antenna specification: " << s << endl ;
1419	ROMSAntennaColumns antennaCols(cPKSMS.antenna());
1420	ROScalarColumn<String> antNames = antennaCols.name();
1421	Int nrow = antNames.nrow() ;
1422	Vector<String> antlist = splitAntennaSelectionString( s ) ;
1423	Int len = antlist.size() ;
1424	Vector<Int> AntId( len ) ;
1425	Regex re( "[0-9]+" ) ;
1426	for ( Int i = 0 ; i < len ; i++ ) {
1427	if ( antlist[i].matches( re ) ) {
1428	AntId[i] = atoi( antlist[i].c_str() ) ;
1429	if ( AntId[i] >= nrow ) {
1430	os << LogIO::SEVERE << "Antenna index out of range: " << AntId[i] << LogIO::EXCEPTION ;
1431	}
1432	}
1433	else {
1434	AntId[i] = -1 ;
1435	for ( uInt j = 0 ; j < antNames.nrow() ; j++ ) {
1436	if ( antlist[i] == antNames(j) ) {
1437	AntId[i] = j ;
1438	break ;
1439	}
1440	}
1441	if ( AntId[i] == -1 ) {
1442	os << LogIO::SEVERE << "Specified antenna name not found: " << antlist[i] << LogIO::EXCEPTION ;
1443	}
1444	}
1445	}
1446	//cerr << "AntId = " << AntId << endl ;
1447	vector<Int> uniqId ;
1448	uniqId.push_back( AntId(0) ) ;
1449	for ( uInt i = 1 ; i < AntId.size() ; i++ ) {
1450	if ( count(uniqId.begin(),uniqId.end(),AntId[i]) == 0 ) {
1451	uniqId.push_back( AntId[i] ) ;
1452	}
1453	}
1454	Vector<Int> newAntId( uniqId ) ;
1455	cAntId.assign( newAntId ) ;
1456	//cerr << "cAntId = " << cAntId << endl ;
1457	}

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

Context Navigation

source: trunk/external-alma/atnf/PKSIO/PKSMS2reader.cc @ 1880

Download in other formats: