source: branches/alma/external-alma/atnf/PKSIO/PKSMS2reader.cc@ 2518

Last change on this file since 2518 was 1818, checked in by Kana Sugimoto, 14 years ago

New Development: Yes

JIRA Issue: No (merge)

Ready for Test: Yes

Interface Changes: No

Description:

Merged changes -r1774:1817 in newfiller branch to alma branch


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