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

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

sync with livedata/implement/atnf
File size: 21.4 KB
Line 
1//#---------------------------------------------------------------------------
2//# PKSMS2reader.cc: Class to read Parkes Multibeam data from a v2 MS.
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2000-2008
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: PKSMS2reader.cc,v 19.13 2008-06-26 02:08:02 cal103 Exp $
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 &obsType,
251 String &bunit,
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 obsType = cObsModeCol(0);
275 if (obsType == "\0") obsType = "RF";
276 } else {
277 obsType = "RF";
278 }
279
280 // Brightness units.
281 bunit = cPKSMS.unit(MSMainEnums::FLOAT_DATA);
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(MBrecord &MBrec)
489{
490 if (!cMSopen) {
491 return 1;
492 }
493
494 // Check for EOF.
495 if (cIdx >= cNRow) {
496 return -1;
497 }
498
499 // Find the next selected beam and IF.
500 Int ibeam;
501 Int iIF;
502 while (True) {
503 ibeam = cBeamNoCol(cIdx);
504 iIF = cDataDescIdCol(cIdx);
505 if (cBeams(ibeam) && cIFs(iIF)) {
506 break;
507 }
508
509 // Check for EOF.
510 if (++cIdx >= cNRow) {
511 return -1;
512 }
513 }
514
515 // Renumerate scan no. Here still is 1-based
516 MBrec.scanNo = cScanNoCol(cIdx) - cScanNoCol(0) + 1;
517
518 if (MBrec.scanNo != cScanNo) {
519 // Start of new scan.
520 cScanNo = MBrec.scanNo;
521 cCycleNo = 1;
522 cTime = cTimeCol(cIdx);
523 }
524
525 Double time = cTimeCol(cIdx);
526 MBrec.mjd = time/86400.0;
527 MBrec.interval = cIntervalCol(cIdx);
528
529 // Reconstruct the integration cycle number; due to small latencies the
530 // integration time is usually slightly less than the time between cycles,
531 // resetting cTime will prevent the difference from accumulating.
532 cCycleNo += nint((time - cTime)/MBrec.interval);
533 MBrec.cycleNo = cCycleNo;
534 cTime = time;
535
536 Int fieldId = cFieldIdCol(cIdx);
537 MBrec.fieldName = cFieldNameCol(fieldId);
538
539 Int srcId = cSrcIdCol(fieldId);
540 MBrec.srcName = cSrcNameCol(srcId);
541 MBrec.srcDir = cSrcDirCol(srcId);
542 MBrec.srcPM = cSrcPMCol(srcId);
543
544 // Systemic velocity.
545 if (!cHaveSrcVel) {
546 MBrec.srcVel = 0.0f;
547 } else {
548 MBrec.srcVel = cSrcVelCol(srcId)(IPosition(1,0));
549 }
550
551 // Observation type.
552 Int stateId = cStateIdCol(cIdx);
553 MBrec.obsType = cObsModeCol(stateId);
554
555 MBrec.IFno = iIF + 1;
556 Int nChan = abs(cEndChan(iIF) - cStartChan(iIF)) + 1;
557
558 // Minimal handling on continuum data.
559 Vector<Double> chanFreq = cChanFreqCol(iIF);
560 if (nChan == 1) {
561 cout << "The input is continuum data. "<< endl;
562 MBrec.freqInc = chanFreq(0);
563 MBrec.refFreq = chanFreq(0);
564 MBrec.restFreq = 0.0f;
565 } else {
566 if (cStartChan(iIF) <= cEndChan(iIF)) {
567 MBrec.freqInc = chanFreq(1) - chanFreq(0);
568 } else {
569 MBrec.freqInc = chanFreq(0) - chanFreq(1);
570 }
571
572 MBrec.refFreq = chanFreq(cRefChan(iIF)-1);
573 MBrec.restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0));
574 }
575 MBrec.bandwidth = abs(MBrec.freqInc * nChan);
576
577 MBrec.tcal.resize(cNPol(iIF));
578 MBrec.tcal = 0.0f;
579 MBrec.tcalTime = "";
580 MBrec.azimuth = 0.0f;
581 MBrec.elevation = 0.0f;
582 MBrec.parAngle = 0.0f;
583 MBrec.focusAxi = 0.0f;
584 MBrec.focusTan = 0.0f;
585 MBrec.focusRot = 0.0f;
586
587 // Find the appropriate entry in the WEATHER subtable.
588 Vector<Double> wTimes = cWeatherTimeCol.getColumn();
589 Int weatherIdx;
590 for (weatherIdx = wTimes.nelements()-1; weatherIdx >= 0; weatherIdx--) {
591 if (cWeatherTimeCol(weatherIdx) <= time) {
592 break;
593 }
594 }
595
596 if (weatherIdx < 0) {
597 // No appropriate WEATHER entry.
598 MBrec.temperature = 0.0f;
599 MBrec.pressure = 0.0f;
600 MBrec.humidity = 0.0f;
601 } else {
602 MBrec.temperature = cTemperatureCol(weatherIdx);
603 MBrec.pressure = cPressureCol(weatherIdx);
604 MBrec.humidity = cHumidityCol(weatherIdx);
605 }
606
607 MBrec.windSpeed = 0.0f;
608 MBrec.windAz = 0.0f;
609
610 MBrec.refBeam = 0;
611 MBrec.beamNo = ibeam + 1;
612
613 Matrix<Double> pointingDir = cPointingCol(fieldId);
614 MBrec.direction = pointingDir.column(0);
615 uInt ncols = pointingDir.ncolumn();
616 if (ncols == 1) {
617 MBrec.scanRate = 0.0f;
618 } else {
619 MBrec.scanRate = pointingDir.column(1);
620 }
621 MBrec.rateAge = 0;
622 MBrec.rateson = 0;
623
624 // Get Tsys assuming that entries in the SYSCAL table match the main table.
625 if (cHaveTsys) {
626 Int nTsysColRow = cTsysCol.nrow();
627 if (nTsysColRow != cNRow) {
628 cHaveTsys=0;
629 }
630 }
631 if (cHaveTsys) {
632 cTsysCol.get(cIdx, MBrec.tsys, True);
633 } else {
634 Int numReceptor;
635 cNumReceptorCol.get(0, numReceptor);
636 MBrec.tsys.resize(numReceptor);
637 MBrec.tsys = 1.0f;
638 }
639 cSigmaCol.get(cIdx, MBrec.sigma, True);
640
641 // Calibration factors (if available).
642 MBrec.calFctr.resize(cNPol(iIF));
643 if (cHaveCalFctr) {
644 cCalFctrCol.get(cIdx, MBrec.calFctr);
645 } else {
646 MBrec.calFctr = 0.0f;
647 }
648
649 // Baseline parameters (if available).
650 if (cHaveBaseLin) {
651 MBrec.baseLin.resize(2,cNPol(iIF));
652 cBaseLinCol.get(cIdx, MBrec.baseLin);
653
654 MBrec.baseSub.resize(9,cNPol(iIF));
655 cBaseSubCol.get(cIdx, MBrec.baseSub);
656
657 } else {
658 MBrec.baseLin.resize(0,0);
659 MBrec.baseSub.resize(0,0);
660 }
661
662
663 // Get spectral data.
664 if (cGetSpectra) {
665 Matrix<Float> tmpData;
666 Matrix<Bool> tmpFlag;
667 cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
668 cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
669
670 // Transpose spectra.
671 Int nPol = tmpData.nrow();
672 MBrec.spectra.resize(nChan, nPol);
673 MBrec.flagged.resize(nChan, nPol);
674 if (cEndChan(iIF) >= cStartChan(iIF)) {
675 // Simple transposition.
676 for (Int ipol = 0; ipol < nPol; ipol++) {
677 for (Int ichan = 0; ichan < nChan; ichan++) {
678 MBrec.spectra(ichan,ipol) = tmpData(ipol,ichan);
679 MBrec.flagged(ichan,ipol) = tmpFlag(ipol,ichan);
680 }
681 }
682
683 } else {
684 // Transpose with inversion.
685 Int jchan = nChan - 1;
686 for (Int ipol = 0; ipol < nPol; ipol++) {
687 for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
688 MBrec.spectra(ichan,ipol) = tmpData(ipol,jchan);
689 MBrec.flagged(ichan,ipol) = tmpFlag(ipol,jchan);
690 }
691 }
692 }
693 }
694
695 // Get cross-polarization data.
696 if (cGetXPol) {
697 if (cHaveXCalFctr) {
698 cXCalFctrCol.get(cIdx, MBrec.xCalFctr);
699 } else {
700 MBrec.xCalFctr = Complex(0.0f, 0.0f);
701 }
702
703 cDataCol.get(cIdx, MBrec.xPol, True);
704
705 if (cEndChan(iIF) < cStartChan(iIF)) {
706 Complex ctmp;
707 Int jchan = nChan - 1;
708 for (Int ichan = 0; ichan < nChan/2; ichan++, jchan--) {
709 ctmp = MBrec.xPol(ichan);
710 MBrec.xPol(ichan) = MBrec.xPol(jchan);
711 MBrec.xPol(jchan) = ctmp;
712 }
713 }
714 }
715
716 cIdx++;
717
718 return 0;
719}
720
721//--------------------------------------------------------- PKSMS2reader::read
722
723// Read the next data record, just the basics.
724
725Int PKSMS2reader::read(
726 Int &IFno,
727 Vector<Float> &tsys,
728 Vector<Float> &calFctr,
729 Matrix<Float> &baseLin,
730 Matrix<Float> &baseSub,
731 Matrix<Float> &spectra,
732 Matrix<uChar> &flagged)
733{
734 if (!cMSopen) {
735 return 1;
736 }
737
738 // Check for EOF.
739 if (cIdx >= cNRow) {
740 return -1;
741 }
742
743 // Find the next selected beam and IF.
744 Int ibeam;
745 Int iIF;
746 while (True) {
747 ibeam = cBeamNoCol(cIdx);
748 iIF = cDataDescIdCol(cIdx);
749 if (cBeams(ibeam) && cIFs(iIF)) {
750 break;
751 }
752
753 // Check for EOF.
754 if (++cIdx >= cNRow) {
755 return -1;
756 }
757 }
758
759 IFno = iIF + 1;
760
761 // Get Tsys assuming that entries in the SYSCAL table match the main table.
762 cTsysCol.get(cIdx, tsys, True);
763
764 // Calibration factors (if available).
765 if (cHaveCalFctr) {
766 cCalFctrCol.get(cIdx, calFctr, True);
767 } else {
768 calFctr.resize(cNPol(iIF));
769 calFctr = 0.0f;
770 }
771
772 // Baseline parameters (if available).
773 if (cHaveBaseLin) {
774 baseLin.resize(2,cNPol(iIF));
775 cBaseLinCol.get(cIdx, baseLin);
776
777 baseSub.resize(9,cNPol(iIF));
778 cBaseSubCol.get(cIdx, baseSub);
779
780 } else {
781 baseLin.resize(0,0);
782 baseSub.resize(0,0);
783 }
784
785 if (cGetSpectra) {
786 // Get spectral data.
787 Matrix<Float> tmpData;
788 Matrix<Bool> tmpFlag;
789 cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
790 cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
791
792 // Transpose spectra.
793 Int nChan = tmpData.ncolumn();
794 Int nPol = tmpData.nrow();
795 spectra.resize(nChan, nPol);
796 flagged.resize(nChan, nPol);
797 if (cEndChan(iIF) >= cStartChan(iIF)) {
798 // Simple transposition.
799 for (Int ipol = 0; ipol < nPol; ipol++) {
800 for (Int ichan = 0; ichan < nChan; ichan++) {
801 spectra(ichan,ipol) = tmpData(ipol,ichan);
802 flagged(ichan,ipol) = tmpFlag(ipol,ichan);
803 }
804 }
805
806 } else {
807 // Transpose with inversion.
808 Int jchan = nChan - 1;
809 for (Int ipol = 0; ipol < nPol; ipol++) {
810 for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
811 spectra(ichan,ipol) = tmpData(ipol,jchan);
812 flagged(ichan,ipol) = tmpFlag(ipol,jchan);
813 }
814 }
815 }
816 }
817
818 cIdx++;
819
820 return 0;
821}
822
823//-------------------------------------------------------- PKSMS2reader::close
824
825// Close the MS.
826
827void PKSMS2reader::close()
828{
829 cPKSMS = MeasurementSet();
830 cMSopen = False;
831}
Note: See TracBrowser for help on using the repository browser.