source: branches/newfiller/external/atnf/PKSIO/PKSMS2reader.cc@ 1792

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

New Development: Yes

JIRA Issue: Yes (CAS-2211)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: ASAP 3.0.0 interface changes

Test Programs:

Put in Release Notes: Yes

Module(s): all the CASA sd tools and tasks are affected.

Description: Merged ATNF-ASAP 3.0.0 developments to CASA (alma) branch.

Note you also need to update casa/code/atnf.
File size: 42.4 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 pksrec.nchan = nChan;
873 if (nChan == 1) {
874 //pksrec.freqInc = chanFreq(0);
875 pksrec.freqInc = cTotBWCol(iIF);
876 pksrec.refFreq = chanFreq(0);
877 pksrec.restFreq.resize(1);
878 pksrec.restFreq[0] = 0.0f;
879 } else {
880
881 if (cStartChan(iIF) <= cEndChan(iIF)) {
882 pksrec.freqInc = chanFreq(1) - chanFreq(0);
883 } else {
884 pksrec.freqInc = chanFreq(0) - chanFreq(1);
885 }
886
887 pksrec.refFreq = chanFreq(cRefChan(iIF)-1);
888
889 Bool HaveSrcRestFreq= cSrcRestFrqCol.isDefined(srcId);
890 if (HaveSrcRestFreq) {
891 //restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0));
892 //restFreq = cSrcRestFrqCol(srcId);
893 pksrec.restFreq = cSrcRestFrqCol(srcId);
894 } else {
895 pksrec.restFreq.resize(1);
896 pksrec.restFreq[0] = 0.0f;
897 }
898 }
899 //pksrec.bandwidth = abs(pksrec.freqInc * nChan);
900 pksrec.bandwidth = abs(cTotBWCol(0));
901
902 pksrec.tcal.resize(cNPol(iIF));
903 pksrec.tcal = 0.0f;
904 pksrec.tcalTime = "";
905// pksrec.azimuth = 0.0f;
906// pksrec.elevation = 0.0f;
907 pksrec.parAngle = 0.0f;
908
909 pksrec.focusAxi = 0.0f;
910 pksrec.focusTan = 0.0f;
911 pksrec.focusRot = 0.0f;
912
913 // Find the appropriate entry in the WEATHER subtable.
914 //Bool cHaveStateTab=Table::isReadable(cPKSMS.stateTableName());
915 Bool cHaveWeatherTab = Table::isReadable(cPKSMS.weatherTableName());
916 Int weatherIdx=-1;
917 if (cHaveWeatherTab) {
918 Vector<Double> wTimes = cWeatherTimeCol.getColumn();
919 for (weatherIdx = wTimes.nelements()-1; weatherIdx >= 0; weatherIdx--) {
920 if (cWeatherTimeCol(weatherIdx) <= time) {
921 break;
922 }
923 }
924 }
925
926 if (weatherIdx < 0 || !cHaveWeatherTab) {
927 // No appropriate WEATHER entry.
928 pksrec.temperature = 0.0f;
929 pksrec.pressure = 0.0f;
930 pksrec.humidity = 0.0f;
931 } else {
932 pksrec.temperature = cTemperatureCol(weatherIdx);
933 pksrec.pressure = cPressureCol(weatherIdx);
934 pksrec.humidity = cHumidityCol(weatherIdx);
935 }
936
937 pksrec.windSpeed = 0.0f;
938 pksrec.windAz = 0.0f;
939
940 pksrec.refBeam = 0;
941 pksrec.beamNo = ibeam + 1;
942
943 //pointing/azel
944 MVPosition mvpos(antennaCols.position()(cAntId[0]));
945 MPosition mp(mvpos);
946 Quantum<Double> qt(time,"s");
947 MVEpoch mvt(qt);
948 MEpoch me(mvt);
949 MeasFrame frame(mp, me);
950 MDirection md;
951 pksrec.pCode = 0;
952 pksrec.rateAge = 0.0f;
953 pksrec.paRate = 0.0f;
954 if (cGetPointing) {
955 //cerr << "get pointing data ...." << endl;
956 ROScalarColumn<Int> pAntIdCol ;
957 ROScalarColumn<Double> psTimeCol ;
958 Table ptTable = cPKSMS.pointing() ;
959 MSPointing selPtTab( ptTable( ptTable.col("ANTENNA_ID") == cAntId[0] ) ) ;
960 pAntIdCol.attach( selPtTab, "ANTENNA_ID" ) ;
961 Vector<Int> antIds = pAntIdCol.getColumn() ;
962 psTimeCol.attach( selPtTab, "TIME" ) ;
963 Vector<Double> pTimes = psTimeCol.getColumn();
964 Bool doInterp = False ;
965 Int PtIdx=-1;
966 for (PtIdx = pTimes.nelements()-1; PtIdx >= 0; PtIdx--) {
967 if ( pTimes[PtIdx] == time ) {
968 break ;
969 }
970 else if ( pTimes[PtIdx] < time ) {
971 if ( PtIdx != pTimes.nelements()-1 ) {
972 doInterp = True ;
973 }
974 break ;
975 }
976 }
977 if ( PtIdx == -1 ) {
978 PtIdx = 0 ;
979 }
980 //cerr << "got index=" << PtIdx << endl;
981 Matrix<Double> pointingDir = cPointingCol(PtIdx);
982 ROMSPointingColumns PtCols( selPtTab ) ;
983 Vector<Double> pointingDirVec ;
984 if ( doInterp ) {
985 Double dt1 = time - pTimes[PtIdx] ;
986 Double dt2 = pTimes[PtIdx+1] - time ;
987 Vector<Double> dirVec1 = pointingDir.column(0) ;
988 Matrix<Double> pointingDir2 = cPointingCol(PtIdx+1) ;
989 Vector<Double> dirVec2 = pointingDir2.column(0) ;
990 pointingDirVec = (dt1*dirVec2+dt2*dirVec1)/(dt1+dt2) ;
991 Vector<MDirection> vmd1(1) ;
992 Vector<MDirection> vmd2(1) ;
993 PtCols.directionMeasCol().get(PtIdx,vmd1) ;
994 Vector<Double> angle1 = vmd1(0).getAngle().getValue("rad") ;
995 PtCols.directionMeasCol().get(PtIdx+1,vmd2) ;
996 Vector<Double> angle2 = vmd2(0).getAngle().getValue("rad") ;
997 Vector<Double> angle = (dt1*angle2+dt2*angle1)/(dt1+dt2) ;
998 Quantum< Vector<Double> > qangle( angle, "rad" ) ;
999 String typeStr = vmd1(0).getRefString() ;
1000 //cerr << "vmd1.getRefString()=" << typeStr << endl ;
1001 MDirection::Types mdType ;
1002 MDirection::getType( mdType, typeStr ) ;
1003 //cerr << "mdType=" << mdType << endl ;
1004 md = MDirection( qangle, mdType ) ;
1005 //cerr << "md=" << md.getAngle().getValue("rad") << endl ;
1006 }
1007 else {
1008 pointingDirVec = pointingDir.column(0) ;
1009 Vector<MDirection> vmd(1);
1010 PtCols.directionMeasCol().get(PtIdx,vmd);
1011 md = vmd[0];
1012 }
1013 // put J2000 coordinates in "direction"
1014 if (cDirRef =="J2000") {
1015 pksrec.direction = pointingDirVec ;
1016 }
1017 else {
1018 pksrec.direction =
1019 MDirection::Convert(md, MDirection::Ref(MDirection::J2000,
1020 frame)
1021 )().getAngle("rad").getValue();
1022
1023 }
1024 uInt ncols = pointingDir.ncolumn();
1025 pksrec.scanRate.resize(2);
1026 if (ncols == 1) {
1027 pksrec.scanRate = 0.0f;
1028 } else {
1029 pksrec.scanRate(0) = pointingDir.column(1)(0);
1030 pksrec.scanRate(1) = pointingDir.column(1)(1);
1031 }
1032 }
1033 else {
1034 // Get direction from FIELD table
1035 // here, assume direction to be the field direction not pointing
1036 Matrix<Double> delayDir = cFieldDelayDirCol(fieldId);
1037 pksrec.direction = delayDir.column(0);
1038 uInt ncols = delayDir.ncolumn();
1039 pksrec.scanRate.resize(2);
1040 if (ncols == 1) {
1041 pksrec.scanRate = 0.0f;
1042 } else {
1043 pksrec.scanRate(0) = delayDir.column(1)(0);
1044 pksrec.scanRate(1) = delayDir.column(1)(1);
1045 }
1046 }
1047 // caluculate azimuth and elevation
1048 // first, get the reference frame
1049 /**
1050 MVPosition mvpos(antennaCols.position()(0));
1051 MPosition mp(mvpos);
1052 Quantum<Double> qt(time,"s");
1053 MVEpoch mvt(qt);
1054 MEpoch me(mvt);
1055 MeasFrame frame(mp, me);
1056 **/
1057 //
1058 ROMSFieldColumns fldCols(cPKSMS.field());
1059 Vector<MDirection> vmd(1);
1060 //MDirection md;
1061 fldCols.delayDirMeasCol().get(fieldId,vmd);
1062 md = vmd[0];
1063 //Vector<Double> dircheck = md.getAngle("rad").getValue();
1064 //cerr<<"dircheck="<<dircheck<<endl;
1065
1066 Vector<Double> azel =
1067 MDirection::Convert(md, MDirection::Ref(MDirection::AZEL,
1068 frame)
1069 )().getAngle("rad").getValue();
1070 //cerr<<"azel="<<azel<<endl;
1071 pksrec.azimuth = azel[0];
1072 pksrec.elevation = azel[1];
1073
1074 // Get Tsys assuming that entries in the SYSCAL table match the main table.
1075 if (cHaveTsys) {
1076 Int nTsysColRow = cTsysCol.nrow();
1077 if (nTsysColRow != cNRow) {
1078 cHaveTsys=0;
1079 }
1080 }
1081 if (cHaveTsys) {
1082 cTsysCol.get(cIdx, pksrec.tsys, True);
1083 } else {
1084 Int numReceptor;
1085 cNumReceptorCol.get(0, numReceptor);
1086 pksrec.tsys.resize(numReceptor);
1087 pksrec.tsys = 1.0f;
1088 }
1089 cSigmaCol.get(cIdx, pksrec.sigma, True);
1090
1091 //get Tcal if available
1092 if (cHaveTcal) {
1093 Int nTcalColRow = cTcalCol.nrow();
1094 uInt nBeam = cBeams.nelements();
1095 uInt nIF = cIFs.nelements();
1096 uInt nrws = nBeam * nIF;
1097 if (nTcalColRow > 0) {
1098 // find tcal match with the data with the data time stamp
1099 Double mjds = pksrec.mjd*(24*3600);
1100 Double dtcalTime;
1101 if ( pksrec.mjd > lastmjd || cIdx==0 ) {
1102 //Table tmptab = cSysCalTab(near(cSysCalTab.col("TIME"),mjds));
1103 tmptab = cSysCalTab(near(cSysCalTab.col("TIME"),mjds), nrws);
1104 //DEBUG
1105 //if (cIdx == 0) {
1106 // cerr<<"inital table retrieved"<<endl;
1107 //}
1108
1109 }
1110
1111 if (nBeam == 1) {
1112 tmptab2 = tmptab( tmptab.col("SPECTRAL_WINDOW_ID") == iIF, 1);
1113 } else {
1114 tmptab2 = tmptab( tmptab.col("SPECTRAL_WINDOW_ID") == iIF &&
1115 tmptab.col("FEED_ID") == ibeam , 1);
1116 }
1117 //cerr<<"first subtab rows="<<tmptab.nrow()<<endl;
1118 int syscalrow = tmptab2.nrow();
1119 ROArrayColumn<Float> tcalCol(tmptab2, "TCAL");
1120 ROScalarColumn<Double> tcalTimeCol(tmptab2, "TIME");
1121 if (syscalrow==0) {
1122 os << LogIO::NORMAL
1123 <<"Cannot find any matching Tcal at/near the data timestamp."
1124 << " Set Tcal=0.0" << LogIO::POST ;
1125 } else {
1126 tcalCol.get(0, pksrec.tcal);
1127 tcalTimeCol.get(0,dtcalTime);
1128 pksrec.tcalTime = MVTime(dtcalTime/(24*3600)).string(MVTime::YMD);
1129 //DEBUG
1130 //cerr<<"cIdx:"<<cIdx<<" tcal="<<tcal<<" tcalTime="<<tcalTime<<endl;
1131 tmptab.markForDelete();
1132 tmptab2.markForDelete();
1133 }
1134 }
1135 lastmjd = pksrec.mjd;
1136 }
1137
1138 // Calibration factors (if available).
1139 pksrec.calFctr.resize(cNPol(iIF));
1140 if (cHaveCalFctr) {
1141 cCalFctrCol.get(cIdx, pksrec.calFctr);
1142 } else {
1143 pksrec.calFctr = 0.0f;
1144 }
1145
1146 // Baseline parameters (if available).
1147 if (cHaveBaseLin) {
1148 pksrec.baseLin.resize(2,cNPol(iIF));
1149 cBaseLinCol.get(cIdx, pksrec.baseLin);
1150
1151 pksrec.baseSub.resize(24,cNPol(iIF));
1152 cBaseSubCol.get(cIdx, pksrec.baseSub);
1153
1154 } else {
1155 pksrec.baseLin.resize(0,0);
1156 pksrec.baseSub.resize(0,0);
1157 }
1158
1159
1160 // Get spectral data.
1161 if (cGetSpectra) {
1162 Matrix<Float> tmpData;
1163 Matrix<Bool> tmpFlag;
1164 if (cHaveDataCol) {
1165 Matrix<Complex> tmpCmplxData;
1166 Matrix<Float> tmpReData;
1167 Matrix<Float> tmpImData;
1168 //cerr<<"reading spectra..."<<endl;
1169 //# TODO - should have a flag to user to select DATA or CORRECTED_DATA
1170 //# currently just automatically determined, --- read CORRECTED one
1171 //# if the column exist.
1172 if (cHaveCorrectedDataCol) {
1173 cCorrectedDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True);
1174 } else {
1175 cDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True);
1176 }
1177 tmpReData = real(tmpCmplxData);
1178 tmpImData = imag(tmpCmplxData);
1179 tmpData = sqrt(tmpReData*tmpReData + tmpImData*tmpImData);
1180 } else {
1181 cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
1182 }
1183 cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
1184
1185 // Transpose spectra.
1186 Int nPol = tmpData.nrow();
1187 pksrec.spectra.resize(nChan, nPol);
1188 pksrec.flagged.resize(nChan, nPol);
1189 if (cEndChan(iIF) >= cStartChan(iIF)) {
1190 // Simple transposition.
1191 for (Int ipol = 0; ipol < nPol; ipol++) {
1192 for (Int ichan = 0; ichan < nChan; ichan++) {
1193 pksrec.spectra(ichan,ipol) = tmpData(ipol,ichan);
1194 pksrec.flagged(ichan,ipol) = tmpFlag(ipol,ichan);
1195 }
1196 }
1197
1198 } else {
1199 // Transpose with inversion.
1200 Int jchan = nChan - 1;
1201 for (Int ipol = 0; ipol < nPol; ipol++) {
1202 for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
1203 pksrec.spectra(ichan,ipol) = tmpData(ipol,jchan);
1204 pksrec.flagged(ichan,ipol) = tmpFlag(ipol,jchan);
1205 }
1206 }
1207 }
1208
1209 // Row-based flagging info. (True:1, False:0)
1210 pksrec.flagrow = (cFlagRowCol(cIdx) ? 1 : 0);
1211 }
1212
1213 // Get cross-polarization data.
1214 if (cGetXPol) {
1215 //cerr<<"cGetXPol="<<cGetXPol<<endl;
1216 //cerr<<"cHaveXCalFctr="<<cHaveXCalFctr<<endl;
1217
1218 if (cHaveXCalFctr) {
1219 cXCalFctrCol.get(cIdx, pksrec.xCalFctr);
1220 } else {
1221 pksrec.xCalFctr = Complex(0.0f, 0.0f);
1222 }
1223
1224 if(!cALMA) {
1225 cDataCol.get(cIdx, pksrec.xPol, True);
1226
1227 if (cEndChan(iIF) < cStartChan(iIF)) {
1228 Complex ctmp;
1229 Int jchan = nChan - 1;
1230 for (Int ichan = 0; ichan < nChan/2; ichan++, jchan--) {
1231 ctmp = pksrec.xPol(ichan);
1232 pksrec.xPol(ichan) = pksrec.xPol(jchan);
1233 pksrec.xPol(jchan) = ctmp;
1234 }
1235 }
1236 }
1237 }
1238 /**
1239 cerr<<"scanNo="<<scanNo<<endl;
1240 cerr<<"cycleNo="<<cycleNo<<endl;
1241 cerr<<"mjd="<<mjd<<endl;
1242 cerr<<"interval="<<interval<<endl;
1243 cerr<<"fieldName="<<fieldName<<endl;
1244 cerr<<"srcNmae="<<srcName<<endl;
1245 cerr<<"srcDir="<<srcDir<<endl;
1246 cerr<<"srcPM="<<srcPM<<endl;
1247 cerr<<"srcVel="<<srcVel<<endl;
1248 cerr<<"obsMode="<<obsMode<<endl;
1249 cerr<<"IFno="<<IFno<<endl;
1250 cerr<<"refFreq="<<refFreq<<endl;
1251 cerr<<"tcal="<<tcal<<endl;
1252 cerr<<"direction="<<direction<<endl;
1253 cerr<<"scanRate="<<scanRate<<endl;
1254 cerr<<"tsys="<<tsys<<endl;
1255 cerr<<"sigma="<<sigma<<endl;
1256 cerr<<"calFctr="<<calFctr<<endl;
1257 cerr<<"baseLin="<<baseLin<<endl;
1258 cerr<<"baseSub="<<baseSub<<endl;
1259 cerr<<"spectra="<<spectra.shape()<<endl;
1260 cerr<<"flagged="<<flagged.shape()<<endl;
1261 cerr<<"xCalFctr="<<xCalFctr<<endl;
1262 cerr<<"xPol="<<xPol<<endl;
1263 **/
1264 cIdx++;
1265
1266 return 0;
1267}
1268
1269//--------------------------------------------------------- PKSMS2reader::read
1270
1271// Read the next data record, just the basics.
1272
1273Int PKSMS2reader::read(
1274 Int &IFno,
1275 Vector<Float> &tsys,
1276 Vector<Float> &calFctr,
1277 Matrix<Float> &baseLin,
1278 Matrix<Float> &baseSub,
1279 Matrix<Float> &spectra,
1280 Matrix<uChar> &flagged)
1281{
1282 if (!cMSopen) {
1283 return 1;
1284 }
1285
1286 // Check for EOF.
1287 if (cIdx >= cNRow) {
1288 return -1;
1289 }
1290
1291 // Find the next selected beam and IF.
1292 Int ibeam;
1293 Int iIF;
1294 Int iDataDesc;
1295 while (True) {
1296 ibeam = cBeamNoCol(cIdx);
1297 //iIF = cDataDescIdCol(cIdx);
1298 iDataDesc = cDataDescIdCol(cIdx);
1299 iIF = cSpWinIdCol(iDataDesc);
1300 if (cBeams(ibeam) && cIFs(iIF)) {
1301 break;
1302 }
1303
1304 // Check for EOF.
1305 if (++cIdx >= cNRow) {
1306 return -1;
1307 }
1308 }
1309
1310 IFno = iIF + 1;
1311 // Get Tsys assuming that entries in the SYSCAL table match the main table.
1312 cTsysCol.get(cIdx, tsys, True);
1313
1314 // Calibration factors (if available).
1315 if (cHaveCalFctr) {
1316 cCalFctrCol.get(cIdx, calFctr, True);
1317 } else {
1318 calFctr.resize(cNPol(iIF));
1319 calFctr = 0.0f;
1320 }
1321
1322 // Baseline parameters (if available).
1323 if (cHaveBaseLin) {
1324 baseLin.resize(2,cNPol(iIF));
1325 cBaseLinCol.get(cIdx, baseLin);
1326
1327 baseSub.resize(24,cNPol(iIF));
1328 cBaseSubCol.get(cIdx, baseSub);
1329
1330 } else {
1331 baseLin.resize(0,0);
1332 baseSub.resize(0,0);
1333 }
1334
1335 if (cGetSpectra) {
1336 // Get spectral data.
1337 Matrix<Float> tmpData;
1338 Matrix<Bool> tmpFlag;
1339 if (cHaveDataCol) {
1340 Matrix<Complex> tmpCmplxData;
1341 cDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True);
1342 tmpData = real(tmpCmplxData);
1343 } else {
1344 cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
1345 }
1346 cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
1347
1348 // Transpose spectra.
1349 Int nChan = tmpData.ncolumn();
1350 Int nPol = tmpData.nrow();
1351 spectra.resize(nChan, nPol);
1352 flagged.resize(nChan, nPol);
1353 if (cEndChan(iIF) >= cStartChan(iIF)) {
1354 // Simple transposition.
1355 for (Int ipol = 0; ipol < nPol; ipol++) {
1356 for (Int ichan = 0; ichan < nChan; ichan++) {
1357 spectra(ichan,ipol) = tmpData(ipol,ichan);
1358 flagged(ichan,ipol) = tmpFlag(ipol,ichan);
1359 }
1360 }
1361
1362 } else {
1363 // Transpose with inversion.
1364 Int jchan = nChan - 1;
1365 for (Int ipol = 0; ipol < nPol; ipol++) {
1366 for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
1367 spectra(ichan,ipol) = tmpData(ipol,jchan);
1368 flagged(ichan,ipol) = tmpFlag(ipol,jchan);
1369 }
1370 }
1371 }
1372 }
1373
1374 cIdx++;
1375
1376 return 0;
1377}
1378
1379//-------------------------------------------------------- PKSMS2reader::close
1380
1381// Close the MS.
1382
1383void PKSMS2reader::close()
1384{
1385 cPKSMS = MeasurementSet();
1386 cMSopen = False;
1387}
1388
1389//-------------------------------------------------------- PKSMS2reader::splitAntenanSelectionString
1390
1391// split antenna selection string
1392// delimiter is ','
1393
1394Vector<String> PKSMS2reader::splitAntennaSelectionString( const String s )
1395{
1396 Char delim = ',' ;
1397 Int n = s.freq( delim ) + 1 ;
1398 Vector<String> antlist ;
1399 string sl[n] ;
1400 Int numSubstr = split( s, sl, n, "," );
1401 antlist.resize( numSubstr ) ;
1402 for ( Int i = 0 ; i < numSubstr ; i++ ) {
1403 antlist[i] = String( sl[i] ) ;
1404 antlist[i].trim() ;
1405 }
1406 //cerr << "antlist = " << antlist << endl ;
1407 return antlist ;
1408}
1409
1410//-------------------------------------------------------- PKSMS2reader::setupAntennaList
1411
1412// Fill cAntenna and cAntId
1413
1414void PKSMS2reader::setupAntennaList( const String s )
1415{
1416 LogIO os( LogOrigin( "PKSMS2reader", "setupAntennaList()", WHERE ) ) ;
1417 //cerr << "antenna specification: " << s << endl ;
1418 ROMSAntennaColumns antennaCols(cPKSMS.antenna());
1419 ROScalarColumn<String> antNames = antennaCols.name();
1420 Int nrow = antNames.nrow() ;
1421 Vector<String> antlist = splitAntennaSelectionString( s ) ;
1422 Int len = antlist.size() ;
1423 Vector<Int> AntId( len ) ;
1424 Regex re( "[0-9]+" ) ;
1425 for ( Int i = 0 ; i < len ; i++ ) {
1426 if ( antlist[i].matches( re ) ) {
1427 AntId[i] = atoi( antlist[i].c_str() ) ;
1428 if ( AntId[i] >= nrow ) {
1429 os << LogIO::SEVERE << "Antenna index out of range: " << AntId[i] << LogIO::EXCEPTION ;
1430 }
1431 }
1432 else {
1433 AntId[i] = -1 ;
1434 for ( uInt j = 0 ; j < antNames.nrow() ; j++ ) {
1435 if ( antlist[i] == antNames(j) ) {
1436 AntId[i] = j ;
1437 break ;
1438 }
1439 }
1440 if ( AntId[i] == -1 ) {
1441 os << LogIO::SEVERE << "Specified antenna name not found: " << antlist[i] << LogIO::EXCEPTION ;
1442 }
1443 }
1444 }
1445 //cerr << "AntId = " << AntId << endl ;
1446 vector<Int> uniqId ;
1447 uniqId.push_back( AntId(0) ) ;
1448 for ( uInt i = 1 ; i < AntId.size() ; i++ ) {
1449 if ( count(uniqId.begin(),uniqId.end(),AntId[i]) == 0 ) {
1450 uniqId.push_back( AntId[i] ) ;
1451 }
1452 }
1453 Vector<Int> newAntId( uniqId ) ;
1454 cAntId.assign( newAntId ) ;
1455 //cerr << "cAntId = " << cAntId << endl ;
1456}
Note: See TracBrowser for help on using the repository browser.