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

Last change on this file since 2909 was 1720, checked in by Malte Marquarding, 15 years ago

Update from livedata CVS repository

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