Context Navigation

source: branches/alma/external/atnf/PKSIO/PKSMS2reader.cc@ 1488

Visit:

Last change on this file since 1488 was 1453, checked in by TakTsutsumi, 16 years ago

New Development: No

JIRA Issue: No

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: many

Test Programs: sd.scantable(), sd.scantable.save()

Put in Release Notes: N/A

Description: copied from current casapy code tree

File size: 31.9 KB

Line
1	//#---------------------------------------------------------------------------
2	//# PKSMS2reader.cc: Class to read Parkes Multibeam data from a v2 MS.
3	//#---------------------------------------------------------------------------
4	//# Copyright (C) 2000-2006
5	//# Associated Universities, Inc. Washington DC, USA.
6	//#
7	//# This library is free software; you can redistribute it and/or modify it
8	//# under the terms of the GNU Library General Public License as published by
9	//# the Free Software Foundation; either version 2 of the License, or (at your
10	//# option) any later version.
11	//#
12	//# This library is distributed in the hope that it will be useful, but
13	//# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14	//# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
15	//# License for more details.
16	//#
17	//# You should have received a copy of the GNU Library General Public License
18	//# along with this library; if not, write to the Free Software Foundation,
19	//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20	//#
21	//# Correspondence concerning AIPS++ should be addressed as follows:
22	//# Internet email: aips2-request@nrao.edu.
23	//# Postal address: AIPS++ Project Office
24	//# National Radio Astronomy Observatory
25	//# 520 Edgemont Road
26	//# Charlottesville, VA 22903-2475 USA
27	//#
28	//# $Id$
29	//#---------------------------------------------------------------------------
30	//# Original: 2000/08/03, Mark Calabretta, ATNF
31	//#---------------------------------------------------------------------------
32
33
34	// AIPS++ includes.
35	#include <casa/stdio.h>
36	#include <casa/Arrays/ArrayMath.h>
37	#include <casa/Arrays/Slice.h>
38	#include <ms/MeasurementSets/MSColumns.h>
39	#include <tables/Tables.h>
40	#include <casa/Quanta/MVTime.h>
41	#include <casa/Quanta/MVAngle.h>
42	#include <casa/BasicMath/Math.h>
43	#include <measures/Measures/MeasConvert.h>
44	#include <measures/Measures/MEpoch.h>
45	#include <measures/Measures/MeasRef.h>
46
47
48	// Parkes includes.
49	#include <atnf/pks/pks_maths.h>
50	#include <atnf/PKSIO/PKSMS2reader.h>
51
52
53	//------------------------------------------------- PKSMS2reader::PKSMS2reader
54
55	// Default constructor.
56
57	PKSMS2reader::PKSMS2reader()
58	{
59	cMSopen = False;
60	}
61
62	//------------------------------------------------ PKSMS2reader::~PKSMS2reader
63
64	PKSMS2reader::~PKSMS2reader()
65	{
66	close();
67	}
68
69	//--------------------------------------------------------- PKSMS2reader::open
70
71	// Open the MS for reading.
72
73	Int PKSMS2reader::open(
74	const String msName,
75	Vector<Bool> &beams,
76	Vector<Bool> &IFs,
77	Vector<uInt> &nChan,
78	Vector<uInt> &nPol,
79	Vector<Bool> &haveXPol,
80	Bool &haveBase,
81	Bool &haveSpectra)
82	{
83	// Check that MS is readable.
84	if (!MS::isReadable(msName)) {
85	return 1;
86	}
87
88	if (cMSopen) {
89	close();
90	}
91
92	cPKSMS = MeasurementSet(msName);
93	// taql access to the syscal table
94	cHaveSysCal = False;
95	if (cHaveSysCal=Table::isReadable(cPKSMS.sysCalTableName())) {
96	cSysCalTab = Table(cPKSMS.sysCalTableName());
97	}
98
99	cIdx = 0;
100	lastmjd = 0.0;
101	cNRow = cPKSMS.nrow();
102	cMSopen = True;
103
104	// Lock the table for read access.
105	cPKSMS.lock(False);
106
107	// Main MS table and subtable column access.
108	ROMSMainColumns msCols(cPKSMS);
109	ROMSDataDescColumns dataDescCols(cPKSMS.dataDescription());
110	ROMSFeedColumns feedCols(cPKSMS.feed());
111	ROMSFieldColumns fieldCols(cPKSMS.field());
112	ROMSPointingColumns pointingCols(cPKSMS.pointing());
113	ROMSPolarizationColumns polarizationCols(cPKSMS.polarization());
114	ROMSSourceColumns sourceCols(cPKSMS.source());
115	ROMSSpWindowColumns spWinCols(cPKSMS.spectralWindow());
116	ROMSStateColumns stateCols(cPKSMS.state());
117	ROMSSysCalColumns sysCalCols(cPKSMS.sysCal());
118	ROMSWeatherColumns weatherCols(cPKSMS.weather());
119	ROMSAntennaColumns antennaCols(cPKSMS.antenna());
120
121	// Column accessors for required columns.
122	cScanNoCol.reference(msCols.scanNumber());
123	cTimeCol.reference(msCols.time());
124	cIntervalCol.reference(msCols.interval());
125
126	cFieldIdCol.reference(msCols.fieldId());
127	cFieldNameCol.reference(fieldCols.name());
128	cFieldDelayDirCol.reference(fieldCols.delayDir());
129
130	cSrcIdCol.reference(fieldCols.sourceId());
131	cSrcId2Col.reference(sourceCols.sourceId());
132	cSrcNameCol.reference(sourceCols.name());
133	cSrcDirCol.reference(sourceCols.direction());
134	cSrcPMCol.reference(sourceCols.properMotion());
135	cSrcRestFrqCol.reference(sourceCols.restFrequency());
136
137	cStateIdCol.reference(msCols.stateId());
138	cObsModeCol.reference(stateCols.obsMode());
139	cCalCol.reference(stateCols.cal());
140	cSigStateCol.reference(stateCols.sig());
141	cRefStateCol.reference(stateCols.ref());
142	cDataDescIdCol.reference(msCols.dataDescId());
143	cSpWinIdCol.reference(dataDescCols.spectralWindowId());
144	cChanFreqCol.reference(spWinCols.chanFreq());
145
146	cWeatherTimeCol.reference(weatherCols.time());
147	cTemperatureCol.reference(weatherCols.temperature());
148	cPressureCol.reference(weatherCols.pressure());
149	cHumidityCol.reference(weatherCols.relHumidity());
150
151	cBeamNoCol.reference(msCols.feed1());
152	cPointingCol.reference(pointingCols.direction());
153	cSigmaCol.reference(msCols.sigma());
154	cNumReceptorCol.reference(feedCols.numReceptors());
155
156	// Optional columns.
157	cHaveTsys = False;
158	cHaveTcal = False;
159	if ((cHaveSrcVel = cPKSMS.source().tableDesc().isColumn("SYSVEL"))) {
160	cSrcVelCol.attach(cPKSMS.source(), "SYSVEL");
161	}
162
163	if (cHaveSysCal && (cHaveTsys = cPKSMS.sysCal().tableDesc().isColumn("TSYS"))) {
164	cTsysCol.attach(cPKSMS.sysCal(), "TSYS");
165	}
166
167	if (cHaveSysCal && (cHaveTcal = cPKSMS.sysCal().tableDesc().isColumn("TCAL"))) {
168	cTcalCol.attach(cPKSMS.sysCal(), "TCAL");
169	}
170
171	if ((cHaveCalFctr = cPKSMS.tableDesc().isColumn("CALFCTR"))) {
172	cCalFctrCol.attach(cPKSMS, "CALFCTR");
173	}
174
175	if ((cHaveBaseLin = cPKSMS.tableDesc().isColumn("BASELIN"))) {
176	cBaseLinCol.attach(cPKSMS, "BASELIN");
177	cBaseSubCol.attach(cPKSMS, "BASESUB");
178	}
179
180	// Spectral data should always be present.
181	haveSpectra = True;
182	cHaveDataCol = False;
183	cHaveCorrectedDataCol = False;
184	//String telName = antennaCols.name()(0);
185	ROMSObservationColumns observationCols(cPKSMS.observation());
186	String telName = observationCols.telescopeName()(0);
187	//cATF = (telName.contains("DA41") \|\| telName.contains("DV01"));
188	cATF = telName.contains("ATF");
189
190	if (cHaveDataCol = cPKSMS.isColumn(MSMainEnums::DATA)) {
191	if (cATF) {
192	//try to read a single baseline interferometeric data
193	//and treat it as single dish data
194	//maybe extended for ALMA commissioning later
195	cDataCol.reference(msCols.data());
196	if (cHaveCorrectedDataCol = cPKSMS.isColumn(MSMainEnums::CORRECTED_DATA)) {
197	//cerr<<"Do have CORRECTED_DATA column"<<endl;
198	cCorrectedDataCol.reference(msCols.correctedData());
199	}
200	}
201	}
202	else {
203	cFloatDataCol.reference(msCols.floatData());
204	}
205	cFlagCol.reference(msCols.flag());
206
207
208	if (cGetXPol = (cPKSMS.isColumn(MSMainEnums::DATA) && (!cATF))) {
209	if ((cHaveXCalFctr = cPKSMS.tableDesc().isColumn("XCALFCTR"))) {
210	cXCalFctrCol.attach(cPKSMS, "XCALFCTR");
211	}
212
213	cDataCol.reference(msCols.data());
214	}
215
216	// Find which beams are present in the data.
217	Vector<Int> beamNos = cBeamNoCol.getColumn();
218	Int maxBeamNo = max(beamNos) + 1;
219	beams.resize(maxBeamNo);
220
221	beams = False;
222	for (uInt irow = 0; irow < beamNos.nelements(); irow++) {
223	beams(beamNos(irow)) = True;
224	}
225
226
227	// Number of IFs.
228	//uInt nIF = dataDescCols.nrow();
229	uInt nIF =spWinCols.nrow();
230	IFs.resize(nIF);
231	IFs = True;
232
233	// Number of polarizations and channels in each IF.
234	ROScalarColumn<Int> spWinIdCol(dataDescCols.spectralWindowId());
235	ROScalarColumn<Int> numChanCol(spWinCols.numChan());
236
237	ROScalarColumn<Int> polIdCol(dataDescCols.polarizationId());
238	ROScalarColumn<Int> numPolCol(polarizationCols.numCorr());
239
240	nChan.resize(nIF);
241	nPol.resize(nIF);
242	for (uInt iIF = 0; iIF < nIF; iIF++) {
243	nChan(iIF) = numChanCol(spWinIdCol(iIF));
244	nPol(iIF) = numPolCol(polIdCol(iIF));
245	}
246
247	// Cross-polarization data present?
248	haveXPol.resize(nIF);
249	haveXPol = False;
250
251	if (cGetXPol && !(cATF)) {
252	for (Int irow = 0; irow < cNRow; irow++) {
253	if (cDataCol.isDefined(irow)) {
254	Int iIF = cDataDescIdCol(irow);
255	haveXPol(iIF) = True;
256	}
257	}
258	}
259
260
261	// Initialize member data.
262	cBeams.assign(beams);
263	cIFs.assign(IFs);
264	cNChan.assign(nChan);
265	cNPol.assign(nPol);
266	cHaveXPol.assign(haveXPol);
267
268
269	// Default channel range selection.
270	cStartChan.resize(nIF);
271	cEndChan.resize(nIF);
272	cRefChan.resize(nIF);
273
274	for (uInt iIF = 0; iIF < nIF; iIF++) {
275	cStartChan(iIF) = 1;
276	cEndChan(iIF) = cNChan(iIF);
277	cRefChan(iIF) = cNChan(iIF)/2 + 1;
278	}
279
280	Slice all;
281	cDataSel.resize(nIF);
282	cDataSel = Slicer(all, all);
283
284	cScanNo = 0;
285	cCycleNo = 1;
286	cTime = cTimeCol(0);
287
288	return 0;
289	}
290
291	//---------------------------------------------------- PKSMS2reader::getHeader
292
293	// Get parameters describing the data.
294
295	Int PKSMS2reader::getHeader(
296	String &observer,
297	String &project,
298	String &antName,
299	Vector<Double> &antPosition,
300	String &obsMode,
301	Float &equinox,
302	String &dopplerFrame,
303	Double &mjd,
304	Double &refFreq,
305	Double &bandwidth,
306	String &fluxunit)
307	{
308	if (!cMSopen) {
309	return 1;
310	}
311
312	// Observer and project.
313	ROMSObservationColumns observationCols(cPKSMS.observation());
314	observer = observationCols.observer()(0);
315	project = observationCols.project()(0);
316
317	// Antenna name and ITRF coordinates.
318	ROMSAntennaColumns antennaCols(cPKSMS.antenna());
319	antName = antennaCols.name()(0);
320	antPosition = antennaCols.position()(0);
321
322	// Observation type.
323	if (cObsModeCol.nrow()) {
324	obsMode = cObsModeCol(0);
325	if (obsMode == "\0") obsMode = "RF";
326	} else {
327	obsMode = "RF";
328	}
329
330	fluxunit = "";
331	if (cHaveDataCol) {
332	const TableRecord& keywordSet2
333	= cDataCol.columnDesc().keywordSet();
334	if(keywordSet2.isDefined("UNIT")) {
335	fluxunit = keywordSet2.asString("UNIT");
336	}
337	} else {
338	const TableRecord& keywordSet
339	= cFloatDataCol.columnDesc().keywordSet();
340	if(keywordSet.isDefined("UNIT")) {
341	fluxunit = keywordSet.asString("UNIT");
342	}
343	}
344
345	/***
346	const TableRecord& keywordSet
347	= cFloatDataCol.columnDesc().keywordSet();
348	if(keywordSet.isDefined("UNIT")) {
349	fluxunit = keywordSet.asString("UNIT");
350	}
351	***/
352	// Coordinate equinox.
353	ROMSPointingColumns pointingCols(cPKSMS.pointing());
354	String dirref = pointingCols.direction().keywordSet().asRecord("MEASINFO").
355	asString("Ref");
356	sscanf(dirref.chars()+1, "%f", &equinox);
357
358	// Frequency/velocity reference frame.
359	ROMSSpWindowColumns spWinCols(cPKSMS.spectralWindow());
360	dopplerFrame = MFrequency::showType(spWinCols.measFreqRef()(0));
361
362	// Translate to FITS standard names.
363	if (dopplerFrame == "TOPO") {
364	dopplerFrame = "TOPOCENT";
365	} else if (dopplerFrame == "GEO") {
366	dopplerFrame = "GEOCENTR";
367	} else if (dopplerFrame == "BARY") {
368	dopplerFrame = "BARYCENT";
369	} else if (dopplerFrame == "GALACTO") {
370	dopplerFrame = "GALACTOC";
371	} else if (dopplerFrame == "LGROUP") {
372	dopplerFrame = "LOCALGRP";
373	} else if (dopplerFrame == "CMB") {
374	dopplerFrame = "CMBDIPOL";
375	} else if (dopplerFrame == "REST") {
376	dopplerFrame = "SOURCE";
377	}
378
379	// MJD at start of observation.
380	mjd = cTimeCol(0)/86400.0;
381
382	// Reference frequency and bandwidth.
383	refFreq = spWinCols.refFrequency()(0);
384	bandwidth = spWinCols.totalBandwidth()(0);
385
386	return 0;
387	}
388
389	//-------------------------------------------------- PKSMS2reader::getFreqInfo
390
391	// Get frequency parameters for each IF.
392
393	Int PKSMS2reader::getFreqInfo(
394	Vector<Double> &startFreq,
395	Vector<Double> &endFreq)
396	{
397	uInt nIF = cIFs.nelements();
398	startFreq.resize(nIF);
399	endFreq.resize(nIF);
400
401	for (uInt iIF = 0; iIF < nIF; iIF++) {
402	Vector<Double> chanFreq = cChanFreqCol(iIF);
403
404	Int nChan = chanFreq.nelements();
405	startFreq(iIF) = chanFreq(0);
406	endFreq(iIF) = chanFreq(nChan-1);
407	}
408
409	return 0;
410	}
411
412	//------------------------------------------------------- PKSMS2reader::select
413
414	// Set data selection by beam number and channel.
415
416	uInt PKSMS2reader::select(
417	const Vector<Bool> beamSel,
418	const Vector<Bool> IFsel,
419	const Vector<Int> startChan,
420	const Vector<Int> endChan,
421	const Vector<Int> refChan,
422	const Bool getSpectra,
423	const Bool getXPol,
424	const Bool getFeedPos)
425	{
426	if (!cMSopen) {
427	return 1;
428	}
429
430	// Beam selection.
431	uInt nBeam = cBeams.nelements();
432	uInt nBeamSel = beamSel.nelements();
433	for (uInt ibeam = 0; ibeam < nBeam; ibeam++) {
434	if (ibeam < nBeamSel) {
435	cBeams(ibeam) = beamSel(ibeam);
436	} else {
437	cBeams(ibeam) = False;
438	}
439	}
440
441	uInt nIF = cIFs.nelements();
442	uInt maxNChan = 0;
443	for (uInt iIF = 0; iIF < nIF; iIF++) {
444	// IF selection.
445	if (iIF < IFsel.nelements()) {
446	cIFs(iIF) = IFsel(iIF);
447	} else {
448	cIFs(iIF) = False;
449	}
450
451	if (!cIFs(iIF)) continue;
452
453
454	// Channel selection.
455	if (iIF < startChan.nelements()) {
456	cStartChan(iIF) = startChan(iIF);
457
458	if (cStartChan(iIF) <= 0) {
459	cStartChan(iIF) += cNChan(iIF);
460	} else if (cStartChan(iIF) > Int(cNChan(iIF))) {
461	cStartChan(iIF) = cNChan(iIF);
462	}
463	}
464
465	if (iIF < endChan.nelements()) {
466	cEndChan(iIF) = endChan(iIF);
467
468	if (cEndChan(iIF) <= 0) {
469	cEndChan(iIF) += cNChan(iIF);
470	} else if (cEndChan(iIF) > Int(cNChan(iIF))) {
471	cEndChan(iIF) = cNChan(iIF);
472	}
473	}
474
475	if (iIF < refChan.nelements()) {
476	cRefChan(iIF) = refChan(iIF);
477	} else {
478	cRefChan(iIF) = cStartChan(iIF);
479	if (cStartChan(iIF) <= cEndChan(iIF)) {
480	cRefChan(iIF) += (cEndChan(iIF) - cStartChan(iIF) + 1)/2;
481	} else {
482	cRefChan(iIF) -= (cStartChan(iIF) - cEndChan(iIF) + 1)/2;
483	}
484	}
485
486	uInt nChan = abs(cEndChan(iIF) - cStartChan(iIF)) + 1;
487	if (maxNChan < nChan) {
488	maxNChan = nChan;
489	}
490
491	// Inverted Slices are not allowed.
492	Slice outPols;
493	Slice outChans(min(cStartChan(iIF),cEndChan(iIF))-1, nChan);
494	cDataSel(iIF) = Slicer(outPols, outChans);
495	}
496
497	// Get spectral data?
498	cGetSpectra = getSpectra;
499
500	// Get cross-polarization data?
501	cGetXPol = cGetXPol && getXPol;
502
503	// Get feed positions? (Not available.)
504	cGetFeedPos = False;
505
506	return maxNChan;
507	}
508
509	//---------------------------------------------------- PKSMS2reader::findRange
510
511	// Find the range of the data in time and position.
512
513	Int PKSMS2reader::findRange(
514	Int &nRow,
515	Int &nSel,
516	Vector<Double> &timeSpan,
517	Matrix<Double> &positions)
518	{
519	if (!cMSopen) {
520	return 1;
521	}
522
523	nRow = cNRow;
524
525	// Find the number of rows selected.
526	nSel = 0;
527	Vector<Bool> sel(nRow);
528	for (Int irow = 0; irow < nRow; irow++) {
529	if ((sel(irow) = cBeams(cBeamNoCol(irow)) &&
530	cIFs(cDataDescIdCol(irow)))) {
531	nSel++;
532	}
533	}
534
535	// Find the time range (s).
536	timeSpan.resize(2);
537	timeSpan(0) = cTimeCol(0);
538	timeSpan(1) = cTimeCol(nRow-1);
539
540	// Retrieve positions for selected data.
541	Int isel = 0;
542	positions.resize(2,nSel);
543	for (Int irow = 0; irow < nRow; irow++) {
544	if (sel(irow)) {
545	Matrix<Double> pointingDir = cPointingCol(cFieldIdCol(irow));
546	positions.column(isel++) = pointingDir.column(0);
547	}
548	}
549
550	return 0;
551	}
552
553	//--------------------------------------------------------- PKSMS2reader::read
554
555	// Read the next data record.
556
557	Int PKSMS2reader::read(
558	Int &scanNo,
559	Int &cycleNo,
560	Double &mjd,
561	Double &interval,
562	String &fieldName,
563	String &srcName,
564	Vector<Double> &srcDir,
565	Vector<Double> &srcPM,
566	Double &srcVel,
567	String &obsMode,
568	Int &IFno,
569	Double &refFreq,
570	Double &bandwidth,
571	Double &freqInc,
572	Vector<Double> &restFreq,
573	Vector<Float> &tcal,
574	String &tcalTime,
575	Float &azimuth,
576	Float &elevation,
577	Float &parAngle,
578	Float &focusAxi,
579	Float &focusTan,
580	Float &focusRot,
581	Float &temperature,
582	Float &pressure,
583	Float &humidity,
584	Float &windSpeed,
585	Float &windAz,
586	Int &refBeam,
587	Int &beamNo,
588	Vector<Double> &direction,
589	Vector<Double> &scanRate,
590	Vector<Float> &tsys,
591	Vector<Float> &sigma,
592	Vector<Float> &calFctr,
593	Matrix<Float> &baseLin,
594	Matrix<Float> &baseSub,
595	Matrix<Float> &spectra,
596	Matrix<uChar> &flagged,
597	Complex &xCalFctr,
598	Vector<Complex> &xPol)
599	{
600	if (!cMSopen) {
601	return 1;
602	}
603
604	// Check for EOF.
605	if (cIdx >= cNRow) {
606	return -1;
607	}
608
609	// Find the next selected beam and IF.
610	Int ibeam;
611	Int iIF;
612	Int iDataDesc;
613
614	while (True) {
615	ibeam = cBeamNoCol(cIdx);
616	iDataDesc = cDataDescIdCol(cIdx);
617	iIF =cSpWinIdCol(iDataDesc);
618	if (cBeams(ibeam) && cIFs(iIF)) {
619	break;
620	}
621
622	// Check for EOF.
623	if (++cIdx >= cNRow) {
624	return -1;
625	}
626	}
627	// Renumerate scan no. Here still is 1-based
628	//scanNo = cScanNoCol(cIdx) - cScanNoCol(0) + 1;
629	scanNo = cScanNoCol(cIdx);
630
631	if (scanNo != cScanNo) {
632	// Start of new scan.
633	cScanNo = scanNo;
634	cCycleNo = 1;
635	cTime = cTimeCol(cIdx);
636	}
637
638	Double time = cTimeCol(cIdx);
639	mjd = time/86400.0;
640	interval = cIntervalCol(cIdx);
641
642	// Reconstruct the integration cycle number; due to small latencies the
643	// integration time is usually slightly less than the time between cycles,
644	// resetting cTime will prevent the difference from accumulating.
645	cCycleNo += nint((time - cTime)/interval);
646	cycleNo = cCycleNo;
647	cTime = time;
648
649	Int fieldId = cFieldIdCol(cIdx);
650	fieldName = cFieldNameCol(fieldId);
651
652	Int srcId = cSrcIdCol(fieldId);
653	//For source with multiple spectral window setting, this is not
654	// correct. Source name of srcId may not be at 'srcId'th row of SrcNameCol
655	//srcName = cSrcNameCol(srcId);
656	for (uInt irow = 0; irow < cSrcId2Col.nrow(); irow++) {
657	if (cSrcId2Col(irow) == srcId) {
658	srcName = cSrcNameCol(irow);
659	}
660	}
661
662	srcDir = cSrcDirCol(srcId);
663	srcPM = cSrcPMCol(srcId);
664
665	// Systemic velocity.
666	if (!cHaveSrcVel \|\| cATF) {
667	srcVel = 0.0f;
668	} else {
669	srcVel = cSrcVelCol(srcId)(IPosition(1,0));
670	}
671
672	ROMSAntennaColumns antennaCols(cPKSMS.antenna());
673	String telescope = antennaCols.name()(0);
674	Bool cGBT = telescope.contains("GBT");
675	// Observation type.
676	// check if State Table exist
677	//Bool cHaveStateTab=Table::isReadable(cPKSMS.stateTableName());
678	Int stateId = 0;
679	Int StateNRow = 0;
680	StateNRow=cObsModeCol.nrow();
681	if (Table::isReadable(cPKSMS.stateTableName())) {
682	obsMode = " ";
683	if (StateNRow > 0) {
684	stateId = cStateIdCol(cIdx);
685	if (stateId == -1) {
686	//obsMode = " ";
687	} else {
688	obsMode = cObsModeCol(stateId);
689	Bool sigState =cSigStateCol(stateId);
690	Bool refState =cRefStateCol(stateId);
691	//DEBUG
692	//cerr <<"stateid="<<stateId<<" obsmode="<<obsMode<<endl;
693	if (cGBT) {
694	// split the obsMode string and append a proper label
695	// (these are GBT specific)
696	int epos = obsMode.find_first_of(':');
697	int nextpos = obsMode.find_first_of(':',epos+1);
698	string obsMode1 = obsMode.substr(0,epos);
699	string obsMode2 = obsMode.substr(epos+1,nextpos-epos-1);
700
701	//cerr <<"obsMode2= "<<obsMode2<<endl;
702	if (!srcName.contains("_ps")
703	&&!srcName.contains("_psr")
704	&&!srcName.contains("_nod")
705	&&!srcName.contains("_fs")
706	&&!srcName.contains("_fsr")) {
707	// if Nod mode observation , append '_nod'
708	if (obsMode1 == "Nod") {
709	srcName.append("_nod");
710	} else if (obsMode1 == "OffOn") {
711	// for GBT position switch observations (OffOn or OnOff)
712	if (obsMode2 == "PSWITCHON") srcName.append("_ps");
713	if (obsMode2 == "PSWITCHOFF") srcName.append("_psr");
714	} else {
715	if (obsMode2 == "FSWITCH") {
716	// for GBT frequency switch mode
717	if (sigState) srcName.append("_fs");
718	if (refState) srcName.append("_fsr");
719	}
720	}
721	}
722	}
723	}
724	}
725	}
726	// CAL state
727	// this should be apply just for GBT data?
728	Double Cal;
729	if (stateId==-1 \|\| StateNRow==0) {
730	Cal = 0;
731	} else {
732	Cal = cCalCol(stateId);
733	}
734	if (cGBT) {
735	if (Cal > 0 && !srcName.contains("_calon")) {
736	srcName.append("_calon");
737	}
738	}
739
740	IFno = iIF + 1;
741	Int nChan = abs(cEndChan(iIF) - cStartChan(iIF)) + 1;
742
743	// Minimal handling on continuum data.
744	Vector<Double> chanFreq = cChanFreqCol(iIF);
745	if (nChan == 1) {
746	freqInc = chanFreq(0);
747	refFreq = chanFreq(0);
748	restFreq = 0.0f;
749	} else {
750
751	if (cStartChan(iIF) <= cEndChan(iIF)) {
752	freqInc = chanFreq(1) - chanFreq(0);
753	} else {
754	freqInc = chanFreq(0) - chanFreq(1);
755	}
756	refFreq = chanFreq(cRefChan(iIF)-1);
757	Bool HaveSrcRestFreq= cSrcRestFrqCol.isDefined(srcId);
758	if (HaveSrcRestFreq) {
759	//restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0));
760	restFreq = cSrcRestFrqCol(srcId);
761	} else {
762	restFreq = 0.0f;
763	}
764	}
765	bandwidth = abs(freqInc * nChan);
766
767	tcal.resize(cNPol(iIF));
768	tcal = 0.0f;
769	tcalTime = "";
770	//azimuth = 0.0f;
771	//elevation = 0.0f;
772	parAngle = 0.0f;
773	focusAxi = 0.0f;
774	focusTan = 0.0f;
775	focusRot = 0.0f;
776
777	// Find the appropriate entry in the WEATHER subtable.
778	//Bool cHaveStateTab=Table::isReadable(cPKSMS.stateTableName());
779	Bool cHaveWeatherTab = Table::isReadable(cPKSMS.weatherTableName());
780	Int weatherIdx=-1;
781	if (cHaveWeatherTab) {
782	Vector<Double> wTimes = cWeatherTimeCol.getColumn();
783	for (weatherIdx = wTimes.nelements()-1; weatherIdx >= 0; weatherIdx--) {
784	if (cWeatherTimeCol(weatherIdx) <= time) {
785	break;
786	}
787	}
788	}
789	if (weatherIdx < 0 \|\| !cHaveWeatherTab) {
790	// No appropriate WEATHER entry.
791	pressure = 0.0f;
792	humidity = 0.0f;
793	temperature = 0.0f;
794	} else {
795	pressure = cPressureCol(weatherIdx);
796	humidity = cHumidityCol(weatherIdx);
797	temperature = cTemperatureCol(weatherIdx);
798	}
799
800	windSpeed = 0.0f;
801	windAz = 0.0f;
802
803	refBeam = 0;
804	beamNo = ibeam + 1;
805
806	//Matrix<Double> pointingDir = cPointingCol(fieldId);
807	//pointingDir = cPointingCol(fieldId);
808	//direction = pointingDir.column(0);
809	//uInt ncols = pointingDir.ncolumn();
810	//if (ncols == 1) {
811	// scanRate = 0.0f;
812	//} else {
813	// scanRate = pointingDir.column(1);
814	//}
815
816	// Get direction from FIELD table
817	// here, assume direction to be the field direction not pointing
818	Matrix<Double> delayDir = cFieldDelayDirCol(fieldId);
819	direction = delayDir.column(0);
820	uInt ncols = delayDir.ncolumn();
821	if (ncols == 1) {
822	scanRate = 0.0f;
823	} else {
824	scanRate = delayDir.column(1);
825	}
826
827	// caluculate azimuth and elevation
828	// first, get the reference frame
829	MVPosition mvpos(antennaCols.position()(0));
830	MPosition mp(mvpos);
831	Quantum<Double> qt(time,"s");
832	MVEpoch mvt(qt);
833	MEpoch me(mvt);
834	MeasFrame frame(mp, me);
835	//
836	ROMSFieldColumns fldCols(cPKSMS.field());
837	Vector<MDirection> vmd(1);
838	MDirection md;
839	fldCols.delayDirMeasCol().get(fieldId,vmd);
840	md = vmd[0];
841	//Vector<Double> dircheck = md.getAngle("rad").getValue();
842	//cerr<<"dircheck="<<dircheck<<endl;
843
844	Vector<Double> azel =
845	MDirection::Convert(md, MDirection::Ref(MDirection::AZEL,
846	frame)
847	)().getAngle("rad").getValue();
848	//cerr<<"azel="<<azel<<endl;
849	azimuth = azel[0];
850	elevation = azel[1];
851
852	// Get Tsys assuming that entries in the SYSCAL table match the main table.
853	if (cHaveTsys) {
854	Int nTsysColRow = cTsysCol.nrow();
855	if (nTsysColRow != cNRow) {
856	cHaveTsys=0;
857	}
858	}
859	if (cHaveTsys) {
860	cTsysCol.get(cIdx, tsys, True);
861	} else {
862	Int numReceptor;
863	cNumReceptorCol.get(0, numReceptor);
864	tsys.resize(numReceptor);
865	tsys = 1.0f;
866	}
867	cSigmaCol.get(cIdx, sigma, True);
868
869	//get Tcal if available
870	if (cHaveTcal) {
871	Int nTcalColRow = cTcalCol.nrow();
872	uInt nBeam = cBeams.nelements();
873	uInt nIF = cIFs.nelements();
874	uInt nrws = nBeam * nIF;
875	if (nTcalColRow > 0) {
876	// find tcal match with the data with the data time stamp
877	Double mjds = mjd(243600);
878	Double dtcalTime;
879	if ( mjd > lastmjd \|\| cIdx==0 ) {
880	//Table tmptab = cSysCalTab(near(cSysCalTab.col("TIME"),mjds));
881	tmptab = cSysCalTab(near(cSysCalTab.col("TIME"),mjds), nrws);
882	//DEBUG
883	//if (cIdx == 0) {
884	// cerr<<"inital table retrieved"<<endl;
885	//}
886
887	}
888
889	if (nBeam == 1) {
890	tmptab2 = tmptab( tmptab.col("SPECTRAL_WINDOW_ID") == iIF, 1);
891	} else {
892	tmptab2 = tmptab( tmptab.col("SPECTRAL_WINDOW_ID") == iIF &&
893	tmptab.col("FEED_ID") == ibeam , 1);
894	}
895	//cerr<<"first subtab rows="<<tmptab.nrow()<<endl;
896	int syscalrow = tmptab2.nrow();
897	ROArrayColumn<Float> tcalCol(tmptab2, "TCAL");
898	ROScalarColumn<Double> tcalTimeCol(tmptab2, "TIME");
899	if (syscalrow==0) {
900	cerr<<"Cannot find any matching Tcal at/near the data timestamp."
901	<< " Set Tcal=0.0"<<endl;
902	} else {
903	tcalCol.get(0, tcal);
904	tcalTimeCol.get(0,dtcalTime);
905	tcalTime = MVTime(dtcalTime/(24*3600)).string(MVTime::YMD);
906	//DEBUG
907	//cerr<<"cIdx:"<<cIdx<<" tcal="<<tcal<<" tcalTime="<<tcalTime<<endl;
908	tmptab.markForDelete();
909	tmptab2.markForDelete();
910	}
911	}
912	lastmjd = mjd;
913	}
914
915	// Calibration factors (if available).
916	calFctr.resize(cNPol(iIF));
917	if (cHaveCalFctr) {
918	cCalFctrCol.get(cIdx, calFctr);
919	} else {
920	calFctr = 0.0f;
921	}
922
923	// Baseline parameters (if available).
924	if (cHaveBaseLin) {
925	baseLin.resize(2,cNPol(iIF));
926	cBaseLinCol.get(cIdx, baseLin);
927
928	baseSub.resize(9,cNPol(iIF));
929	cBaseSubCol.get(cIdx, baseSub);
930
931	} else {
932	baseLin.resize(0,0);
933	baseSub.resize(0,0);
934	}
935	// Get spectral data.
936	if (cGetSpectra) {
937	Matrix<Float> tmpData;
938	Matrix<Bool> tmpFlag;
939	if (cHaveDataCol) {
940	Matrix<Complex> tmpCmplxData;
941	Matrix<Float> tmpReData;
942	Matrix<Float> tmpImData;
943	//cerr<<"reading spectra..."<<endl;
944	//# TODO - should have a flag to user to select DATA or CORRECTED_DATA
945	//# currently just automatically determined, --- read CORRECTED one
946	//# if the column exist.
947	if (cHaveCorrectedDataCol) {
948	cCorrectedDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True);
949	} else {
950	cDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True);
951	}
952	tmpReData = real(tmpCmplxData);
953	tmpImData = imag(tmpCmplxData);
954	tmpData = sqrt(tmpReDatatmpReData + tmpImDatatmpImData);
955	} else {
956	cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
957	}
958	cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
959
960	// Transpose spectra.
961	Int nPol = tmpData.nrow();
962	spectra.resize(nChan, nPol);
963	flagged.resize(nChan, nPol);
964	if (cEndChan(iIF) >= cStartChan(iIF)) {
965	// Simple transposition.
966	for (Int ipol = 0; ipol < nPol; ipol++) {
967	for (Int ichan = 0; ichan < nChan; ichan++) {
968	spectra(ichan,ipol) = tmpData(ipol,ichan);
969	flagged(ichan,ipol) = tmpFlag(ipol,ichan);
970	}
971	}
972
973	} else {
974	// Transpose with inversion.
975	Int jchan = nChan - 1;
976	for (Int ipol = 0; ipol < nPol; ipol++) {
977	for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
978	spectra(ichan,ipol) = tmpData(ipol,jchan);
979	flagged(ichan,ipol) = tmpFlag(ipol,jchan);
980	}
981	}
982	}
983	}
984
985	// Get cross-polarization data.
986	if (cGetXPol) {
987	//cerr<<"cGetXPol="<<cGetXPol<<endl;
988	//cerr<<"cHaveXCalFctr="<<cHaveXCalFctr<<endl;
989
990	if (cHaveXCalFctr) {
991	cXCalFctrCol.get(cIdx, xCalFctr);
992	} else {
993	xCalFctr = Complex(0.0f, 0.0f);
994	}
995
996	if(!cATF) {
997	cDataCol.get(cIdx, xPol, True);
998
999	if (cEndChan(iIF) < cStartChan(iIF)) {
1000	Complex ctmp;
1001	Int jchan = nChan - 1;
1002	for (Int ichan = 0; ichan < nChan/2; ichan++, jchan--) {
1003	ctmp = xPol(ichan);
1004	xPol(ichan) = xPol(jchan);
1005	xPol(jchan) = ctmp;
1006	}
1007	}
1008	}
1009	}
1010	/**
1011	cerr<<"scanNo="<<scanNo<<endl;
1012	cerr<<"cycleNo="<<cycleNo<<endl;
1013	cerr<<"mjd="<<mjd<<endl;
1014	cerr<<"interval="<<interval<<endl;
1015	cerr<<"fieldName="<<fieldName<<endl;
1016	cerr<<"srcNmae="<<srcName<<endl;
1017	cerr<<"srcDir="<<srcDir<<endl;
1018	cerr<<"srcPM="<<srcPM<<endl;
1019	cerr<<"srcVel="<<srcVel<<endl;
1020	cerr<<"obsMode="<<obsMode<<endl;
1021	cerr<<"IFno="<<IFno<<endl;
1022	cerr<<"refFreq="<<refFreq<<endl;
1023	cerr<<"tcal="<<tcal<<endl;
1024	cerr<<"direction="<<direction<<endl;
1025	cerr<<"scanRate="<<scanRate<<endl;
1026	cerr<<"tsys="<<tsys<<endl;
1027	cerr<<"sigma="<<sigma<<endl;
1028	cerr<<"calFctr="<<calFctr<<endl;
1029	cerr<<"baseLin="<<baseLin<<endl;
1030	cerr<<"baseSub="<<baseSub<<endl;
1031	cerr<<"spectra="<<spectra.shape()<<endl;
1032	cerr<<"flagged="<<flagged.shape()<<endl;
1033	cerr<<"xCalFctr="<<xCalFctr<<endl;
1034	cerr<<"xPol="<<xPol<<endl;
1035	**/
1036	cIdx++;
1037
1038	return 0;
1039	}
1040
1041	//--------------------------------------------------------- PKSMS2reader::read
1042
1043	// Read the next data record, just the basics.
1044
1045	Int PKSMS2reader::read(
1046	Int &IFno,
1047	Vector<Float> &tsys,
1048	Vector<Float> &calFctr,
1049	Matrix<Float> &baseLin,
1050	Matrix<Float> &baseSub,
1051	Matrix<Float> &spectra,
1052	Matrix<uChar> &flagged)
1053	{
1054	if (!cMSopen) {
1055	return 1;
1056	}
1057
1058	// Check for EOF.
1059	if (cIdx >= cNRow) {
1060	return -1;
1061	}
1062
1063	// Find the next selected beam and IF.
1064	Int ibeam;
1065	Int iIF;
1066	Int iDataDesc;
1067	while (True) {
1068	ibeam = cBeamNoCol(cIdx);
1069	//iIF = cDataDescIdCol(cIdx);
1070	iDataDesc = cDataDescIdCol(cIdx);
1071	iIF = cSpWinIdCol(iDataDesc);
1072	if (cBeams(ibeam) && cIFs(iIF)) {
1073	break;
1074	}
1075
1076	// Check for EOF.
1077	if (++cIdx >= cNRow) {
1078	return -1;
1079	}
1080	}
1081
1082	IFno = iIF + 1;
1083	// Get Tsys assuming that entries in the SYSCAL table match the main table.
1084	cTsysCol.get(cIdx, tsys, True);
1085
1086	// Calibration factors (if available).
1087	if (cHaveCalFctr) {
1088	cCalFctrCol.get(cIdx, calFctr, True);
1089	} else {
1090	calFctr.resize(cNPol(iIF));
1091	calFctr = 0.0f;
1092	}
1093
1094	// Baseline parameters (if available).
1095	if (cHaveBaseLin) {
1096	baseLin.resize(2,cNPol(iIF));
1097	cBaseLinCol.get(cIdx, baseLin);
1098
1099	baseSub.resize(9,cNPol(iIF));
1100	cBaseSubCol.get(cIdx, baseSub);
1101
1102	} else {
1103	baseLin.resize(0,0);
1104	baseSub.resize(0,0);
1105	}
1106
1107	if (cGetSpectra) {
1108	// Get spectral data.
1109	Matrix<Float> tmpData;
1110	Matrix<Bool> tmpFlag;
1111	if (cHaveDataCol) {
1112	Matrix<Complex> tmpCmplxData;
1113	cDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True);
1114	tmpData = real(tmpCmplxData);
1115	} else {
1116	cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
1117	}
1118	cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
1119
1120	// Transpose spectra.
1121	Int nChan = tmpData.ncolumn();
1122	Int nPol = tmpData.nrow();
1123	spectra.resize(nChan, nPol);
1124	flagged.resize(nChan, nPol);
1125	if (cEndChan(iIF) >= cStartChan(iIF)) {
1126	// Simple transposition.
1127	for (Int ipol = 0; ipol < nPol; ipol++) {
1128	for (Int ichan = 0; ichan < nChan; ichan++) {
1129	spectra(ichan,ipol) = tmpData(ipol,ichan);
1130	flagged(ichan,ipol) = tmpFlag(ipol,ichan);
1131	}
1132	}
1133
1134	} else {
1135	// Transpose with inversion.
1136	Int jchan = nChan - 1;
1137	for (Int ipol = 0; ipol < nPol; ipol++) {
1138	for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
1139	spectra(ichan,ipol) = tmpData(ipol,jchan);
1140	flagged(ichan,ipol) = tmpFlag(ipol,jchan);
1141	}
1142	}
1143	}
1144	}
1145
1146	cIdx++;
1147
1148	return 0;
1149	}
1150
1151	//-------------------------------------------------------- PKSMS2reader::close
1152
1153	// Close the MS.
1154
1155	void PKSMS2reader::close()
1156	{
1157	cPKSMS = MeasurementSet();
1158	cMSopen = False;
1159	}

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

Download in other formats: