source: trunk/external/atnf/PKSIO/PKSMS2reader.cc@ 1423

Last change on this file since 1423 was 1399, checked in by Malte Marquarding, 17 years ago

Mark C added brightness unit to getHeader()
File size: 22.3 KB
RevLine 
[1325]1//#---------------------------------------------------------------------------
2//# PKSMS2reader.cc: Class to read Parkes Multibeam data from a v2 MS.
3//#---------------------------------------------------------------------------
[1399]4//# Copyright (C) 2000-2007
[1325]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//#
[1399]28//# $Id: PKSMS2reader.cc,v 19.12 2007/11/12 03:37:56 cal103 Exp $
[1325]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
41// Parkes includes.
42#include <atnf/pks/pks_maths.h>
43#include <atnf/PKSIO/PKSMS2reader.h>
44
45
46//------------------------------------------------- PKSMS2reader::PKSMS2reader
47
48// Default constructor.
49
50PKSMS2reader::PKSMS2reader()
51{
52 cMSopen = False;
53}
54
55//------------------------------------------------ PKSMS2reader::~PKSMS2reader
56
57PKSMS2reader::~PKSMS2reader()
58{
59 close();
60}
61
62//--------------------------------------------------------- PKSMS2reader::open
63
64// Open the MS for reading.
65
66Int PKSMS2reader::open(
67 const String msName,
68 Vector<Bool> &beams,
69 Vector<Bool> &IFs,
70 Vector<uInt> &nChan,
71 Vector<uInt> &nPol,
72 Vector<Bool> &haveXPol,
73 Bool &haveBase,
74 Bool &haveSpectra)
75{
76 // Check that MS is readable.
77 if (!MS::isReadable(msName)) {
78 return 1;
79 }
80
81 if (cMSopen) {
82 close();
83 }
84
85 cPKSMS = MeasurementSet(msName);
86 cIdx = 0;
87 cNRow = cPKSMS.nrow();
88 cMSopen = True;
89
90 // Lock the table for read access.
91 cPKSMS.lock(False);
92
93 // Main MS table and subtable column access.
94 ROMSMainColumns msCols(cPKSMS);
95 ROMSDataDescColumns dataDescCols(cPKSMS.dataDescription());
96 ROMSFeedColumns feedCols(cPKSMS.feed());
97 ROMSFieldColumns fieldCols(cPKSMS.field());
98 ROMSPointingColumns pointingCols(cPKSMS.pointing());
99 ROMSPolarizationColumns polarizationCols(cPKSMS.polarization());
100 ROMSSourceColumns sourceCols(cPKSMS.source());
101 ROMSSpWindowColumns spWinCols(cPKSMS.spectralWindow());
102 ROMSStateColumns stateCols(cPKSMS.state());
103 ROMSSysCalColumns sysCalCols(cPKSMS.sysCal());
104 ROMSWeatherColumns weatherCols(cPKSMS.weather());
105
106 // Column accessors for required columns.
107 cScanNoCol.reference(msCols.scanNumber());
108 cTimeCol.reference(msCols.time());
109 cIntervalCol.reference(msCols.interval());
110
111 cFieldIdCol.reference(msCols.fieldId());
112 cFieldNameCol.reference(fieldCols.name());
113
114 cSrcIdCol.reference(fieldCols.sourceId());
115 cSrcNameCol.reference(sourceCols.name());
116 cSrcDirCol.reference(sourceCols.direction());
117 cSrcPMCol.reference(sourceCols.properMotion());
118 cSrcRestFrqCol.reference(sourceCols.restFrequency());
119
120 cStateIdCol.reference(msCols.stateId());
121 cObsModeCol.reference(stateCols.obsMode());
122
123 cDataDescIdCol.reference(msCols.dataDescId());
124 cChanFreqCol.reference(spWinCols.chanFreq());
125
126 cWeatherTimeCol.reference(weatherCols.time());
127 cTemperatureCol.reference(weatherCols.temperature());
128 cPressureCol.reference(weatherCols.pressure());
129 cHumidityCol.reference(weatherCols.relHumidity());
130
131 cBeamNoCol.reference(msCols.feed1());
132 cPointingCol.reference(pointingCols.direction());
133 cSigmaCol.reference(msCols.sigma());
134 cNumReceptorCol.reference(feedCols.numReceptors());
135
136 // Optional columns.
137 if ((cHaveSrcVel = cPKSMS.source().tableDesc().isColumn("SYSVEL"))) {
138 cSrcVelCol.attach(cPKSMS.source(), "SYSVEL");
139 }
140
141 if ((cHaveTsys = cPKSMS.sysCal().tableDesc().isColumn("TSYS"))) {
142 cTsysCol.attach(cPKSMS.sysCal(), "TSYS");
143 }
144
145 if ((cHaveCalFctr = cPKSMS.tableDesc().isColumn("CALFCTR"))) {
146 cCalFctrCol.attach(cPKSMS, "CALFCTR");
147 }
148
149 if ((cHaveBaseLin = cPKSMS.tableDesc().isColumn("BASELIN"))) {
150 cBaseLinCol.attach(cPKSMS, "BASELIN");
151 cBaseSubCol.attach(cPKSMS, "BASESUB");
152 }
153
154 // Spectral data should always be present.
155 haveSpectra = True;
156 cFloatDataCol.reference(msCols.floatData());
157 cFlagCol.reference(msCols.flag());
158
159 if ((cGetXPol = cPKSMS.isColumn(MSMainEnums::DATA))) {
160 if ((cHaveXCalFctr = cPKSMS.tableDesc().isColumn("XCALFCTR"))) {
161 cXCalFctrCol.attach(cPKSMS, "XCALFCTR");
162 }
163
164 cDataCol.reference(msCols.data());
165 }
166
167 // Find which beams are present in the data.
168 Vector<Int> beamNos = cBeamNoCol.getColumn();
169 Int maxBeamNo = max(beamNos) + 1;
170 beams.resize(maxBeamNo);
171
172 beams = False;
173 for (uInt irow = 0; irow < beamNos.nelements(); irow++) {
174 beams(beamNos(irow)) = True;
175 }
176
177
178 // Number of IFs.
179 uInt nIF = dataDescCols.nrow();
180 IFs.resize(nIF);
181 IFs = True;
182
183 // Number of polarizations and channels in each IF.
184 ROScalarColumn<Int> spWinIdCol(dataDescCols.spectralWindowId());
185 ROScalarColumn<Int> numChanCol(spWinCols.numChan());
186
187 ROScalarColumn<Int> polIdCol(dataDescCols.polarizationId());
188 ROScalarColumn<Int> numPolCol(polarizationCols.numCorr());
189
190 nChan.resize(nIF);
191 nPol.resize(nIF);
192 for (uInt iIF = 0; iIF < nIF; iIF++) {
193 nChan(iIF) = numChanCol(spWinIdCol(iIF));
194 nPol(iIF) = numPolCol(polIdCol(iIF));
195 }
196
197 // Cross-polarization data present?
198 haveXPol.resize(nIF);
199 haveXPol = False;
200
201 if (cGetXPol) {
202 for (Int irow = 0; irow < cNRow; irow++) {
203 if (cDataCol.isDefined(irow)) {
204 Int iIF = cDataDescIdCol(irow);
205 haveXPol(iIF) = True;
206 }
207 }
208 }
209
210
211 // Initialize member data.
212 cBeams.assign(beams);
213 cIFs.assign(IFs);
214 cNChan.assign(nChan);
215 cNPol.assign(nPol);
216 cHaveXPol.assign(haveXPol);
217
218
219 // Default channel range selection.
220 cStartChan.resize(nIF);
221 cEndChan.resize(nIF);
222 cRefChan.resize(nIF);
223
224 for (uInt iIF = 0; iIF < nIF; iIF++) {
225 cStartChan(iIF) = 1;
226 cEndChan(iIF) = cNChan(iIF);
227 cRefChan(iIF) = cNChan(iIF)/2 + 1;
228 }
229
230 Slice all;
231 cDataSel.resize(nIF);
232 cDataSel = Slicer(all, all);
233
234 cScanNo = 0;
235 cCycleNo = 1;
236 cTime = cTimeCol(0);
237
238 return 0;
239}
240
241//---------------------------------------------------- PKSMS2reader::getHeader
242
243// Get parameters describing the data.
244
245Int PKSMS2reader::getHeader(
246 String &observer,
247 String &project,
248 String &antName,
249 Vector<Double> &antPosition,
250 String &obsMode,
[1399]251 String &bunit,
[1325]252 Float &equinox,
253 String &dopplerFrame,
254 Double &mjd,
255 Double &refFreq,
256 Double &bandwidth)
257{
258 if (!cMSopen) {
259 return 1;
260 }
261
262 // Observer and project.
263 ROMSObservationColumns observationCols(cPKSMS.observation());
264 observer = observationCols.observer()(0);
265 project = observationCols.project()(0);
266
267 // Antenna name and ITRF coordinates.
268 ROMSAntennaColumns antennaCols(cPKSMS.antenna());
269 antName = antennaCols.name()(0);
270 antPosition = antennaCols.position()(0);
271
272 // Observation type.
273 if (cObsModeCol.nrow()) {
274 obsMode = cObsModeCol(0);
275 if (obsMode == "\0") obsMode = "RF";
276 } else {
277 obsMode = "RF";
278 }
279
[1399]280 // Brightness units.
281 bunit = cPKSMS.unit(MSMainEnums::FLOAT_DATA);
[1325]282
283 // Coordinate equinox.
284 ROMSPointingColumns pointingCols(cPKSMS.pointing());
285 String dirref = pointingCols.direction().keywordSet().asRecord("MEASINFO").
286 asString("Ref");
287 sscanf(dirref.chars()+1, "%f", &equinox);
288
289 // Frequency/velocity reference frame.
290 ROMSSpWindowColumns spWinCols(cPKSMS.spectralWindow());
291 dopplerFrame = MFrequency::showType(spWinCols.measFreqRef()(0));
292
293 // Translate to FITS standard names.
294 if (dopplerFrame == "TOPO") {
295 dopplerFrame = "TOPOCENT";
296 } else if (dopplerFrame == "GEO") {
297 dopplerFrame = "GEOCENTR";
298 } else if (dopplerFrame == "BARY") {
299 dopplerFrame = "BARYCENT";
300 } else if (dopplerFrame == "GALACTO") {
301 dopplerFrame = "GALACTOC";
302 } else if (dopplerFrame == "LGROUP") {
303 dopplerFrame = "LOCALGRP";
304 } else if (dopplerFrame == "CMB") {
305 dopplerFrame = "CMBDIPOL";
306 } else if (dopplerFrame == "REST") {
307 dopplerFrame = "SOURCE";
308 }
309
310 // MJD at start of observation.
311 mjd = cTimeCol(0)/86400.0;
312
313 // Reference frequency and bandwidth.
314 refFreq = spWinCols.refFrequency()(0);
315 bandwidth = spWinCols.totalBandwidth()(0);
316
317 return 0;
318}
319
320//-------------------------------------------------- PKSMS2reader::getFreqInfo
321
322// Get frequency parameters for each IF.
323
324Int PKSMS2reader::getFreqInfo(
325 Vector<Double> &startFreq,
326 Vector<Double> &endFreq)
327{
328 uInt nIF = cIFs.nelements();
329 startFreq.resize(nIF);
330 endFreq.resize(nIF);
331
332 for (uInt iIF = 0; iIF < nIF; iIF++) {
333 Vector<Double> chanFreq = cChanFreqCol(iIF);
334
335 Int nChan = chanFreq.nelements();
336 startFreq(iIF) = chanFreq(0);
337 endFreq(iIF) = chanFreq(nChan-1);
338 }
339
340 return 0;
341}
342
343//------------------------------------------------------- PKSMS2reader::select
344
345// Set data selection by beam number and channel.
346
347uInt PKSMS2reader::select(
348 const Vector<Bool> beamSel,
349 const Vector<Bool> IFsel,
350 const Vector<Int> startChan,
351 const Vector<Int> endChan,
352 const Vector<Int> refChan,
353 const Bool getSpectra,
354 const Bool getXPol,
355 const Bool getFeedPos)
356{
357 if (!cMSopen) {
358 return 1;
359 }
360
361 // Beam selection.
362 uInt nBeam = cBeams.nelements();
363 uInt nBeamSel = beamSel.nelements();
364 for (uInt ibeam = 0; ibeam < nBeam; ibeam++) {
365 if (ibeam < nBeamSel) {
366 cBeams(ibeam) = beamSel(ibeam);
367 } else {
368 cBeams(ibeam) = False;
369 }
370 }
371
372 uInt nIF = cIFs.nelements();
373 uInt maxNChan = 0;
374 for (uInt iIF = 0; iIF < nIF; iIF++) {
375 // IF selection.
376 if (iIF < IFsel.nelements()) {
377 cIFs(iIF) = IFsel(iIF);
378 } else {
379 cIFs(iIF) = False;
380 }
381
382 if (!cIFs(iIF)) continue;
383
384
385 // Channel selection.
386 if (iIF < startChan.nelements()) {
387 cStartChan(iIF) = startChan(iIF);
388
389 if (cStartChan(iIF) <= 0) {
390 cStartChan(iIF) += cNChan(iIF);
391 } else if (cStartChan(iIF) > Int(cNChan(iIF))) {
392 cStartChan(iIF) = cNChan(iIF);
393 }
394 }
395
396 if (iIF < endChan.nelements()) {
397 cEndChan(iIF) = endChan(iIF);
398
399 if (cEndChan(iIF) <= 0) {
400 cEndChan(iIF) += cNChan(iIF);
401 } else if (cEndChan(iIF) > Int(cNChan(iIF))) {
402 cEndChan(iIF) = cNChan(iIF);
403 }
404 }
405
406 if (iIF < refChan.nelements()) {
407 cRefChan(iIF) = refChan(iIF);
408 } else {
409 cRefChan(iIF) = cStartChan(iIF);
410 if (cStartChan(iIF) <= cEndChan(iIF)) {
411 cRefChan(iIF) += (cEndChan(iIF) - cStartChan(iIF) + 1)/2;
412 } else {
413 cRefChan(iIF) -= (cStartChan(iIF) - cEndChan(iIF) + 1)/2;
414 }
415 }
416
417 uInt nChan = abs(cEndChan(iIF) - cStartChan(iIF)) + 1;
418 if (maxNChan < nChan) {
419 maxNChan = nChan;
420 }
421
422 // Inverted Slices are not allowed.
423 Slice outPols;
424 Slice outChans(min(cStartChan(iIF),cEndChan(iIF))-1, nChan);
425 cDataSel(iIF) = Slicer(outPols, outChans);
426 }
427
428 // Get spectral data?
429 cGetSpectra = getSpectra;
430
431 // Get cross-polarization data?
432 cGetXPol = cGetXPol && getXPol;
433
434 // Get feed positions? (Not available.)
435 cGetFeedPos = False;
436
437 return maxNChan;
438}
439
440//---------------------------------------------------- PKSMS2reader::findRange
441
442// Find the range of the data in time and position.
443
444Int PKSMS2reader::findRange(
445 Int &nRow,
446 Int &nSel,
447 Vector<Double> &timeSpan,
448 Matrix<Double> &positions)
449{
450 if (!cMSopen) {
451 return 1;
452 }
453
454 nRow = cNRow;
455
456 // Find the number of rows selected.
457 nSel = 0;
458 Vector<Bool> sel(nRow);
459 for (Int irow = 0; irow < nRow; irow++) {
460 if ((sel(irow) = cBeams(cBeamNoCol(irow)) &&
461 cIFs(cDataDescIdCol(irow)))) {
462 nSel++;
463 }
464 }
465
466 // Find the time range (s).
467 timeSpan.resize(2);
468 timeSpan(0) = cTimeCol(0);
469 timeSpan(1) = cTimeCol(nRow-1);
470
471 // Retrieve positions for selected data.
472 Int isel = 0;
473 positions.resize(2,nSel);
474 for (Int irow = 0; irow < nRow; irow++) {
475 if (sel(irow)) {
476 Matrix<Double> pointingDir = cPointingCol(cFieldIdCol(irow));
477 positions.column(isel++) = pointingDir.column(0);
478 }
479 }
480
481 return 0;
482}
483
484//--------------------------------------------------------- PKSMS2reader::read
485
486// Read the next data record.
487
488Int PKSMS2reader::read(
489 Int &scanNo,
490 Int &cycleNo,
491 Double &mjd,
492 Double &interval,
493 String &fieldName,
494 String &srcName,
495 Vector<Double> &srcDir,
496 Vector<Double> &srcPM,
497 Double &srcVel,
498 String &obsMode,
499 Int &IFno,
500 Double &refFreq,
501 Double &bandwidth,
502 Double &freqInc,
503 Double &restFreq,
504 Vector<Float> &tcal,
505 String &tcalTime,
506 Float &azimuth,
507 Float &elevation,
508 Float &parAngle,
509 Float &focusAxi,
510 Float &focusTan,
511 Float &focusRot,
512 Float &temperature,
513 Float &pressure,
514 Float &humidity,
515 Float &windSpeed,
516 Float &windAz,
517 Int &refBeam,
518 Int &beamNo,
519 Vector<Double> &direction,
520 Vector<Double> &scanRate,
521 Vector<Float> &tsys,
522 Vector<Float> &sigma,
523 Vector<Float> &calFctr,
524 Matrix<Float> &baseLin,
525 Matrix<Float> &baseSub,
526 Matrix<Float> &spectra,
527 Matrix<uChar> &flagged,
528 Complex &xCalFctr,
529 Vector<Complex> &xPol)
530{
531 if (!cMSopen) {
532 return 1;
533 }
534
535 // Check for EOF.
536 if (cIdx >= cNRow) {
537 return -1;
538 }
539
540 // Find the next selected beam and IF.
541 Int ibeam;
542 Int iIF;
543 while (True) {
544 ibeam = cBeamNoCol(cIdx);
545 iIF = cDataDescIdCol(cIdx);
546 if (cBeams(ibeam) && cIFs(iIF)) {
547 break;
548 }
549
550 // Check for EOF.
551 if (++cIdx >= cNRow) {
552 return -1;
553 }
554 }
555
556 // Renumerate scan no. Here still is 1-based
557 scanNo = cScanNoCol(cIdx) - cScanNoCol(0) + 1;
558
559 if (scanNo != cScanNo) {
560 // Start of new scan.
561 cScanNo = scanNo;
562 cCycleNo = 1;
563 cTime = cTimeCol(cIdx);
564 }
565
566 Double time = cTimeCol(cIdx);
567 mjd = time/86400.0;
568 interval = cIntervalCol(cIdx);
569
570 // Reconstruct the integration cycle number; due to small latencies the
571 // integration time is usually slightly less than the time between cycles,
572 // resetting cTime will prevent the difference from accumulating.
573 cCycleNo += nint((time - cTime)/interval);
574 cycleNo = cCycleNo;
575 cTime = time;
576
577 Int fieldId = cFieldIdCol(cIdx);
578 fieldName = cFieldNameCol(fieldId);
579
580 Int srcId = cSrcIdCol(fieldId);
581 srcName = cSrcNameCol(srcId);
582 srcDir = cSrcDirCol(srcId);
583 srcPM = cSrcPMCol(srcId);
584
585 // Systemic velocity.
586 if (!cHaveSrcVel) {
587 srcVel = 0.0f;
588 } else {
589 srcVel = cSrcVelCol(srcId)(IPosition(1,0));
590 }
591
592 // Observation type.
593 Int stateId = cStateIdCol(cIdx);
594 obsMode = cObsModeCol(stateId);
595
596 IFno = iIF + 1;
597 Int nChan = abs(cEndChan(iIF) - cStartChan(iIF)) + 1;
598
599 // Minimal handling on continuum data.
600 Vector<Double> chanFreq = cChanFreqCol(iIF);
601 if (nChan == 1) {
602 cout << "The input is continuum data. "<< endl;
603 freqInc = chanFreq(0);
604 refFreq = chanFreq(0);
605 restFreq = 0.0f;
606 } else {
607 if (cStartChan(iIF) <= cEndChan(iIF)) {
608 freqInc = chanFreq(1) - chanFreq(0);
609 } else {
610 freqInc = chanFreq(0) - chanFreq(1);
611 }
612
613 refFreq = chanFreq(cRefChan(iIF)-1);
614 restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0));
615 }
616 bandwidth = abs(freqInc * nChan);
617
618 tcal.resize(cNPol(iIF));
619 tcal = 0.0f;
620 tcalTime = "";
621 azimuth = 0.0f;
622 elevation = 0.0f;
623 parAngle = 0.0f;
624 focusAxi = 0.0f;
625 focusTan = 0.0f;
626 focusRot = 0.0f;
627
628 // Find the appropriate entry in the WEATHER subtable.
629 Vector<Double> wTimes = cWeatherTimeCol.getColumn();
630 Int weatherIdx;
631 for (weatherIdx = wTimes.nelements()-1; weatherIdx >= 0; weatherIdx--) {
632 if (cWeatherTimeCol(weatherIdx) <= time) {
633 break;
634 }
635 }
636
637 if (weatherIdx < 0) {
638 // No appropriate WEATHER entry.
639 pressure = 0.0f;
640 humidity = 0.0f;
641 temperature = 0.0f;
642 } else {
643 pressure = cPressureCol(weatherIdx);
644 humidity = cHumidityCol(weatherIdx);
645 temperature = cTemperatureCol(weatherIdx);
646 }
647
648 windSpeed = 0.0f;
649 windAz = 0.0f;
650
651 refBeam = 0;
652 beamNo = ibeam + 1;
653
654 Matrix<Double> pointingDir = cPointingCol(fieldId);
655 direction = pointingDir.column(0);
656 uInt ncols = pointingDir.ncolumn();
657 if (ncols == 1) {
658 scanRate = 0.0f;
659 } else {
660 scanRate = pointingDir.column(1);
661 }
662
663 // Get Tsys assuming that entries in the SYSCAL table match the main table.
664 if (cHaveTsys) {
665 Int nTsysColRow = cTsysCol.nrow();
666 if (nTsysColRow != cNRow) {
667 cHaveTsys=0;
668 }
669 }
670 if (cHaveTsys) {
671 cTsysCol.get(cIdx, tsys, True);
672 } else {
673 Int numReceptor;
674 cNumReceptorCol.get(0, numReceptor);
675 tsys.resize(numReceptor);
676 tsys = 1.0f;
677 }
678 cSigmaCol.get(cIdx, sigma, True);
679
680 // Calibration factors (if available).
681 calFctr.resize(cNPol(iIF));
682 if (cHaveCalFctr) {
683 cCalFctrCol.get(cIdx, calFctr);
684 } else {
685 calFctr = 0.0f;
686 }
687
688 // Baseline parameters (if available).
689 if (cHaveBaseLin) {
690 baseLin.resize(2,cNPol(iIF));
691 cBaseLinCol.get(cIdx, baseLin);
692
693 baseSub.resize(9,cNPol(iIF));
694 cBaseSubCol.get(cIdx, baseSub);
695
696 } else {
697 baseLin.resize(0,0);
698 baseSub.resize(0,0);
699 }
700
701
702 // Get spectral data.
703 if (cGetSpectra) {
704 Matrix<Float> tmpData;
705 Matrix<Bool> tmpFlag;
706 cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
707 cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
708
709 // Transpose spectra.
710 Int nPol = tmpData.nrow();
711 spectra.resize(nChan, nPol);
712 flagged.resize(nChan, nPol);
713 if (cEndChan(iIF) >= cStartChan(iIF)) {
714 // Simple transposition.
715 for (Int ipol = 0; ipol < nPol; ipol++) {
716 for (Int ichan = 0; ichan < nChan; ichan++) {
717 spectra(ichan,ipol) = tmpData(ipol,ichan);
718 flagged(ichan,ipol) = tmpFlag(ipol,ichan);
719 }
720 }
721
722 } else {
723 // Transpose with inversion.
724 Int jchan = nChan - 1;
725 for (Int ipol = 0; ipol < nPol; ipol++) {
726 for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
727 spectra(ichan,ipol) = tmpData(ipol,jchan);
728 flagged(ichan,ipol) = tmpFlag(ipol,jchan);
729 }
730 }
731 }
732 }
733
734 // Get cross-polarization data.
735 if (cGetXPol) {
736 if (cHaveXCalFctr) {
737 cXCalFctrCol.get(cIdx, xCalFctr);
738 } else {
739 xCalFctr = Complex(0.0f, 0.0f);
740 }
741
742 cDataCol.get(cIdx, xPol, True);
743
744 if (cEndChan(iIF) < cStartChan(iIF)) {
745 Complex ctmp;
746 Int jchan = nChan - 1;
747 for (Int ichan = 0; ichan < nChan/2; ichan++, jchan--) {
748 ctmp = xPol(ichan);
749 xPol(ichan) = xPol(jchan);
750 xPol(jchan) = ctmp;
751 }
752 }
753 }
754
755 cIdx++;
756
757 return 0;
758}
759
760//--------------------------------------------------------- PKSMS2reader::read
761
762// Read the next data record, just the basics.
763
764Int PKSMS2reader::read(
765 Int &IFno,
766 Vector<Float> &tsys,
767 Vector<Float> &calFctr,
768 Matrix<Float> &baseLin,
769 Matrix<Float> &baseSub,
770 Matrix<Float> &spectra,
771 Matrix<uChar> &flagged)
772{
773 if (!cMSopen) {
774 return 1;
775 }
776
777 // Check for EOF.
778 if (cIdx >= cNRow) {
779 return -1;
780 }
781
782 // Find the next selected beam and IF.
783 Int ibeam;
784 Int iIF;
785 while (True) {
786 ibeam = cBeamNoCol(cIdx);
787 iIF = cDataDescIdCol(cIdx);
788 if (cBeams(ibeam) && cIFs(iIF)) {
789 break;
790 }
791
792 // Check for EOF.
793 if (++cIdx >= cNRow) {
794 return -1;
795 }
796 }
797
798 IFno = iIF + 1;
799
800 // Get Tsys assuming that entries in the SYSCAL table match the main table.
801 cTsysCol.get(cIdx, tsys, True);
802
803 // Calibration factors (if available).
804 if (cHaveCalFctr) {
805 cCalFctrCol.get(cIdx, calFctr, True);
806 } else {
807 calFctr.resize(cNPol(iIF));
808 calFctr = 0.0f;
809 }
810
811 // Baseline parameters (if available).
812 if (cHaveBaseLin) {
813 baseLin.resize(2,cNPol(iIF));
814 cBaseLinCol.get(cIdx, baseLin);
815
816 baseSub.resize(9,cNPol(iIF));
817 cBaseSubCol.get(cIdx, baseSub);
818
819 } else {
820 baseLin.resize(0,0);
821 baseSub.resize(0,0);
822 }
823
824 if (cGetSpectra) {
825 // Get spectral data.
826 Matrix<Float> tmpData;
827 Matrix<Bool> tmpFlag;
828 cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
829 cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
830
831 // Transpose spectra.
832 Int nChan = tmpData.ncolumn();
833 Int nPol = tmpData.nrow();
834 spectra.resize(nChan, nPol);
835 flagged.resize(nChan, nPol);
836 if (cEndChan(iIF) >= cStartChan(iIF)) {
837 // Simple transposition.
838 for (Int ipol = 0; ipol < nPol; ipol++) {
839 for (Int ichan = 0; ichan < nChan; ichan++) {
840 spectra(ichan,ipol) = tmpData(ipol,ichan);
841 flagged(ichan,ipol) = tmpFlag(ipol,ichan);
842 }
843 }
844
845 } else {
846 // Transpose with inversion.
847 Int jchan = nChan - 1;
848 for (Int ipol = 0; ipol < nPol; ipol++) {
849 for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
850 spectra(ichan,ipol) = tmpData(ipol,jchan);
851 flagged(ichan,ipol) = tmpFlag(ipol,jchan);
852 }
853 }
854 }
855 }
856
857 cIdx++;
858
859 return 0;
860}
861
862//-------------------------------------------------------- PKSMS2reader::close
863
864// Close the MS.
865
866void PKSMS2reader::close()
867{
868 cPKSMS = MeasurementSet();
869 cMSopen = False;
870}
Note: See TracBrowser for help on using the repository browser.