source: trunk/external-alma/atnf/PKSIO/PKSMS2reader.cc @ 3067

Last change on this file since 3067 was 2343, checked in by Takeshi Nakazato, 13 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

Patch from Wes.
This is changes intend to compile with llvm c++ and Intel
and those might not affect actual behavior.
File size: 42.5 KB
Line 
1//#---------------------------------------------------------------------------
2//# PKSMS2reader.cc: Class to read Parkes Multibeam data from a v2 MS.
3//#---------------------------------------------------------------------------
4//# livedata - processing pipeline for single-dish, multibeam spectral data.
5//# Copyright (C) 2000-2009, Australia Telescope National Facility, CSIRO
6//#
7//# This file is part of livedata.
8//#
9//# livedata is free software: you can redistribute it and/or modify it under
10//# the terms of the GNU General Public License as published by the Free
11//# Software Foundation, either version 3 of the License, or (at your option)
12//# any later version.
13//#
14//# livedata is distributed in the hope that it will be useful, but WITHOUT
15//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16//# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
17//# more details.
18//#
19//# You should have received a copy of the GNU General Public License along
20//# with livedata.  If not, see <http://www.gnu.org/licenses/>.
21//#
22//# Correspondence concerning livedata may be directed to:
23//#        Internet email: mcalabre@atnf.csiro.au
24//#        Postal address: Dr. Mark Calabretta
25//#                        Australia Telescope National Facility, CSIRO
26//#                        PO Box 76
27//#                        Epping NSW 1710
28//#                        AUSTRALIA
29//#
30//# http://www.atnf.csiro.au/computing/software/livedata.html
31//# $Id: PKSMS2reader.cc,v 19.23 2009-09-29 07:33:38 cal103 Exp $
32//#---------------------------------------------------------------------------
33//# Original: 2000/08/03, Mark Calabretta, ATNF
34//#---------------------------------------------------------------------------
35
36// AIPS++ includes.
37#include <casa/stdio.h>
38#include <casa/Arrays/ArrayMath.h>
39#include <casa/Arrays/Slice.h>
40#include <ms/MeasurementSets/MSColumns.h>
41#include <tables/Tables.h>
42#include <casa/Quanta/MVTime.h>
43#include <casa/Quanta/MVAngle.h>
44#include <casa/BasicMath/Math.h>
45#include <casa/Logging/LogIO.h>
46#include <casa/Utilities/Sort.h>
47#include <measures/Measures/MeasConvert.h>
48#include <measures/Measures/MEpoch.h>
49#include <measures/Measures/MeasRef.h>
50
51
52// Parkes includes.
53#include <atnf/pks/pks_maths.h>
54#include <atnf/PKSIO/PKSrecord.h>
55#include <atnf/PKSIO/PKSMS2reader.h>
56
57
58//------------------------------------------------- PKSMS2reader::PKSMS2reader
59
60// Default constructor.
61
62PKSMS2reader::PKSMS2reader()
63{
64  cMSopen = False;
65}
66
67//------------------------------------------------ PKSMS2reader::~PKSMS2reader
68
69PKSMS2reader::~PKSMS2reader()
70{
71  close();
72}
73
74//--------------------------------------------------------- PKSMS2reader::open
75
76// Open the MS for reading.
77
78Int PKSMS2reader::open(
79        const String msName,
80        const String antenna,
81        Vector<Bool> &beams,
82        Vector<Bool> &IFs,
83        Vector<uInt> &nChan,
84        Vector<uInt> &nPol,
85        Vector<Bool> &haveXPol,
86        Bool   &haveBase,
87        Bool   &haveSpectra)
88{
89  // Check that MS is readable.
90  if (!MS::isReadable(msName)) {
91    return 1;
92  }
93
94  if (cMSopen) {
95    close();
96  }
97
98  cPKSMS  = MeasurementSet(msName);
99
100  // data selection by antenna
101  if ( antenna.length() == 0 ) {
102    cAntId.resize( 1 ) ;
103    //cAntId[0] = 0 ;
104    ROScalarColumn<Int> ant1Col( cPKSMS, "ANTENNA1" ) ;
105    cAntId[0] = ant1Col(0) ;
106  }
107  else {
108    setupAntennaList( antenna ) ;
109    if ( cAntId.size() > 1 ) {
110      LogIO os( LogOrigin( "PKSMS2reader", "open()", WHERE ) ) ;
111      os << LogIO::WARN << "PKSMS2reader is not ready for multiple antenna selection. Use first antenna id " << cAntId[0] << "."<< LogIO::POST ;
112      Int tmp = cAntId[0] ;
113      cAntId.resize( 1 ) ;
114      cAntId[0] = tmp ;
115    }
116    stringstream ss ;
117    ss << "SELECT FROM $1 WHERE ANTENNA1 == ANTENNA2 && ANTENNA1 IN [" ;
118    for ( uInt i = 0 ; i < cAntId.size() ; i++ ) {
119      ss << cAntId[i] ;
120      if ( i == cAntId.size()-1 ) {
121        ss << "]" ;
122      }
123      else {
124        ss << "," ;
125      }
126    }
127    string taql = ss.str() ;
128    //cerr << "taql = " << taql << endl ;
129    cPKSMS = MeasurementSet( tableCommand( taql, cPKSMS ) ) ;
130  }
131
132  // taql access to the syscal table
133  cHaveSysCal = False;
134  if (cHaveSysCal=Table::isReadable(cPKSMS.sysCalTableName())) {
135    cSysCalTab = Table(cPKSMS.sysCalTableName());
136  }
137
138  // Lock the table for read access.
139  cPKSMS.lock(False);
140
141  cIdx    = 0;
142  lastmjd = 0.0;
143  cNRow   = cPKSMS.nrow();
144  cMSopen = True;
145
146  // Main MS table and subtable column access.
147  ROMSMainColumns         msCols(cPKSMS);
148  ROMSDataDescColumns     dataDescCols(cPKSMS.dataDescription());
149  ROMSFeedColumns         feedCols(cPKSMS.feed());
150  ROMSFieldColumns        fieldCols(cPKSMS.field());
151  ROMSPointingColumns     pointingCols(cPKSMS.pointing());
152  ROMSPolarizationColumns polarizationCols(cPKSMS.polarization());
153  ROMSSourceColumns       sourceCols(cPKSMS.source());
154  ROMSSpWindowColumns     spWinCols(cPKSMS.spectralWindow());
155  ROMSStateColumns        stateCols(cPKSMS.state());
156  ROMSSysCalColumns       sysCalCols(cPKSMS.sysCal());
157  ROMSWeatherColumns      weatherCols(cPKSMS.weather());
158  ROMSAntennaColumns      antennaCols(cPKSMS.antenna());
159
160  // Column accessors for required columns.
161  cScanNoCol.reference(msCols.scanNumber());
162  cTimeCol.reference(msCols.time());
163  cIntervalCol.reference(msCols.interval());
164
165  cFieldIdCol.reference(msCols.fieldId());
166  cFieldNameCol.reference(fieldCols.name());
167  cFieldDelayDirCol.reference(fieldCols.delayDir());
168
169  cSrcIdCol.reference(fieldCols.sourceId());
170  cSrcId2Col.reference(sourceCols.sourceId());
171  cSrcNameCol.reference(sourceCols.name());
172  cSrcDirCol.reference(sourceCols.direction());
173  cSrcPMCol.reference(sourceCols.properMotion());
174  cSrcRestFrqCol.reference(sourceCols.restFrequency());
175
176  cStateIdCol.reference(msCols.stateId());
177  cObsModeCol.reference(stateCols.obsMode());
178  cCalCol.reference(stateCols.cal());
179  cSigStateCol.reference(stateCols.sig());
180  cRefStateCol.reference(stateCols.ref());
181
182  cDataDescIdCol.reference(msCols.dataDescId());
183  cSpWinIdCol.reference(dataDescCols.spectralWindowId());
184  cChanFreqCol.reference(spWinCols.chanFreq());
185  cTotBWCol.reference(spWinCols.totalBandwidth());
186
187  cWeatherTimeCol.reference(weatherCols.time());
188  cTemperatureCol.reference(weatherCols.temperature());
189  cPressureCol.reference(weatherCols.pressure());
190  cHumidityCol.reference(weatherCols.relHumidity());
191
192  cBeamNoCol.reference(msCols.feed1());
193  cPointingCol.reference(pointingCols.direction());
194  cPointingTimeCol.reference(pointingCols.time());
195  cSigmaCol.reference(msCols.sigma());
196  cNumReceptorCol.reference(feedCols.numReceptors());
197
198  // Optional columns.
199  cHaveTsys = False;
200  cHaveTcal = False;
201  if ((cHaveSrcVel = cPKSMS.source().tableDesc().isColumn("SYSVEL"))) {
202    cSrcVelCol.attach(cPKSMS.source(), "SYSVEL");
203  }
204
205  if (cHaveSysCal && (cHaveTsys = cPKSMS.sysCal().tableDesc().isColumn("TSYS"))) {
206    cTsysCol.attach(cPKSMS.sysCal(), "TSYS");
207  }
208 
209  if (cHaveSysCal && (cHaveTcal = cPKSMS.sysCal().tableDesc().isColumn("TCAL"))) {
210    cTcalCol.attach(cPKSMS.sysCal(), "TCAL");
211  }
212
213  if ((cHaveCalFctr = cPKSMS.tableDesc().isColumn("CALFCTR"))) {
214    cCalFctrCol.attach(cPKSMS, "CALFCTR");
215  }
216
217  if ((cHaveBaseLin = cPKSMS.tableDesc().isColumn("BASELIN"))) {
218    cBaseLinCol.attach(cPKSMS, "BASELIN");
219    cBaseSubCol.attach(cPKSMS, "BASESUB");
220  }
221
222  // Spectral data should always be present.
223  haveSpectra = True;
224  cHaveDataCol = False;
225  cHaveCorrectedDataCol = False;
226  ROMSObservationColumns observationCols(cPKSMS.observation());
227  //String telName = observationCols.telescopeName()(0);
228  cTelName = observationCols.telescopeName()(0);
229  //cATF = cTelName.contains("ATF");
230  //cOSF = cTelName.contains("OSF");
231  //cALMA = cTelName.contains("ALMA");
232  cALMA = cTelName.contains("ATF")||cTelName.contains("OSF")||
233           cTelName.contains("ALMA");
234
235  if (cHaveDataCol = cPKSMS.isColumn(MSMainEnums::DATA)) {
236    if (cALMA) {
237      //try to read a single baseline interferometeric data
238      //and treat it as single dish data
239      //maybe extended for ALMA commissioning later
240      cDataCol.reference(msCols.data());
241      if (cHaveCorrectedDataCol = cPKSMS.isColumn(MSMainEnums::CORRECTED_DATA)) {
242        //cerr<<"Do have CORRECTED_DATA column"<<endl;
243        cCorrectedDataCol.reference(msCols.correctedData());
244      }
245    }
246  }
247  else {
248    cFloatDataCol.reference(msCols.floatData());
249  }
250  cFlagCol.reference(msCols.flag());
251  cFlagRowCol.reference(msCols.flagRow());
252
253  if (cGetXPol = (cPKSMS.isColumn(MSMainEnums::DATA) && (!cALMA))) {
254    if ((cHaveXCalFctr = cPKSMS.tableDesc().isColumn("XCALFCTR"))) {
255      cXCalFctrCol.attach(cPKSMS, "XCALFCTR");
256    }
257
258    cDataCol.reference(msCols.data());
259  }
260
261  // Find which beams are present in the data.
262  Vector<Int> beamNos = cBeamNoCol.getColumn();
263  Int maxBeamNo = max(beamNos) + 1;
264  beams.resize(maxBeamNo);
265
266  beams = False;
267  for (uInt irow = 0; irow < beamNos.nelements(); irow++) {
268    beams(beamNos(irow)) = True;
269  }
270
271
272  // Number of IFs.
273  //uInt nIF = dataDescCols.nrow();
274  uInt nIF =spWinCols.nrow();
275  Vector<Int> spWinIds = cSpWinIdCol.getColumn() ;
276  IFs.resize(nIF);
277  IFs = True;
278  for ( Int ispw = 0 ; ispw < nIF ; ispw++ ) {
279    if ( allNE( ispw, spWinIds ) ) {
280      IFs(ispw) = False ;
281    }
282  }
283
284  // Number of polarizations and channels in each IF.
285  ROScalarColumn<Int> numChanCol(spWinCols.numChan());
286
287  ROScalarColumn<Int> polIdCol(dataDescCols.polarizationId());
288  ROScalarColumn<Int> numPolCol(polarizationCols.numCorr());
289
290  nChan.resize(nIF);
291  nPol.resize(nIF);
292  for (uInt iIF = 0; iIF < nIF; iIF++) {
293    if ( IFs(iIF) ) {
294      nChan(iIF) = numChanCol(cSpWinIdCol(iIF)) ;
295      nPol(iIF) = numPolCol(polIdCol(iIF)) ;
296    }
297    else {
298      nChan(iIF) = 0 ;
299      nPol(iIF) = 0 ;
300    }
301  }
302
303  // Cross-polarization data present?
304  haveXPol.resize(nIF);
305  haveXPol = False;
306
307  if (cGetXPol && !(cALMA)) {
308    for (Int irow = 0; irow < cNRow; irow++) {
309      if (cDataCol.isDefined(irow)) {
310        Int iIF = cDataDescIdCol(irow);
311        haveXPol(iIF) = True;
312      }
313    }
314  }
315
316
317  // Initialize member data.
318  cBeams.assign(beams);
319  cIFs.assign(IFs);
320  cNChan.assign(nChan);
321  cNPol.assign(nPol);
322  cHaveXPol.assign(haveXPol);
323
324
325  // Default channel range selection.
326  cStartChan.resize(nIF);
327  cEndChan.resize(nIF);
328  cRefChan.resize(nIF);
329
330  for (uInt iIF = 0; iIF < nIF; iIF++) {
331    cStartChan(iIF) = 1;
332    cEndChan(iIF)   = cNChan(iIF);
333    cRefChan(iIF)   = cNChan(iIF)/2 + 1;
334  }
335
336  Slice all;
337  cDataSel.resize(nIF);
338  cDataSel = Slicer(all, all);
339
340  cScanNo  = 0;
341  cCycleNo = 1;
342  cTime    = cTimeCol(0);
343
344  return 0;
345}
346
347//---------------------------------------------------- PKSMS2reader::getHeader
348
349// Get parameters describing the data.
350
351Int PKSMS2reader::getHeader(
352        String &observer,
353        String &project,
354        String &antName,
355        Vector<Double> &antPosition,
356        // before merge...
357        //String &obsMode,
358        String &obsType,
359        String &bunit,
360        Float  &equinox,
361        String &dopplerFrame,
362        Double &mjd,
363        Double &refFreq,
364        Double &bandwidth)
365{
366  if (!cMSopen) {
367    return 1;
368  }
369
370  // Observer and project.
371  ROMSObservationColumns observationCols(cPKSMS.observation());
372  observer = observationCols.observer()(0);
373  project  = observationCols.project()(0);
374
375  // Antenna name and ITRF coordinates.
376  ROMSAntennaColumns antennaCols(cPKSMS.antenna());
377  //antName = antennaCols.name()(0);
378  antName = antennaCols.name()(cAntId[0]);
379  if (cALMA) {
380     antName = cTelName + "-" + antName;
381  }
382  //antPosition = antennaCols.position()(0);
383  antPosition = antennaCols.position()(cAntId[0]);
384
385  // Observation type.
386  if (cObsModeCol.nrow()) {
387    obsType = cObsModeCol(0);
388    if (obsType == "\0") obsType = "RF";
389  } else {
390    obsType = "RF";
391  }
392
393  bunit = "";
394  if (cHaveDataCol) {
395    const TableRecord& keywordSet2
396       = cDataCol.columnDesc().keywordSet();
397    if(keywordSet2.isDefined("UNIT")) {
398      bunit = keywordSet2.asString("UNIT");
399    }
400  } else {
401    const TableRecord& keywordSet
402       = cFloatDataCol.columnDesc().keywordSet();
403    if(keywordSet.isDefined("UNIT")) {
404      bunit = keywordSet.asString("UNIT");
405    }
406  }
407
408/***
409  const TableRecord& keywordSet
410       = cFloatDataCol.columnDesc().keywordSet();
411  if(keywordSet.isDefined("UNIT")) {
412    fluxunit = keywordSet.asString("UNIT");
413  }
414***/
415  // Coordinate equinox.
416  ROMSPointingColumns pointingCols(cPKSMS.pointing());
417  String dirref = pointingCols.direction().keywordSet().asRecord("MEASINFO").
418                    asString("Ref");
419  cDirRef = dirref;
420  if (dirref =="AZELGEO" || dirref == "AZEL") {
421     dirref = "J2000";
422  }
423  sscanf(dirref.chars()+1, "%f", &equinox);
424
425  // Frequency/velocity reference frame.
426  ROMSSpWindowColumns spWinCols(cPKSMS.spectralWindow());
427  dopplerFrame = MFrequency::showType(spWinCols.measFreqRef()(0));
428
429  // Translate to FITS standard names.
430  if (dopplerFrame == "TOPO") {
431    dopplerFrame = "TOPOCENT";
432  } else if (dopplerFrame == "GEO") {
433    dopplerFrame = "GEOCENTR";
434  } else if (dopplerFrame == "BARY") {
435    dopplerFrame = "BARYCENT";
436  } else if (dopplerFrame == "GALACTO") {
437    dopplerFrame = "GALACTOC";
438  } else if (dopplerFrame == "LGROUP") {
439    dopplerFrame = "LOCALGRP";
440  } else if (dopplerFrame == "CMB") {
441    dopplerFrame = "CMBDIPOL";
442  } else if (dopplerFrame == "REST") {
443    dopplerFrame = "SOURCE";
444  }
445
446  // MJD at start of observation.
447  mjd = cTimeCol(0)/86400.0;
448
449  // Reference frequency and bandwidth.
450  refFreq   = spWinCols.refFrequency()(0);
451  bandwidth = spWinCols.totalBandwidth()(0);
452
453  return 0;
454}
455
456//-------------------------------------------------- PKSMS2reader::getFreqInfo
457
458// Get frequency parameters for each IF.
459
460Int PKSMS2reader::getFreqInfo(
461        Vector<Double> &startFreq,
462        Vector<Double> &endFreq)
463{
464  uInt nIF = cIFs.nelements();
465  startFreq.resize(nIF);
466  endFreq.resize(nIF);
467
468  for (uInt iIF = 0; iIF < nIF; iIF++) {
469    Vector<Double> chanFreq = cChanFreqCol(iIF);
470
471    Int nChan = chanFreq.nelements();
472    startFreq(iIF) = chanFreq(0);
473    endFreq(iIF)   = chanFreq(nChan-1);
474  }
475
476  return 0;
477}
478
479//------------------------------------------------------- PKSMS2reader::select
480
481// Set data selection by beam number and channel.
482
483uInt PKSMS2reader::select(
484        const Vector<Bool> beamSel,
485        const Vector<Bool> IFsel,
486        const Vector<Int>  startChan,
487        const Vector<Int>  endChan,
488        const Vector<Int>  refChan,
489        const Bool getSpectra,
490        const Bool getXPol,
491        const Bool getFeedPos,
492        const Bool getPointing,
493        const Int  coordSys)
494{
495  if (!cMSopen) {
496    return 1;
497  }
498
499  // Beam selection.
500  uInt nBeam = cBeams.nelements();
501  uInt nBeamSel = beamSel.nelements();
502  for (uInt ibeam = 0; ibeam < nBeam; ibeam++) {
503    if (ibeam < nBeamSel) {
504      cBeams(ibeam) = beamSel(ibeam);
505    } else {
506      cBeams(ibeam) = False;
507    }
508  }
509
510  uInt nIF = cIFs.nelements();
511  uInt maxNChan = 0;
512  for (uInt iIF = 0; iIF < nIF; iIF++) {
513    // IF selection.
514    if (iIF < IFsel.nelements()) {
515      cIFs(iIF) = IFsel(iIF);
516    } else {
517      cIFs(iIF) = False;
518    }
519
520    if (!cIFs(iIF)) continue;
521
522
523    // Channel selection.
524    if (iIF < startChan.nelements()) {
525      cStartChan(iIF) = startChan(iIF);
526
527      if (cStartChan(iIF) <= 0) {
528        cStartChan(iIF) += cNChan(iIF);
529      } else if (cStartChan(iIF) > Int(cNChan(iIF))) {
530        cStartChan(iIF)  = cNChan(iIF);
531      }
532    }
533
534    if (iIF < endChan.nelements()) {
535      cEndChan(iIF) = endChan(iIF);
536
537      if (cEndChan(iIF) <= 0) {
538        cEndChan(iIF) += cNChan(iIF);
539      } else if (cEndChan(iIF) > Int(cNChan(iIF))) {
540        cEndChan(iIF)  = cNChan(iIF);
541      }
542    }
543
544    if (iIF < refChan.nelements()) {
545      cRefChan(iIF) = refChan(iIF);
546    } else {
547      cRefChan(iIF) = cStartChan(iIF);
548      if (cStartChan(iIF) <= cEndChan(iIF)) {
549        cRefChan(iIF) += (cEndChan(iIF) - cStartChan(iIF) + 1)/2;
550      } else {
551        cRefChan(iIF) -= (cStartChan(iIF) - cEndChan(iIF) + 1)/2;
552      }
553    }
554
555    uInt nChan = abs(cEndChan(iIF) - cStartChan(iIF)) + 1;
556    if (maxNChan < nChan) {
557      maxNChan = nChan;
558    }
559
560    // Inverted Slices are not allowed.
561    Slice outPols;
562    Slice outChans(min(cStartChan(iIF),cEndChan(iIF))-1, nChan);
563    cDataSel(iIF) = Slicer(outPols, outChans);
564  }
565
566  // Get spectral data?
567  cGetSpectra = getSpectra;
568
569  // Get cross-polarization data?
570  cGetXPol = cGetXPol && getXPol;
571
572  // Get feed positions?  (Not available.)
573  cGetFeedPos = False;
574
575  // Get Pointing data (for MS)
576  cGetPointing = getPointing;
577
578  // Coordinate system?  (Only equatorial available.)
579  cCoordSys = 0;
580
581  return maxNChan;
582}
583
584//---------------------------------------------------- PKSMS2reader::findRange
585
586// Find the range of the data in time and position.
587
588Int PKSMS2reader::findRange(
589        Int    &nRow,
590        Int    &nSel,
591        Vector<Double> &timeSpan,
592        Matrix<Double> &positions)
593{
594  if (!cMSopen) {
595    return 1;
596  }
597
598  nRow = cNRow;
599
600  // Find the number of rows selected.
601  nSel = 0;
602  Vector<Bool> sel(nRow);
603  for (Int irow = 0; irow < nRow; irow++) {
604    if ((sel(irow) = cBeams(cBeamNoCol(irow)) &&
605                     cIFs(cDataDescIdCol(irow)))) {
606      nSel++;
607    }
608  }
609
610  // Find the time range (s).
611  timeSpan.resize(2);
612  timeSpan(0) = cTimeCol(0);
613  timeSpan(1) = cTimeCol(nRow-1);
614
615  // Retrieve positions for selected data.
616  Int isel = 0;
617  positions.resize(2,nSel);
618  for (Int irow = 0; irow < nRow; irow++) {
619    if (sel(irow)) {
620      Matrix<Double> pointingDir = cPointingCol(cFieldIdCol(irow));
621      positions.column(isel++) = pointingDir.column(0);
622    }
623  }
624
625  return 0;
626}
627
628//--------------------------------------------------------- PKSMS2reader::read
629
630// Read the next data record.
631
632/**
633Int PKSMS2reader::read(
634        Int             &scanNo,
635        Int             &cycleNo,
636        Double          &mjd,
637        Double          &interval,
638        String          &fieldName,
639        String          &srcName,
640        Vector<Double>  &srcDir,
641        Vector<Double>  &srcPM,
642        Double          &srcVel,
643        String          &obsMode,
644        Int             &IFno,
645        Double          &refFreq,
646        Double          &bandwidth,
647        Double          &freqInc,
648        Vector<Double>  &restFreq,
649        Vector<Float>   &tcal,
650        String          &tcalTime,
651        Float           &azimuth,
652        Float           &elevation,
653        Float           &parAngle,
654        Float           &focusAxi,
655        Float           &focusTan,
656        Float           &focusRot,
657        Float           &temperature,
658        Float           &pressure,
659        Float           &humidity,
660        Float           &windSpeed,
661        Float           &windAz,
662        Int             &refBeam,
663        Int             &beamNo,
664        Vector<Double>  &direction,
665        Vector<Double>  &scanRate,
666        Vector<Float>   &tsys,
667        Vector<Float>   &sigma,
668        Vector<Float>   &calFctr,
669        Matrix<Float>   &baseLin,
670        Matrix<Float>   &baseSub,
671        Matrix<Float>   &spectra,
672        Matrix<uChar>   &flagged,
673        uInt            &flagrow,
674        Complex         &xCalFctr,
675        Vector<Complex> &xPol)
676**/
677Int PKSMS2reader::read(PKSrecord &pksrec)
678{
679  LogIO os( LogOrigin( "PKSMS2reader", "read()", WHERE ) ) ;
680
681  if (!cMSopen) {
682    return 1;
683  }
684
685  // Check for EOF.
686  if (cIdx >= cNRow) {
687    return -1;
688  }
689
690  // Find the next selected beam and IF.
691  Int ibeam;
692  Int iIF;
693  Int iDataDesc;
694
695  while (True) {
696    ibeam = cBeamNoCol(cIdx);
697    iDataDesc   = cDataDescIdCol(cIdx);
698    iIF   =cSpWinIdCol(iDataDesc);
699    if (cBeams(ibeam) && cIFs(iIF)) {
700      break;
701    }
702
703    // Check for EOF.
704    if (++cIdx >= cNRow) {
705      return -1;
706    }
707  }
708  // Renumerate scan no. Here still is 1-based
709  //scanNo = cScanNoCol(cIdx) - cScanNoCol(0) + 1;
710  //scanNo = cScanNoCol(cIdx);
711  pksrec.scanNo = cScanNoCol(cIdx);
712
713  if (pksrec.scanNo != cScanNo) {
714    // Start of new scan.
715    cScanNo  = pksrec.scanNo;
716    cCycleNo = 1;
717    cTime    = cTimeCol(cIdx);
718  }
719
720  Double time = cTimeCol(cIdx);
721  pksrec.mjd      = time/86400.0;
722  pksrec.interval = cIntervalCol(cIdx);
723
724  // Reconstruct the integration cycle number; due to small latencies the
725  // integration time is usually slightly less than the time between cycles,
726  // resetting cTime will prevent the difference from accumulating.
727  cCycleNo += nint((time - cTime)/pksrec.interval);
728  pksrec.cycleNo = cCycleNo;
729  cTime = time;
730
731  Int fieldId = cFieldIdCol(cIdx);
732  pksrec.fieldName = cFieldNameCol(fieldId);
733
734  Int srcId = cSrcIdCol(fieldId);
735  //For source with multiple spectral window setting, this is not
736  // correct. Source name of srcId may not be at 'srcId'th row of SrcNameCol
737  //srcName = cSrcNameCol(srcId);
738  for (uInt irow = 0; irow < cSrcId2Col.nrow(); irow++) {
739    if (cSrcId2Col(irow) == srcId) {
740      //srcName = cSrcNameCol(irow);
741      pksrec.srcName = cSrcNameCol(irow);
742    }
743  }
744
745  pksrec.srcDir  = cSrcDirCol(srcId);
746  pksrec.srcPM   = cSrcPMCol(srcId);
747
748  // Systemic velocity.
749  if (!cHaveSrcVel || cALMA) {
750    pksrec.srcVel = 0.0f;
751  } else {
752    pksrec.srcVel = cSrcVelCol(srcId)(IPosition(1,0));
753  }
754
755  ROMSAntennaColumns antennaCols(cPKSMS.antenna());
756  //String telescope = antennaCols.name()(0);
757  String telescope = antennaCols.name()(cAntId[0]);
758  Bool cGBT = telescope.contains("GBT");
759  //Bool cPM = telescope.contains("PM"); // ACA TP antenna
760  //Bool cDV = telescope.contains("DV"); // VERTEX
761  //Bool cCM = telescope.contains("CM"); // ACA 7m antenna
762  //Bool cALMA = cPM || cDV || cCM ;
763  // Observation type.
764  // check if State Table exist
765  //Bool cHaveStateTab=Table::isReadable(cPKSMS.stateTableName());
766  Int stateId = 0;
767  Int StateNRow = 0;
768  StateNRow=cObsModeCol.nrow();
769  if (Table::isReadable(cPKSMS.stateTableName())) {
770        pksrec.obsType = " ";
771    if (StateNRow > 0) {
772      stateId = cStateIdCol(cIdx);
773      if (stateId == -1) {
774        //pksrec.obsType = " ";
775      } else {
776        pksrec.obsType = cObsModeCol(stateId);
777        Bool sigState =cSigStateCol(stateId);
778        Bool refState =cRefStateCol(stateId);
779        //DEBUG
780        //cerr <<"stateid="<<stateId<<" obsmode="<<pksrec.obsType<<endl;
781        if (cGBT) {
782          // split the obsType string and append a proper label
783          // (these are GBT specific)
784          int epos = pksrec.obsType.find_first_of(':');
785          int nextpos = pksrec.obsType.find_first_of(':',epos+1);
786          string obsMode1 = pksrec.obsType.substr(0,epos);
787          string obsMode2 = pksrec.obsType.substr(epos+1,nextpos-epos-1);
788     
789          //cerr <<"obsMode2= "<<obsMode2<<endl;
790          if (!pksrec.srcName.contains("_ps")
791              &&!pksrec.srcName.contains("_psr")
792              &&!pksrec.srcName.contains("_nod")
793              &&!pksrec.srcName.contains("_fs")
794              &&!pksrec.srcName.contains("_fsr")) {
795            // if Nod mode observation , append '_nod'
796            if (obsMode1 == "Nod") {
797              //pksrec.srcName.append("_nod");
798              pksrec.srcType = SrcType::NOD ;
799            } else if (obsMode1 == "OffOn") {
800            // for GBT position switch observations (OffOn or OnOff)
801              //if (obsMode2 == "PSWITCHON") pksrec.srcName.append("_ps");
802              //if (obsMode2 == "PSWITCHOFF") pksrec.srcName.append("_psr");
803              if (obsMode2 == "PSWITCHON") pksrec.srcType = SrcType::PSON ;
804              if (obsMode2 == "PSWITCHOFF") pksrec.srcType = SrcType::PSOFF ;
805            } else {
806              if (obsMode2 == "FSWITCH") {
807              // for GBT frequency switch mode
808                //if (sigState) pksrec.srcName.append("_fs");
809                //if (refState) pksrec.srcName.append("_fsr");
810                if (sigState) pksrec.srcType = SrcType::FSON ;
811                if (refState) pksrec.srcType = SrcType::FSOFF ;
812              }
813            }
814          }
815        }
816        else if (cALMA) {
817          // ALMA tag
818          // split the obsType string and append a proper label
819          string substr[1] ;
820          int numSubstr = split( pksrec.obsType, substr, 1, "," );
821          String obsType = String( substr[0] );
822          int epos = obsType.find_first_of('.');
823          int nextpos = obsType.find_first_of('.',epos+1);
824          string obsMode1 = obsType.substr(0,epos);
825          string obsMode2 = obsType.substr(epos+1,nextpos-epos-1);
826     
827          //cerr <<"obsMode2= "<<obsMode2<<endl;
828          // Current OBS_MODE format:
829          //
830          //     ON: OBSERVE_TARGET.ON_SOURCE
831          //    OFF: OBSERVE_TARGET.OFF_SOURCE
832          //
833          if (obsMode1 == "OBSERVE_TARGET") {
834            //if (obsMode2 == "ON_SOURCE") pksrec.srcName.append("_pson");
835            //if (obsMode2 == "OFF_SOURCE") pksrec.srcName.append("_psoff");
836            if (obsMode2 == "ON_SOURCE") pksrec.srcType = SrcType::PSON ;
837            if (obsMode2 == "OFF_SOURCE") pksrec.srcType = SrcType::PSOFF ;
838          }
839        }
840      }
841    }
842  }
843  // CAL state
844  // this should be apply just for GBT data?
845  Double Cal;
846  if (stateId==-1 || StateNRow==0) {
847    Cal = 0;
848  } else {
849    Cal = cCalCol(stateId);
850  }
851  if (cGBT) {
852    if (Cal > 0 && !pksrec.srcName.contains("_calon")) {
853      //pksrec.srcName.append("_calon");
854      if ( pksrec.srcType == SrcType::NOD )
855        pksrec.srcType = SrcType::NODCAL ;
856      else if ( pksrec.srcType == SrcType::PSON )
857        pksrec.srcType = SrcType::PONCAL ;
858      else if ( pksrec.srcType == SrcType::PSOFF )
859        pksrec.srcType = SrcType::POFFCAL ;
860      else if ( pksrec.srcType == SrcType::FSON )
861        pksrec.srcType = SrcType::FONCAL ;
862      else if ( pksrec.srcType == SrcType::FSOFF )
863        pksrec.srcType = SrcType::FOFFCAL ;
864      else
865        pksrec.srcName.append("_calon");
866    }
867  }
868
869  pksrec.IFno = iIF + 1;
870  Int nChan = abs(cEndChan(iIF) - cStartChan(iIF)) + 1;
871 
872  // Minimal handling on continuum data.
873  Vector<Double> chanFreq = cChanFreqCol(iIF);
874  if (nChan == 1) {
875    //pksrec.freqInc  = chanFreq(0);
876    pksrec.freqInc  = cTotBWCol(iIF);
877    pksrec.refFreq  = chanFreq(0);
878    pksrec.restFreq.resize(1);
879    pksrec.restFreq[0] = 0.0f;
880  } else {
881 
882    if (cStartChan(iIF) <= cEndChan(iIF)) {
883      pksrec.freqInc = chanFreq(1) - chanFreq(0);
884    } else {
885      pksrec.freqInc = chanFreq(0) - chanFreq(1);
886    }
887
888    pksrec.refFreq  = chanFreq(cRefChan(iIF)-1);
889
890    Bool HaveSrcRestFreq= cSrcRestFrqCol.isDefined(srcId);
891    if (HaveSrcRestFreq) {
892      //restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0));
893      //restFreq = cSrcRestFrqCol(srcId);
894      pksrec.restFreq = cSrcRestFrqCol(srcId);
895    } else {
896      pksrec.restFreq.resize(1);
897      pksrec.restFreq[0] = 0.0f;
898    }
899  }
900  //pksrec.bandwidth = abs(pksrec.freqInc * nChan);
901  pksrec.bandwidth = abs(cTotBWCol(0));
902
903  pksrec.tcal.resize(cNPol(iIF));
904  pksrec.tcal      = 0.0f;
905  pksrec.tcalTime  = "";
906//  pksrec.azimuth   = 0.0f;
907//  pksrec.elevation = 0.0f;
908  pksrec.parAngle  = 0.0f;
909
910  pksrec.focusAxi  = 0.0f;
911  pksrec.focusTan  = 0.0f;
912  pksrec.focusRot  = 0.0f;
913
914  // Find the appropriate entry in the WEATHER subtable.
915  //Bool cHaveStateTab=Table::isReadable(cPKSMS.stateTableName());
916  Bool cHaveWeatherTab = Table::isReadable(cPKSMS.weatherTableName());
917  Int weatherIdx=-1;
918  if (cHaveWeatherTab) {
919    Vector<Double> wTimes = cWeatherTimeCol.getColumn();
920    for (weatherIdx = wTimes.nelements()-1; weatherIdx >= 0; weatherIdx--) {
921      if (cWeatherTimeCol(weatherIdx) <= time) {
922        break;
923      }
924    }
925  }
926
927  if (weatherIdx < 0 || !cHaveWeatherTab) {
928    // No appropriate WEATHER entry.
929    pksrec.temperature = 0.0f;
930    pksrec.pressure    = 0.0f;
931    pksrec.humidity    = 0.0f;
932  } else {
933    pksrec.temperature = cTemperatureCol(weatherIdx);
934    pksrec.pressure    = cPressureCol(weatherIdx);
935    pksrec.humidity    = cHumidityCol(weatherIdx);
936  }
937
938  pksrec.windSpeed = 0.0f;
939  pksrec.windAz    = 0.0f;
940
941  pksrec.refBeam = 0;
942  pksrec.beamNo  = ibeam + 1;
943
944  //pointing/azel
945  MVPosition mvpos(antennaCols.position()(cAntId[0]));
946  MPosition mp(mvpos);
947  Quantum<Double> qt(time,"s");
948  MVEpoch mvt(qt);
949  MEpoch me(mvt);
950  MeasFrame frame(mp, me);
951  MDirection md;
952  pksrec.pCode = 0;
953  pksrec.rateAge = 0.0f;
954  pksrec.paRate = 0.0f;
955  if (cGetPointing) {
956    //cerr << "get pointing data ...." << endl;
957    ROScalarColumn<Int> pAntIdCol ;
958    ROScalarColumn<Double> psTimeCol ;
959    Table ptTable = cPKSMS.pointing() ;
960    MSPointing selPtTab( ptTable( ptTable.col("ANTENNA_ID") == cAntId[0] ) ) ;
961    pAntIdCol.attach( selPtTab, "ANTENNA_ID" ) ;
962    Vector<Int> antIds = pAntIdCol.getColumn() ;
963    psTimeCol.attach( selPtTab, "TIME" ) ;
964    Vector<Double> pTimes = psTimeCol.getColumn();
965    Bool doInterp = False ;
966    Int PtIdx=-1;
967    for (PtIdx = pTimes.nelements()-1; PtIdx >= 0; PtIdx--) {
968      if ( pTimes[PtIdx] == time ) {
969        break ;
970      }
971      else if ( pTimes[PtIdx] < time ) {
972        if ( PtIdx != pTimes.nelements()-1 ) {
973          doInterp = True ;
974        }
975        break ;
976      }
977    }
978    if ( PtIdx == -1 ) {
979      PtIdx = 0 ;
980    }
981    //cerr << "got index=" << PtIdx << endl;
982    Matrix<Double> pointingDir = cPointingCol(PtIdx);
983    ROMSPointingColumns PtCols( selPtTab ) ;
984    Vector<Double> pointingDirVec ;
985    if ( doInterp ) {
986      Double dt1 = time - pTimes[PtIdx] ;
987      Double dt2 = pTimes[PtIdx+1] - time ;
988      Vector<Double> dirVec1 = pointingDir.column(0) ;
989      Matrix<Double> pointingDir2 = cPointingCol(PtIdx+1) ;
990      Vector<Double> dirVec2 = pointingDir2.column(0) ;
991      pointingDirVec = (dt1*dirVec2+dt2*dirVec1)/(dt1+dt2) ;
992      Vector<MDirection> vmd1(1) ;
993      Vector<MDirection> vmd2(1) ;
994      PtCols.directionMeasCol().get(PtIdx,vmd1) ;
995      Vector<Double> angle1 = vmd1(0).getAngle().getValue("rad") ;
996      PtCols.directionMeasCol().get(PtIdx+1,vmd2) ;
997      Vector<Double> angle2 = vmd2(0).getAngle().getValue("rad") ;
998      Vector<Double> angle = (dt1*angle2+dt2*angle1)/(dt1+dt2) ;
999      Quantum< Vector<Double> > qangle( angle, "rad" ) ;
1000      String typeStr = vmd1(0).getRefString() ;
1001      //cerr << "vmd1.getRefString()=" << typeStr << endl ;
1002      MDirection::Types mdType ;
1003      MDirection::getType( mdType, typeStr ) ;
1004      //cerr << "mdType=" << mdType << endl ;
1005      md = MDirection( qangle, mdType ) ;
1006      //cerr << "md=" << md.getAngle().getValue("rad") << endl ;
1007    }
1008    else {
1009      pointingDirVec = pointingDir.column(0) ;
1010      Vector<MDirection> vmd(1);
1011      PtCols.directionMeasCol().get(PtIdx,vmd);
1012      md = vmd[0];
1013    }
1014    // put J2000 coordinates in "direction"
1015    if (cDirRef =="J2000") {
1016      pksrec.direction = pointingDirVec ;
1017    }
1018    else {
1019      pksrec.direction =
1020        MDirection::Convert(md, MDirection::Ref(MDirection::J2000,
1021                                                frame)
1022                            )().getAngle("rad").getValue();
1023     
1024    }
1025    uInt ncols = pointingDir.ncolumn();
1026    pksrec.scanRate.resize(2);
1027    if (ncols == 1) {
1028      pksrec.scanRate = 0.0f;
1029    } else {
1030      pksrec.scanRate(0) = pointingDir.column(1)(0);
1031      pksrec.scanRate(1) = pointingDir.column(1)(1);
1032    }
1033  }
1034  else {
1035  // Get direction from FIELD table
1036  // here, assume direction to be the field direction not pointing
1037    Matrix<Double> delayDir = cFieldDelayDirCol(fieldId);
1038    pksrec.direction = delayDir.column(0);
1039    uInt ncols = delayDir.ncolumn();
1040    pksrec.scanRate.resize(2);
1041    if (ncols == 1) {
1042      pksrec.scanRate = 0.0f;
1043    } else {
1044      pksrec.scanRate(0)  = delayDir.column(1)(0);
1045      pksrec.scanRate(1)  = delayDir.column(1)(1);
1046    }
1047  }
1048  // caluculate azimuth and elevation
1049  // first, get the reference frame
1050 /**
1051  MVPosition mvpos(antennaCols.position()(0));
1052  MPosition mp(mvpos);
1053  Quantum<Double> qt(time,"s");
1054  MVEpoch mvt(qt);
1055  MEpoch me(mvt);
1056  MeasFrame frame(mp, me);
1057  **/
1058  //
1059  ROMSFieldColumns fldCols(cPKSMS.field());
1060  Vector<MDirection> vmd(1);
1061  //MDirection md;
1062  fldCols.delayDirMeasCol().get(fieldId,vmd);
1063  md = vmd[0];
1064  //Vector<Double> dircheck = md.getAngle("rad").getValue();
1065  //cerr<<"dircheck="<<dircheck<<endl;
1066
1067  Vector<Double> azel =
1068        MDirection::Convert(md, MDirection::Ref(MDirection::AZEL,
1069                                                frame)
1070                            )().getAngle("rad").getValue();
1071  //cerr<<"azel="<<azel<<endl;
1072  pksrec.azimuth = azel[0];
1073  pksrec.elevation = azel[1];
1074
1075  // Get Tsys assuming that entries in the SYSCAL table match the main table.
1076  if (cHaveTsys) {
1077    Int nTsysColRow = cTsysCol.nrow();
1078    if (nTsysColRow != cNRow) {
1079      cHaveTsys=0;
1080    }
1081  }
1082  if (cHaveTsys) {
1083    cTsysCol.get(cIdx, pksrec.tsys, True);
1084  } else {
1085    Int numReceptor;
1086    cNumReceptorCol.get(0, numReceptor);
1087    pksrec.tsys.resize(numReceptor);
1088    pksrec.tsys = 1.0f;
1089  }
1090  cSigmaCol.get(cIdx, pksrec.sigma, True);
1091
1092  //get Tcal if available
1093  if (cHaveTcal) {
1094    Int nTcalColRow = cTcalCol.nrow();
1095    uInt nBeam = cBeams.nelements();
1096    uInt nIF = cIFs.nelements();
1097    uInt nrws = nBeam * nIF;
1098    if (nTcalColRow > 0) { 
1099    // find tcal match with the data with the data time stamp
1100      Double mjds = pksrec.mjd*(24*3600);
1101      Double dtcalTime;
1102      if ( pksrec.mjd > lastmjd || cIdx==0 ) {
1103        //Table tmptab = cSysCalTab(near(cSysCalTab.col("TIME"),mjds));
1104        tmptab = cSysCalTab(near(cSysCalTab.col("TIME"),mjds), nrws);
1105        //DEBUG
1106        //if (cIdx == 0) {
1107        //  cerr<<"inital table retrieved"<<endl;
1108        //}
1109       
1110      }
1111
1112      if (nBeam == 1) {
1113        tmptab2 = tmptab( tmptab.col("SPECTRAL_WINDOW_ID") == iIF, 1);
1114      } else {
1115        tmptab2 = tmptab( tmptab.col("SPECTRAL_WINDOW_ID") == iIF &&
1116                              tmptab.col("FEED_ID") == ibeam , 1);
1117      }
1118      //cerr<<"first subtab rows="<<tmptab.nrow()<<endl;
1119      int syscalrow = tmptab2.nrow();
1120      ROArrayColumn<Float> tcalCol(tmptab2, "TCAL");
1121      ROScalarColumn<Double> tcalTimeCol(tmptab2, "TIME");
1122      if (syscalrow==0) {
1123        os << LogIO::NORMAL
1124           <<"Cannot find any matching Tcal at/near the data timestamp."
1125           << " Set Tcal=0.0" << LogIO::POST ;
1126      } else {
1127        tcalCol.get(0, pksrec.tcal);
1128        tcalTimeCol.get(0,dtcalTime);
1129        pksrec.tcalTime = MVTime(dtcalTime/(24*3600)).string(MVTime::YMD);
1130        //DEBUG
1131        //cerr<<"cIdx:"<<cIdx<<" tcal="<<tcal<<" tcalTime="<<tcalTime<<endl;
1132        tmptab.markForDelete();
1133        tmptab2.markForDelete();
1134      }
1135    }
1136    lastmjd = pksrec.mjd;
1137  }
1138
1139  // Calibration factors (if available).
1140  pksrec.calFctr.resize(cNPol(iIF));
1141  if (cHaveCalFctr) {
1142    cCalFctrCol.get(cIdx, pksrec.calFctr);
1143  } else {
1144    pksrec.calFctr = 0.0f;
1145  }
1146
1147  // Baseline parameters (if available).
1148  if (cHaveBaseLin) {
1149    pksrec.baseLin.resize(2,cNPol(iIF));
1150    cBaseLinCol.get(cIdx, pksrec.baseLin);
1151
1152    pksrec.baseSub.resize(24,cNPol(iIF));
1153    cBaseSubCol.get(cIdx, pksrec.baseSub);
1154
1155  } else {
1156    pksrec.baseLin.resize(0,0);
1157    pksrec.baseSub.resize(0,0);
1158  }
1159
1160
1161  // Get spectral data.
1162  if (cGetSpectra) {
1163    Matrix<Float> tmpData;
1164    Matrix<Bool>  tmpFlag;
1165    if (cHaveDataCol) {
1166      Matrix<Complex> tmpCmplxData;
1167      Matrix<Float> tmpReData;
1168      Matrix<Float> tmpImData;
1169      //cerr<<"reading spectra..."<<endl;
1170      //# TODO - should have a flag to user to select DATA or CORRECTED_DATA
1171      //# currently just automatically determined, --- read CORRECTED one
1172      //# if the column exist.
1173      if (cHaveCorrectedDataCol) {
1174        cCorrectedDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True);
1175      } else {
1176        cDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True);
1177      }
1178      tmpReData = real(tmpCmplxData);
1179      tmpImData = imag(tmpCmplxData);
1180      tmpData = sqrt(tmpReData*tmpReData + tmpImData*tmpImData);
1181    } else {
1182      cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
1183    }
1184    cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
1185
1186    // Transpose spectra.
1187    Int nPol = tmpData.nrow();
1188    pksrec.spectra.resize(nChan, nPol);
1189    pksrec.flagged.resize(nChan, nPol);
1190    if (cEndChan(iIF) >= cStartChan(iIF)) {
1191      // Simple transposition.
1192      for (Int ipol = 0; ipol < nPol; ipol++) {
1193        for (Int ichan = 0; ichan < nChan; ichan++) {
1194          pksrec.spectra(ichan,ipol) = tmpData(ipol,ichan);
1195          pksrec.flagged(ichan,ipol) = tmpFlag(ipol,ichan);
1196        }
1197      }
1198
1199    } else {
1200      // Transpose with inversion.
1201      Int jchan = nChan - 1;
1202      for (Int ipol = 0; ipol < nPol; ipol++) {
1203        for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
1204          pksrec.spectra(ichan,ipol) = tmpData(ipol,jchan);
1205          pksrec.flagged(ichan,ipol) = tmpFlag(ipol,jchan);
1206        }
1207      }
1208    }
1209
1210    // Row-based flagging info. (True:1, False:0)
1211    pksrec.flagrow = (cFlagRowCol(cIdx) ? 1 : 0);
1212  }
1213
1214  // Get cross-polarization data.
1215  if (cGetXPol) {
1216    //cerr<<"cGetXPol="<<cGetXPol<<endl;
1217    //cerr<<"cHaveXCalFctr="<<cHaveXCalFctr<<endl;
1218
1219    if (cHaveXCalFctr) {
1220      cXCalFctrCol.get(cIdx, pksrec.xCalFctr);
1221    } else {
1222      pksrec.xCalFctr = Complex(0.0f, 0.0f);
1223    }
1224
1225    if(!cALMA) {
1226      cDataCol.get(cIdx, pksrec.xPol, True);
1227
1228      if (cEndChan(iIF) < cStartChan(iIF)) {
1229        Complex ctmp;
1230        Int jchan = nChan - 1;
1231        for (Int ichan = 0; ichan < nChan/2; ichan++, jchan--) {
1232          ctmp = pksrec.xPol(ichan);
1233          pksrec.xPol(ichan) = pksrec.xPol(jchan);
1234          pksrec.xPol(jchan) = ctmp;
1235        }
1236      }
1237    }
1238  }
1239  /**
1240  cerr<<"scanNo="<<scanNo<<endl;
1241  cerr<<"cycleNo="<<cycleNo<<endl;
1242  cerr<<"mjd="<<mjd<<endl;
1243  cerr<<"interval="<<interval<<endl;
1244  cerr<<"fieldName="<<fieldName<<endl;
1245  cerr<<"srcNmae="<<srcName<<endl;
1246  cerr<<"srcDir="<<srcDir<<endl;
1247  cerr<<"srcPM="<<srcPM<<endl;
1248  cerr<<"srcVel="<<srcVel<<endl;
1249  cerr<<"obsMode="<<obsMode<<endl;
1250  cerr<<"IFno="<<IFno<<endl;
1251  cerr<<"refFreq="<<refFreq<<endl;
1252  cerr<<"tcal="<<tcal<<endl;
1253  cerr<<"direction="<<direction<<endl;
1254  cerr<<"scanRate="<<scanRate<<endl;
1255  cerr<<"tsys="<<tsys<<endl;
1256  cerr<<"sigma="<<sigma<<endl;
1257  cerr<<"calFctr="<<calFctr<<endl;
1258  cerr<<"baseLin="<<baseLin<<endl;
1259  cerr<<"baseSub="<<baseSub<<endl;
1260  cerr<<"spectra="<<spectra.shape()<<endl;
1261  cerr<<"flagged="<<flagged.shape()<<endl;
1262  cerr<<"xCalFctr="<<xCalFctr<<endl;
1263  cerr<<"xPol="<<xPol<<endl;
1264  **/
1265  cIdx++;
1266
1267  return 0;
1268}
1269
1270//--------------------------------------------------------- PKSMS2reader::read
1271
1272// Read the next data record, just the basics.
1273
1274Int PKSMS2reader::read(
1275        Int           &IFno,
1276        Vector<Float> &tsys,
1277        Vector<Float> &calFctr,
1278        Matrix<Float> &baseLin,
1279        Matrix<Float> &baseSub,
1280        Matrix<Float> &spectra,
1281        Matrix<uChar> &flagged)
1282{
1283  if (!cMSopen) {
1284    return 1;
1285  }
1286
1287  // Check for EOF.
1288  if (cIdx >= cNRow) {
1289    return -1;
1290  }
1291
1292  // Find the next selected beam and IF.
1293  Int ibeam;
1294  Int iIF;
1295  Int iDataDesc;
1296  while (True) {
1297    ibeam = cBeamNoCol(cIdx);
1298    //iIF   = cDataDescIdCol(cIdx);
1299    iDataDesc   = cDataDescIdCol(cIdx);
1300    iIF   = cSpWinIdCol(iDataDesc);
1301    if (cBeams(ibeam) && cIFs(iIF)) {
1302      break;
1303    }
1304
1305    // Check for EOF.
1306    if (++cIdx >= cNRow) {
1307      return -1;
1308    }
1309  }
1310
1311  IFno = iIF + 1;
1312  // Get Tsys assuming that entries in the SYSCAL table match the main table.
1313  cTsysCol.get(cIdx, tsys, True);
1314
1315  // Calibration factors (if available).
1316  if (cHaveCalFctr) {
1317    cCalFctrCol.get(cIdx, calFctr, True);
1318  } else {
1319    calFctr.resize(cNPol(iIF));
1320    calFctr = 0.0f;
1321  }
1322
1323  // Baseline parameters (if available).
1324  if (cHaveBaseLin) {
1325    baseLin.resize(2,cNPol(iIF));
1326    cBaseLinCol.get(cIdx, baseLin);
1327
1328    baseSub.resize(24,cNPol(iIF));
1329    cBaseSubCol.get(cIdx, baseSub);
1330
1331  } else {
1332    baseLin.resize(0,0);
1333    baseSub.resize(0,0);
1334  }
1335
1336  if (cGetSpectra) {
1337    // Get spectral data.
1338    Matrix<Float> tmpData;
1339    Matrix<Bool>  tmpFlag;
1340    if (cHaveDataCol) {
1341      Matrix<Complex> tmpCmplxData;
1342      cDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True);
1343      tmpData = real(tmpCmplxData);
1344    } else {
1345      cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
1346    }
1347    cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
1348
1349    // Transpose spectra.
1350    Int nChan = tmpData.ncolumn();
1351    Int nPol  = tmpData.nrow();
1352    spectra.resize(nChan, nPol);
1353    flagged.resize(nChan, nPol);
1354    if (cEndChan(iIF) >= cStartChan(iIF)) {
1355      // Simple transposition.
1356      for (Int ipol = 0; ipol < nPol; ipol++) {
1357        for (Int ichan = 0; ichan < nChan; ichan++) {
1358          spectra(ichan,ipol) = tmpData(ipol,ichan);
1359          flagged(ichan,ipol) = tmpFlag(ipol,ichan);
1360        }
1361      }
1362
1363    } else {
1364      // Transpose with inversion.
1365      Int jchan = nChan - 1;
1366      for (Int ipol = 0; ipol < nPol; ipol++) {
1367        for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
1368          spectra(ichan,ipol) = tmpData(ipol,jchan);
1369          flagged(ichan,ipol) = tmpFlag(ipol,jchan);
1370        }
1371      }
1372    }
1373  }
1374
1375  cIdx++;
1376
1377  return 0;
1378}
1379
1380//-------------------------------------------------------- PKSMS2reader::close
1381
1382// Close the MS.
1383
1384void PKSMS2reader::close()
1385{
1386  cPKSMS = MeasurementSet();
1387  cMSopen = False;
1388}
1389
1390//-------------------------------------------------------- PKSMS2reader::splitAntenanSelectionString
1391
1392// split antenna selection string
1393// delimiter is ','
1394
1395Vector<String> PKSMS2reader::splitAntennaSelectionString( const String s )
1396{
1397  Char delim = ',' ;
1398  Int n = s.freq( delim ) + 1 ;
1399  Vector<String> antlist ;
1400  string *sl = new string[n] ;
1401  Int numSubstr = split( s, sl, n, "," );
1402  antlist.resize( numSubstr ) ;
1403  for ( Int i = 0 ; i < numSubstr ; i++ ) {
1404    antlist[i] = String( sl[i] ) ;
1405    antlist[i].trim() ;
1406  }
1407  //cerr << "antlist = " << antlist << endl ;
1408  delete [] sl;
1409  return antlist ;
1410}
1411
1412//-------------------------------------------------------- PKSMS2reader::setupAntennaList
1413
1414// Fill cAntenna and cAntId
1415
1416void PKSMS2reader::setupAntennaList( const String s )
1417{
1418  LogIO os( LogOrigin( "PKSMS2reader", "setupAntennaList()", WHERE ) ) ;
1419  //cerr << "antenna specification: " << s << endl ;
1420  ROMSAntennaColumns antennaCols(cPKSMS.antenna());
1421  ROScalarColumn<String> antNames = antennaCols.name();
1422  Int nrow = antNames.nrow() ;
1423  Vector<String> antlist = splitAntennaSelectionString( s ) ;
1424  Int len = antlist.size() ;
1425  Vector<Int> AntId( len ) ;
1426  Regex re( "[0-9]+" ) ;
1427  for ( Int i = 0 ; i < len ; i++ ) {
1428    if ( antlist[i].matches( re ) ) {
1429      AntId[i] = atoi( antlist[i].c_str() ) ;
1430      if ( AntId[i] >= nrow ) {
1431        os << LogIO::SEVERE << "Antenna index out of range: " << AntId[i] << LogIO::EXCEPTION ;
1432      }
1433    }
1434    else {
1435      AntId[i] = -1 ;
1436      for ( uInt j = 0 ; j < antNames.nrow() ; j++ ) {
1437        if ( antlist[i] == antNames(j) ) {
1438          AntId[i] = j ;
1439          break ;
1440        }
1441      }
1442      if ( AntId[i] == -1 ) {
1443        os << LogIO::SEVERE << "Specified antenna name not found: " << antlist[i] << LogIO::EXCEPTION ;
1444      }
1445    }
1446  }
1447  //cerr << "AntId = " << AntId << endl ;
1448  vector<Int> uniqId ;
1449  uniqId.push_back( AntId(0) ) ;
1450  for ( uInt i = 1 ; i < AntId.size() ; i++ ) {
1451    if ( count(uniqId.begin(),uniqId.end(),AntId[i]) == 0 ) {
1452      uniqId.push_back( AntId[i] ) ;
1453    }
1454  }
1455  Vector<Int> newAntId( uniqId ) ;
1456  cAntId.assign( newAntId ) ;
1457  //cerr << "cAntId = " << cAntId << endl ;
1458}
Note: See TracBrowser for help on using the repository browser.