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

Last change on this file since 2474 was 1880, checked in by Takeshi Nakazato, 14 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...

An implementation of asapmath.splitant() is changed since
the new filler doesn't support antenna parameter at the
moment.

The default behavior of MS filler also changed.
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 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[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 return antlist ;
1409}
1410
1411//-------------------------------------------------------- PKSMS2reader::setupAntennaList
1412
1413// Fill cAntenna and cAntId
1414
1415void PKSMS2reader::setupAntennaList( const String s )
1416{
1417 LogIO os( LogOrigin( "PKSMS2reader", "setupAntennaList()", WHERE ) ) ;
1418 //cerr << "antenna specification: " << s << endl ;
1419 ROMSAntennaColumns antennaCols(cPKSMS.antenna());
1420 ROScalarColumn<String> antNames = antennaCols.name();
1421 Int nrow = antNames.nrow() ;
1422 Vector<String> antlist = splitAntennaSelectionString( s ) ;
1423 Int len = antlist.size() ;
1424 Vector<Int> AntId( len ) ;
1425 Regex re( "[0-9]+" ) ;
1426 for ( Int i = 0 ; i < len ; i++ ) {
1427 if ( antlist[i].matches( re ) ) {
1428 AntId[i] = atoi( antlist[i].c_str() ) ;
1429 if ( AntId[i] >= nrow ) {
1430 os << LogIO::SEVERE << "Antenna index out of range: " << AntId[i] << LogIO::EXCEPTION ;
1431 }
1432 }
1433 else {
1434 AntId[i] = -1 ;
1435 for ( uInt j = 0 ; j < antNames.nrow() ; j++ ) {
1436 if ( antlist[i] == antNames(j) ) {
1437 AntId[i] = j ;
1438 break ;
1439 }
1440 }
1441 if ( AntId[i] == -1 ) {
1442 os << LogIO::SEVERE << "Specified antenna name not found: " << antlist[i] << LogIO::EXCEPTION ;
1443 }
1444 }
1445 }
1446 //cerr << "AntId = " << AntId << endl ;
1447 vector<Int> uniqId ;
1448 uniqId.push_back( AntId(0) ) ;
1449 for ( uInt i = 1 ; i < AntId.size() ; i++ ) {
1450 if ( count(uniqId.begin(),uniqId.end(),AntId[i]) == 0 ) {
1451 uniqId.push_back( AntId[i] ) ;
1452 }
1453 }
1454 Vector<Int> newAntId( uniqId ) ;
1455 cAntId.assign( newAntId ) ;
1456 //cerr << "cAntId = " << cAntId << endl ;
1457}
Note: See TracBrowser for help on using the repository browser.