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

Last change on this file since 1427 was 1427, checked in by Malte Marquarding, 16 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.