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

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

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