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

Last change on this file since 1375 was 1325, checked in by mar637, 18 years ago

Changes to use casacore instead of casa_asap/aips++\nAdded atnf PKSIO library snapshot to external and linking against this local copy
File size: 22.2 KB
RevLine 
[1325]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: PKSMS2reader.cc,v 19.11 2006/07/05 04:59:20 mcalabre 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 &obsMode,
251 Float &equinox,
252 String &dopplerFrame,
253 Double &mjd,
254 Double &refFreq,
255 Double &bandwidth)
256{
257 if (!cMSopen) {
258 return 1;
259 }
260
261 // Observer and project.
262 ROMSObservationColumns observationCols(cPKSMS.observation());
263 observer = observationCols.observer()(0);
264 project = observationCols.project()(0);
265
266 // Antenna name and ITRF coordinates.
267 ROMSAntennaColumns antennaCols(cPKSMS.antenna());
268 antName = antennaCols.name()(0);
269 antPosition = antennaCols.position()(0);
270
271 // Observation type.
272 if (cObsModeCol.nrow()) {
273 obsMode = cObsModeCol(0);
274 if (obsMode == "\0") obsMode = "RF";
275 } else {
276 obsMode = "RF";
277 }
278
279
280 // Coordinate equinox.
281 ROMSPointingColumns pointingCols(cPKSMS.pointing());
282 String dirref = pointingCols.direction().keywordSet().asRecord("MEASINFO").
283 asString("Ref");
284 sscanf(dirref.chars()+1, "%f", &equinox);
285
286 // Frequency/velocity reference frame.
287 ROMSSpWindowColumns spWinCols(cPKSMS.spectralWindow());
288 dopplerFrame = MFrequency::showType(spWinCols.measFreqRef()(0));
289
290 // Translate to FITS standard names.
291 if (dopplerFrame == "TOPO") {
292 dopplerFrame = "TOPOCENT";
293 } else if (dopplerFrame == "GEO") {
294 dopplerFrame = "GEOCENTR";
295 } else if (dopplerFrame == "BARY") {
296 dopplerFrame = "BARYCENT";
297 } else if (dopplerFrame == "GALACTO") {
298 dopplerFrame = "GALACTOC";
299 } else if (dopplerFrame == "LGROUP") {
300 dopplerFrame = "LOCALGRP";
301 } else if (dopplerFrame == "CMB") {
302 dopplerFrame = "CMBDIPOL";
303 } else if (dopplerFrame == "REST") {
304 dopplerFrame = "SOURCE";
305 }
306
307 // MJD at start of observation.
308 mjd = cTimeCol(0)/86400.0;
309
310 // Reference frequency and bandwidth.
311 refFreq = spWinCols.refFrequency()(0);
312 bandwidth = spWinCols.totalBandwidth()(0);
313
314 return 0;
315}
316
317//-------------------------------------------------- PKSMS2reader::getFreqInfo
318
319// Get frequency parameters for each IF.
320
321Int PKSMS2reader::getFreqInfo(
322 Vector<Double> &startFreq,
323 Vector<Double> &endFreq)
324{
325 uInt nIF = cIFs.nelements();
326 startFreq.resize(nIF);
327 endFreq.resize(nIF);
328
329 for (uInt iIF = 0; iIF < nIF; iIF++) {
330 Vector<Double> chanFreq = cChanFreqCol(iIF);
331
332 Int nChan = chanFreq.nelements();
333 startFreq(iIF) = chanFreq(0);
334 endFreq(iIF) = chanFreq(nChan-1);
335 }
336
337 return 0;
338}
339
340//------------------------------------------------------- PKSMS2reader::select
341
342// Set data selection by beam number and channel.
343
344uInt PKSMS2reader::select(
345 const Vector<Bool> beamSel,
346 const Vector<Bool> IFsel,
347 const Vector<Int> startChan,
348 const Vector<Int> endChan,
349 const Vector<Int> refChan,
350 const Bool getSpectra,
351 const Bool getXPol,
352 const Bool getFeedPos)
353{
354 if (!cMSopen) {
355 return 1;
356 }
357
358 // Beam selection.
359 uInt nBeam = cBeams.nelements();
360 uInt nBeamSel = beamSel.nelements();
361 for (uInt ibeam = 0; ibeam < nBeam; ibeam++) {
362 if (ibeam < nBeamSel) {
363 cBeams(ibeam) = beamSel(ibeam);
364 } else {
365 cBeams(ibeam) = False;
366 }
367 }
368
369 uInt nIF = cIFs.nelements();
370 uInt maxNChan = 0;
371 for (uInt iIF = 0; iIF < nIF; iIF++) {
372 // IF selection.
373 if (iIF < IFsel.nelements()) {
374 cIFs(iIF) = IFsel(iIF);
375 } else {
376 cIFs(iIF) = False;
377 }
378
379 if (!cIFs(iIF)) continue;
380
381
382 // Channel selection.
383 if (iIF < startChan.nelements()) {
384 cStartChan(iIF) = startChan(iIF);
385
386 if (cStartChan(iIF) <= 0) {
387 cStartChan(iIF) += cNChan(iIF);
388 } else if (cStartChan(iIF) > Int(cNChan(iIF))) {
389 cStartChan(iIF) = cNChan(iIF);
390 }
391 }
392
393 if (iIF < endChan.nelements()) {
394 cEndChan(iIF) = endChan(iIF);
395
396 if (cEndChan(iIF) <= 0) {
397 cEndChan(iIF) += cNChan(iIF);
398 } else if (cEndChan(iIF) > Int(cNChan(iIF))) {
399 cEndChan(iIF) = cNChan(iIF);
400 }
401 }
402
403 if (iIF < refChan.nelements()) {
404 cRefChan(iIF) = refChan(iIF);
405 } else {
406 cRefChan(iIF) = cStartChan(iIF);
407 if (cStartChan(iIF) <= cEndChan(iIF)) {
408 cRefChan(iIF) += (cEndChan(iIF) - cStartChan(iIF) + 1)/2;
409 } else {
410 cRefChan(iIF) -= (cStartChan(iIF) - cEndChan(iIF) + 1)/2;
411 }
412 }
413
414 uInt nChan = abs(cEndChan(iIF) - cStartChan(iIF)) + 1;
415 if (maxNChan < nChan) {
416 maxNChan = nChan;
417 }
418
419 // Inverted Slices are not allowed.
420 Slice outPols;
421 Slice outChans(min(cStartChan(iIF),cEndChan(iIF))-1, nChan);
422 cDataSel(iIF) = Slicer(outPols, outChans);
423 }
424
425 // Get spectral data?
426 cGetSpectra = getSpectra;
427
428 // Get cross-polarization data?
429 cGetXPol = cGetXPol && getXPol;
430
431 // Get feed positions? (Not available.)
432 cGetFeedPos = False;
433
434 return maxNChan;
435}
436
437//---------------------------------------------------- PKSMS2reader::findRange
438
439// Find the range of the data in time and position.
440
441Int PKSMS2reader::findRange(
442 Int &nRow,
443 Int &nSel,
444 Vector<Double> &timeSpan,
445 Matrix<Double> &positions)
446{
447 if (!cMSopen) {
448 return 1;
449 }
450
451 nRow = cNRow;
452
453 // Find the number of rows selected.
454 nSel = 0;
455 Vector<Bool> sel(nRow);
456 for (Int irow = 0; irow < nRow; irow++) {
457 if ((sel(irow) = cBeams(cBeamNoCol(irow)) &&
458 cIFs(cDataDescIdCol(irow)))) {
459 nSel++;
460 }
461 }
462
463 // Find the time range (s).
464 timeSpan.resize(2);
465 timeSpan(0) = cTimeCol(0);
466 timeSpan(1) = cTimeCol(nRow-1);
467
468 // Retrieve positions for selected data.
469 Int isel = 0;
470 positions.resize(2,nSel);
471 for (Int irow = 0; irow < nRow; irow++) {
472 if (sel(irow)) {
473 Matrix<Double> pointingDir = cPointingCol(cFieldIdCol(irow));
474 positions.column(isel++) = pointingDir.column(0);
475 }
476 }
477
478 return 0;
479}
480
481//--------------------------------------------------------- PKSMS2reader::read
482
483// Read the next data record.
484
485Int PKSMS2reader::read(
486 Int &scanNo,
487 Int &cycleNo,
488 Double &mjd,
489 Double &interval,
490 String &fieldName,
491 String &srcName,
492 Vector<Double> &srcDir,
493 Vector<Double> &srcPM,
494 Double &srcVel,
495 String &obsMode,
496 Int &IFno,
497 Double &refFreq,
498 Double &bandwidth,
499 Double &freqInc,
500 Double &restFreq,
501 Vector<Float> &tcal,
502 String &tcalTime,
503 Float &azimuth,
504 Float &elevation,
505 Float &parAngle,
506 Float &focusAxi,
507 Float &focusTan,
508 Float &focusRot,
509 Float &temperature,
510 Float &pressure,
511 Float &humidity,
512 Float &windSpeed,
513 Float &windAz,
514 Int &refBeam,
515 Int &beamNo,
516 Vector<Double> &direction,
517 Vector<Double> &scanRate,
518 Vector<Float> &tsys,
519 Vector<Float> &sigma,
520 Vector<Float> &calFctr,
521 Matrix<Float> &baseLin,
522 Matrix<Float> &baseSub,
523 Matrix<Float> &spectra,
524 Matrix<uChar> &flagged,
525 Complex &xCalFctr,
526 Vector<Complex> &xPol)
527{
528 if (!cMSopen) {
529 return 1;
530 }
531
532 // Check for EOF.
533 if (cIdx >= cNRow) {
534 return -1;
535 }
536
537 // Find the next selected beam and IF.
538 Int ibeam;
539 Int iIF;
540 while (True) {
541 ibeam = cBeamNoCol(cIdx);
542 iIF = cDataDescIdCol(cIdx);
543 if (cBeams(ibeam) && cIFs(iIF)) {
544 break;
545 }
546
547 // Check for EOF.
548 if (++cIdx >= cNRow) {
549 return -1;
550 }
551 }
552
553 // Renumerate scan no. Here still is 1-based
554 scanNo = cScanNoCol(cIdx) - cScanNoCol(0) + 1;
555
556 if (scanNo != cScanNo) {
557 // Start of new scan.
558 cScanNo = scanNo;
559 cCycleNo = 1;
560 cTime = cTimeCol(cIdx);
561 }
562
563 Double time = cTimeCol(cIdx);
564 mjd = time/86400.0;
565 interval = cIntervalCol(cIdx);
566
567 // Reconstruct the integration cycle number; due to small latencies the
568 // integration time is usually slightly less than the time between cycles,
569 // resetting cTime will prevent the difference from accumulating.
570 cCycleNo += nint((time - cTime)/interval);
571 cycleNo = cCycleNo;
572 cTime = time;
573
574 Int fieldId = cFieldIdCol(cIdx);
575 fieldName = cFieldNameCol(fieldId);
576
577 Int srcId = cSrcIdCol(fieldId);
578 srcName = cSrcNameCol(srcId);
579 srcDir = cSrcDirCol(srcId);
580 srcPM = cSrcPMCol(srcId);
581
582 // Systemic velocity.
583 if (!cHaveSrcVel) {
584 srcVel = 0.0f;
585 } else {
586 srcVel = cSrcVelCol(srcId)(IPosition(1,0));
587 }
588
589 // Observation type.
590 Int stateId = cStateIdCol(cIdx);
591 obsMode = cObsModeCol(stateId);
592
593 IFno = iIF + 1;
594 Int nChan = abs(cEndChan(iIF) - cStartChan(iIF)) + 1;
595
596 // Minimal handling on continuum data.
597 Vector<Double> chanFreq = cChanFreqCol(iIF);
598 if (nChan == 1) {
599 cout << "The input is continuum data. "<< endl;
600 freqInc = chanFreq(0);
601 refFreq = chanFreq(0);
602 restFreq = 0.0f;
603 } else {
604 if (cStartChan(iIF) <= cEndChan(iIF)) {
605 freqInc = chanFreq(1) - chanFreq(0);
606 } else {
607 freqInc = chanFreq(0) - chanFreq(1);
608 }
609
610 refFreq = chanFreq(cRefChan(iIF)-1);
611 restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0));
612 }
613 bandwidth = abs(freqInc * nChan);
614
615 tcal.resize(cNPol(iIF));
616 tcal = 0.0f;
617 tcalTime = "";
618 azimuth = 0.0f;
619 elevation = 0.0f;
620 parAngle = 0.0f;
621 focusAxi = 0.0f;
622 focusTan = 0.0f;
623 focusRot = 0.0f;
624
625 // Find the appropriate entry in the WEATHER subtable.
626 Vector<Double> wTimes = cWeatherTimeCol.getColumn();
627 Int weatherIdx;
628 for (weatherIdx = wTimes.nelements()-1; weatherIdx >= 0; weatherIdx--) {
629 if (cWeatherTimeCol(weatherIdx) <= time) {
630 break;
631 }
632 }
633
634 if (weatherIdx < 0) {
635 // No appropriate WEATHER entry.
636 pressure = 0.0f;
637 humidity = 0.0f;
638 temperature = 0.0f;
639 } else {
640 pressure = cPressureCol(weatherIdx);
641 humidity = cHumidityCol(weatherIdx);
642 temperature = cTemperatureCol(weatherIdx);
643 }
644
645 windSpeed = 0.0f;
646 windAz = 0.0f;
647
648 refBeam = 0;
649 beamNo = ibeam + 1;
650
651 Matrix<Double> pointingDir = cPointingCol(fieldId);
652 direction = pointingDir.column(0);
653 uInt ncols = pointingDir.ncolumn();
654 if (ncols == 1) {
655 scanRate = 0.0f;
656 } else {
657 scanRate = pointingDir.column(1);
658 }
659
660 // Get Tsys assuming that entries in the SYSCAL table match the main table.
661 if (cHaveTsys) {
662 Int nTsysColRow = cTsysCol.nrow();
663 if (nTsysColRow != cNRow) {
664 cHaveTsys=0;
665 }
666 }
667 if (cHaveTsys) {
668 cTsysCol.get(cIdx, tsys, True);
669 } else {
670 Int numReceptor;
671 cNumReceptorCol.get(0, numReceptor);
672 tsys.resize(numReceptor);
673 tsys = 1.0f;
674 }
675 cSigmaCol.get(cIdx, sigma, True);
676
677 // Calibration factors (if available).
678 calFctr.resize(cNPol(iIF));
679 if (cHaveCalFctr) {
680 cCalFctrCol.get(cIdx, calFctr);
681 } else {
682 calFctr = 0.0f;
683 }
684
685 // Baseline parameters (if available).
686 if (cHaveBaseLin) {
687 baseLin.resize(2,cNPol(iIF));
688 cBaseLinCol.get(cIdx, baseLin);
689
690 baseSub.resize(9,cNPol(iIF));
691 cBaseSubCol.get(cIdx, baseSub);
692
693 } else {
694 baseLin.resize(0,0);
695 baseSub.resize(0,0);
696 }
697
698
699 // Get spectral data.
700 if (cGetSpectra) {
701 Matrix<Float> tmpData;
702 Matrix<Bool> tmpFlag;
703 cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
704 cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
705
706 // Transpose spectra.
707 Int nPol = tmpData.nrow();
708 spectra.resize(nChan, nPol);
709 flagged.resize(nChan, nPol);
710 if (cEndChan(iIF) >= cStartChan(iIF)) {
711 // Simple transposition.
712 for (Int ipol = 0; ipol < nPol; ipol++) {
713 for (Int ichan = 0; ichan < nChan; ichan++) {
714 spectra(ichan,ipol) = tmpData(ipol,ichan);
715 flagged(ichan,ipol) = tmpFlag(ipol,ichan);
716 }
717 }
718
719 } else {
720 // Transpose with inversion.
721 Int jchan = nChan - 1;
722 for (Int ipol = 0; ipol < nPol; ipol++) {
723 for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
724 spectra(ichan,ipol) = tmpData(ipol,jchan);
725 flagged(ichan,ipol) = tmpFlag(ipol,jchan);
726 }
727 }
728 }
729 }
730
731 // Get cross-polarization data.
732 if (cGetXPol) {
733 if (cHaveXCalFctr) {
734 cXCalFctrCol.get(cIdx, xCalFctr);
735 } else {
736 xCalFctr = Complex(0.0f, 0.0f);
737 }
738
739 cDataCol.get(cIdx, xPol, True);
740
741 if (cEndChan(iIF) < cStartChan(iIF)) {
742 Complex ctmp;
743 Int jchan = nChan - 1;
744 for (Int ichan = 0; ichan < nChan/2; ichan++, jchan--) {
745 ctmp = xPol(ichan);
746 xPol(ichan) = xPol(jchan);
747 xPol(jchan) = ctmp;
748 }
749 }
750 }
751
752 cIdx++;
753
754 return 0;
755}
756
757//--------------------------------------------------------- PKSMS2reader::read
758
759// Read the next data record, just the basics.
760
761Int PKSMS2reader::read(
762 Int &IFno,
763 Vector<Float> &tsys,
764 Vector<Float> &calFctr,
765 Matrix<Float> &baseLin,
766 Matrix<Float> &baseSub,
767 Matrix<Float> &spectra,
768 Matrix<uChar> &flagged)
769{
770 if (!cMSopen) {
771 return 1;
772 }
773
774 // Check for EOF.
775 if (cIdx >= cNRow) {
776 return -1;
777 }
778
779 // Find the next selected beam and IF.
780 Int ibeam;
781 Int iIF;
782 while (True) {
783 ibeam = cBeamNoCol(cIdx);
784 iIF = cDataDescIdCol(cIdx);
785 if (cBeams(ibeam) && cIFs(iIF)) {
786 break;
787 }
788
789 // Check for EOF.
790 if (++cIdx >= cNRow) {
791 return -1;
792 }
793 }
794
795 IFno = iIF + 1;
796
797 // Get Tsys assuming that entries in the SYSCAL table match the main table.
798 cTsysCol.get(cIdx, tsys, True);
799
800 // Calibration factors (if available).
801 if (cHaveCalFctr) {
802 cCalFctrCol.get(cIdx, calFctr, True);
803 } else {
804 calFctr.resize(cNPol(iIF));
805 calFctr = 0.0f;
806 }
807
808 // Baseline parameters (if available).
809 if (cHaveBaseLin) {
810 baseLin.resize(2,cNPol(iIF));
811 cBaseLinCol.get(cIdx, baseLin);
812
813 baseSub.resize(9,cNPol(iIF));
814 cBaseSubCol.get(cIdx, baseSub);
815
816 } else {
817 baseLin.resize(0,0);
818 baseSub.resize(0,0);
819 }
820
821 if (cGetSpectra) {
822 // Get spectral data.
823 Matrix<Float> tmpData;
824 Matrix<Bool> tmpFlag;
825 cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
826 cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
827
828 // Transpose spectra.
829 Int nChan = tmpData.ncolumn();
830 Int nPol = tmpData.nrow();
831 spectra.resize(nChan, nPol);
832 flagged.resize(nChan, nPol);
833 if (cEndChan(iIF) >= cStartChan(iIF)) {
834 // Simple transposition.
835 for (Int ipol = 0; ipol < nPol; ipol++) {
836 for (Int ichan = 0; ichan < nChan; ichan++) {
837 spectra(ichan,ipol) = tmpData(ipol,ichan);
838 flagged(ichan,ipol) = tmpFlag(ipol,ichan);
839 }
840 }
841
842 } else {
843 // Transpose with inversion.
844 Int jchan = nChan - 1;
845 for (Int ipol = 0; ipol < nPol; ipol++) {
846 for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
847 spectra(ichan,ipol) = tmpData(ipol,jchan);
848 flagged(ichan,ipol) = tmpFlag(ipol,jchan);
849 }
850 }
851 }
852 }
853
854 cIdx++;
855
856 return 0;
857}
858
859//-------------------------------------------------------- PKSMS2reader::close
860
861// Close the MS.
862
863void PKSMS2reader::close()
864{
865 cPKSMS = MeasurementSet();
866 cMSopen = False;
867}
Note: See TracBrowser for help on using the repository browser.