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

Last change on this file since 1325 was 1325, checked in by mar637, 17 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
Line 
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.