source: branches/casa-release-4_3-test/external-alma/atnf/PKSIO/PKSMS2reader.cc@ 3037

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

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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