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

Last change on this file since 1488 was 1453, checked in by TakTsutsumi, 16 years ago

New Development: No

JIRA Issue: No

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: many

Test Programs: sd.scantable(), sd.scantable.save()

Put in Release Notes: N/A

Description: copied from current casapy code tree


File size: 31.9 KB
Line 
1//#---------------------------------------------------------------------------
2//# PKSMS2reader.cc: Class to read Parkes Multibeam data from a v2 MS.
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2000-2006
5//# Associated Universities, Inc. Washington DC, USA.
6//#
7//# This library is free software; you can redistribute it and/or modify it
8//# under the terms of the GNU Library General Public License as published by
9//# the Free Software Foundation; either version 2 of the License, or (at your
10//# option) any later version.
11//#
12//# This library is distributed in the hope that it will be useful, but
13//# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14//# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
15//# License for more details.
16//#
17//# You should have received a copy of the GNU Library General Public License
18//# along with this library; if not, write to the Free Software Foundation,
19//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20//#
21//# Correspondence concerning AIPS++ should be addressed as follows:
22//# Internet email: aips2-request@nrao.edu.
23//# Postal address: AIPS++ Project Office
24//# National Radio Astronomy Observatory
25//# 520 Edgemont Road
26//# Charlottesville, VA 22903-2475 USA
27//#
28//# $Id$
29//#---------------------------------------------------------------------------
30//# Original: 2000/08/03, Mark Calabretta, ATNF
31//#---------------------------------------------------------------------------
32
33
34// AIPS++ includes.
35#include <casa/stdio.h>
36#include <casa/Arrays/ArrayMath.h>
37#include <casa/Arrays/Slice.h>
38#include <ms/MeasurementSets/MSColumns.h>
39#include <tables/Tables.h>
40#include <casa/Quanta/MVTime.h>
41#include <casa/Quanta/MVAngle.h>
42#include <casa/BasicMath/Math.h>
43#include <measures/Measures/MeasConvert.h>
44#include <measures/Measures/MEpoch.h>
45#include <measures/Measures/MeasRef.h>
46
47
48// Parkes includes.
49#include <atnf/pks/pks_maths.h>
50#include <atnf/PKSIO/PKSMS2reader.h>
51
52
53//------------------------------------------------- PKSMS2reader::PKSMS2reader
54
55// Default constructor.
56
57PKSMS2reader::PKSMS2reader()
58{
59 cMSopen = False;
60}
61
62//------------------------------------------------ PKSMS2reader::~PKSMS2reader
63
64PKSMS2reader::~PKSMS2reader()
65{
66 close();
67}
68
69//--------------------------------------------------------- PKSMS2reader::open
70
71// Open the MS for reading.
72
73Int PKSMS2reader::open(
74 const String msName,
75 Vector<Bool> &beams,
76 Vector<Bool> &IFs,
77 Vector<uInt> &nChan,
78 Vector<uInt> &nPol,
79 Vector<Bool> &haveXPol,
80 Bool &haveBase,
81 Bool &haveSpectra)
82{
83 // Check that MS is readable.
84 if (!MS::isReadable(msName)) {
85 return 1;
86 }
87
88 if (cMSopen) {
89 close();
90 }
91
92 cPKSMS = MeasurementSet(msName);
93 // taql access to the syscal table
94 cHaveSysCal = False;
95 if (cHaveSysCal=Table::isReadable(cPKSMS.sysCalTableName())) {
96 cSysCalTab = Table(cPKSMS.sysCalTableName());
97 }
98
99 cIdx = 0;
100 lastmjd = 0.0;
101 cNRow = cPKSMS.nrow();
102 cMSopen = True;
103
104 // Lock the table for read access.
105 cPKSMS.lock(False);
106
107 // Main MS table and subtable column access.
108 ROMSMainColumns msCols(cPKSMS);
109 ROMSDataDescColumns dataDescCols(cPKSMS.dataDescription());
110 ROMSFeedColumns feedCols(cPKSMS.feed());
111 ROMSFieldColumns fieldCols(cPKSMS.field());
112 ROMSPointingColumns pointingCols(cPKSMS.pointing());
113 ROMSPolarizationColumns polarizationCols(cPKSMS.polarization());
114 ROMSSourceColumns sourceCols(cPKSMS.source());
115 ROMSSpWindowColumns spWinCols(cPKSMS.spectralWindow());
116 ROMSStateColumns stateCols(cPKSMS.state());
117 ROMSSysCalColumns sysCalCols(cPKSMS.sysCal());
118 ROMSWeatherColumns weatherCols(cPKSMS.weather());
119 ROMSAntennaColumns antennaCols(cPKSMS.antenna());
120
121 // Column accessors for required columns.
122 cScanNoCol.reference(msCols.scanNumber());
123 cTimeCol.reference(msCols.time());
124 cIntervalCol.reference(msCols.interval());
125
126 cFieldIdCol.reference(msCols.fieldId());
127 cFieldNameCol.reference(fieldCols.name());
128 cFieldDelayDirCol.reference(fieldCols.delayDir());
129
130 cSrcIdCol.reference(fieldCols.sourceId());
131 cSrcId2Col.reference(sourceCols.sourceId());
132 cSrcNameCol.reference(sourceCols.name());
133 cSrcDirCol.reference(sourceCols.direction());
134 cSrcPMCol.reference(sourceCols.properMotion());
135 cSrcRestFrqCol.reference(sourceCols.restFrequency());
136
137 cStateIdCol.reference(msCols.stateId());
138 cObsModeCol.reference(stateCols.obsMode());
139 cCalCol.reference(stateCols.cal());
140 cSigStateCol.reference(stateCols.sig());
141 cRefStateCol.reference(stateCols.ref());
142 cDataDescIdCol.reference(msCols.dataDescId());
143 cSpWinIdCol.reference(dataDescCols.spectralWindowId());
144 cChanFreqCol.reference(spWinCols.chanFreq());
145
146 cWeatherTimeCol.reference(weatherCols.time());
147 cTemperatureCol.reference(weatherCols.temperature());
148 cPressureCol.reference(weatherCols.pressure());
149 cHumidityCol.reference(weatherCols.relHumidity());
150
151 cBeamNoCol.reference(msCols.feed1());
152 cPointingCol.reference(pointingCols.direction());
153 cSigmaCol.reference(msCols.sigma());
154 cNumReceptorCol.reference(feedCols.numReceptors());
155
156 // Optional columns.
157 cHaveTsys = False;
158 cHaveTcal = False;
159 if ((cHaveSrcVel = cPKSMS.source().tableDesc().isColumn("SYSVEL"))) {
160 cSrcVelCol.attach(cPKSMS.source(), "SYSVEL");
161 }
162
163 if (cHaveSysCal && (cHaveTsys = cPKSMS.sysCal().tableDesc().isColumn("TSYS"))) {
164 cTsysCol.attach(cPKSMS.sysCal(), "TSYS");
165 }
166
167 if (cHaveSysCal && (cHaveTcal = cPKSMS.sysCal().tableDesc().isColumn("TCAL"))) {
168 cTcalCol.attach(cPKSMS.sysCal(), "TCAL");
169 }
170
171 if ((cHaveCalFctr = cPKSMS.tableDesc().isColumn("CALFCTR"))) {
172 cCalFctrCol.attach(cPKSMS, "CALFCTR");
173 }
174
175 if ((cHaveBaseLin = cPKSMS.tableDesc().isColumn("BASELIN"))) {
176 cBaseLinCol.attach(cPKSMS, "BASELIN");
177 cBaseSubCol.attach(cPKSMS, "BASESUB");
178 }
179
180 // Spectral data should always be present.
181 haveSpectra = True;
182 cHaveDataCol = False;
183 cHaveCorrectedDataCol = False;
184 //String telName = antennaCols.name()(0);
185 ROMSObservationColumns observationCols(cPKSMS.observation());
186 String telName = observationCols.telescopeName()(0);
187 //cATF = (telName.contains("DA41") || telName.contains("DV01"));
188 cATF = telName.contains("ATF");
189
190 if (cHaveDataCol = cPKSMS.isColumn(MSMainEnums::DATA)) {
191 if (cATF) {
192 //try to read a single baseline interferometeric data
193 //and treat it as single dish data
194 //maybe extended for ALMA commissioning later
195 cDataCol.reference(msCols.data());
196 if (cHaveCorrectedDataCol = cPKSMS.isColumn(MSMainEnums::CORRECTED_DATA)) {
197 //cerr<<"Do have CORRECTED_DATA column"<<endl;
198 cCorrectedDataCol.reference(msCols.correctedData());
199 }
200 }
201 }
202 else {
203 cFloatDataCol.reference(msCols.floatData());
204 }
205 cFlagCol.reference(msCols.flag());
206
207
208 if (cGetXPol = (cPKSMS.isColumn(MSMainEnums::DATA) && (!cATF))) {
209 if ((cHaveXCalFctr = cPKSMS.tableDesc().isColumn("XCALFCTR"))) {
210 cXCalFctrCol.attach(cPKSMS, "XCALFCTR");
211 }
212
213 cDataCol.reference(msCols.data());
214 }
215
216 // Find which beams are present in the data.
217 Vector<Int> beamNos = cBeamNoCol.getColumn();
218 Int maxBeamNo = max(beamNos) + 1;
219 beams.resize(maxBeamNo);
220
221 beams = False;
222 for (uInt irow = 0; irow < beamNos.nelements(); irow++) {
223 beams(beamNos(irow)) = True;
224 }
225
226
227 // Number of IFs.
228 //uInt nIF = dataDescCols.nrow();
229 uInt nIF =spWinCols.nrow();
230 IFs.resize(nIF);
231 IFs = True;
232
233 // Number of polarizations and channels in each IF.
234 ROScalarColumn<Int> spWinIdCol(dataDescCols.spectralWindowId());
235 ROScalarColumn<Int> numChanCol(spWinCols.numChan());
236
237 ROScalarColumn<Int> polIdCol(dataDescCols.polarizationId());
238 ROScalarColumn<Int> numPolCol(polarizationCols.numCorr());
239
240 nChan.resize(nIF);
241 nPol.resize(nIF);
242 for (uInt iIF = 0; iIF < nIF; iIF++) {
243 nChan(iIF) = numChanCol(spWinIdCol(iIF));
244 nPol(iIF) = numPolCol(polIdCol(iIF));
245 }
246
247 // Cross-polarization data present?
248 haveXPol.resize(nIF);
249 haveXPol = False;
250
251 if (cGetXPol && !(cATF)) {
252 for (Int irow = 0; irow < cNRow; irow++) {
253 if (cDataCol.isDefined(irow)) {
254 Int iIF = cDataDescIdCol(irow);
255 haveXPol(iIF) = True;
256 }
257 }
258 }
259
260
261 // Initialize member data.
262 cBeams.assign(beams);
263 cIFs.assign(IFs);
264 cNChan.assign(nChan);
265 cNPol.assign(nPol);
266 cHaveXPol.assign(haveXPol);
267
268
269 // Default channel range selection.
270 cStartChan.resize(nIF);
271 cEndChan.resize(nIF);
272 cRefChan.resize(nIF);
273
274 for (uInt iIF = 0; iIF < nIF; iIF++) {
275 cStartChan(iIF) = 1;
276 cEndChan(iIF) = cNChan(iIF);
277 cRefChan(iIF) = cNChan(iIF)/2 + 1;
278 }
279
280 Slice all;
281 cDataSel.resize(nIF);
282 cDataSel = Slicer(all, all);
283
284 cScanNo = 0;
285 cCycleNo = 1;
286 cTime = cTimeCol(0);
287
288 return 0;
289}
290
291//---------------------------------------------------- PKSMS2reader::getHeader
292
293// Get parameters describing the data.
294
295Int PKSMS2reader::getHeader(
296 String &observer,
297 String &project,
298 String &antName,
299 Vector<Double> &antPosition,
300 String &obsMode,
301 Float &equinox,
302 String &dopplerFrame,
303 Double &mjd,
304 Double &refFreq,
305 Double &bandwidth,
306 String &fluxunit)
307{
308 if (!cMSopen) {
309 return 1;
310 }
311
312 // Observer and project.
313 ROMSObservationColumns observationCols(cPKSMS.observation());
314 observer = observationCols.observer()(0);
315 project = observationCols.project()(0);
316
317 // Antenna name and ITRF coordinates.
318 ROMSAntennaColumns antennaCols(cPKSMS.antenna());
319 antName = antennaCols.name()(0);
320 antPosition = antennaCols.position()(0);
321
322 // Observation type.
323 if (cObsModeCol.nrow()) {
324 obsMode = cObsModeCol(0);
325 if (obsMode == "\0") obsMode = "RF";
326 } else {
327 obsMode = "RF";
328 }
329
330 fluxunit = "";
331 if (cHaveDataCol) {
332 const TableRecord& keywordSet2
333 = cDataCol.columnDesc().keywordSet();
334 if(keywordSet2.isDefined("UNIT")) {
335 fluxunit = keywordSet2.asString("UNIT");
336 }
337 } else {
338 const TableRecord& keywordSet
339 = cFloatDataCol.columnDesc().keywordSet();
340 if(keywordSet.isDefined("UNIT")) {
341 fluxunit = keywordSet.asString("UNIT");
342 }
343 }
344
345/***
346 const TableRecord& keywordSet
347 = cFloatDataCol.columnDesc().keywordSet();
348 if(keywordSet.isDefined("UNIT")) {
349 fluxunit = keywordSet.asString("UNIT");
350 }
351***/
352 // Coordinate equinox.
353 ROMSPointingColumns pointingCols(cPKSMS.pointing());
354 String dirref = pointingCols.direction().keywordSet().asRecord("MEASINFO").
355 asString("Ref");
356 sscanf(dirref.chars()+1, "%f", &equinox);
357
358 // Frequency/velocity reference frame.
359 ROMSSpWindowColumns spWinCols(cPKSMS.spectralWindow());
360 dopplerFrame = MFrequency::showType(spWinCols.measFreqRef()(0));
361
362 // Translate to FITS standard names.
363 if (dopplerFrame == "TOPO") {
364 dopplerFrame = "TOPOCENT";
365 } else if (dopplerFrame == "GEO") {
366 dopplerFrame = "GEOCENTR";
367 } else if (dopplerFrame == "BARY") {
368 dopplerFrame = "BARYCENT";
369 } else if (dopplerFrame == "GALACTO") {
370 dopplerFrame = "GALACTOC";
371 } else if (dopplerFrame == "LGROUP") {
372 dopplerFrame = "LOCALGRP";
373 } else if (dopplerFrame == "CMB") {
374 dopplerFrame = "CMBDIPOL";
375 } else if (dopplerFrame == "REST") {
376 dopplerFrame = "SOURCE";
377 }
378
379 // MJD at start of observation.
380 mjd = cTimeCol(0)/86400.0;
381
382 // Reference frequency and bandwidth.
383 refFreq = spWinCols.refFrequency()(0);
384 bandwidth = spWinCols.totalBandwidth()(0);
385
386 return 0;
387}
388
389//-------------------------------------------------- PKSMS2reader::getFreqInfo
390
391// Get frequency parameters for each IF.
392
393Int PKSMS2reader::getFreqInfo(
394 Vector<Double> &startFreq,
395 Vector<Double> &endFreq)
396{
397 uInt nIF = cIFs.nelements();
398 startFreq.resize(nIF);
399 endFreq.resize(nIF);
400
401 for (uInt iIF = 0; iIF < nIF; iIF++) {
402 Vector<Double> chanFreq = cChanFreqCol(iIF);
403
404 Int nChan = chanFreq.nelements();
405 startFreq(iIF) = chanFreq(0);
406 endFreq(iIF) = chanFreq(nChan-1);
407 }
408
409 return 0;
410}
411
412//------------------------------------------------------- PKSMS2reader::select
413
414// Set data selection by beam number and channel.
415
416uInt PKSMS2reader::select(
417 const Vector<Bool> beamSel,
418 const Vector<Bool> IFsel,
419 const Vector<Int> startChan,
420 const Vector<Int> endChan,
421 const Vector<Int> refChan,
422 const Bool getSpectra,
423 const Bool getXPol,
424 const Bool getFeedPos)
425{
426 if (!cMSopen) {
427 return 1;
428 }
429
430 // Beam selection.
431 uInt nBeam = cBeams.nelements();
432 uInt nBeamSel = beamSel.nelements();
433 for (uInt ibeam = 0; ibeam < nBeam; ibeam++) {
434 if (ibeam < nBeamSel) {
435 cBeams(ibeam) = beamSel(ibeam);
436 } else {
437 cBeams(ibeam) = False;
438 }
439 }
440
441 uInt nIF = cIFs.nelements();
442 uInt maxNChan = 0;
443 for (uInt iIF = 0; iIF < nIF; iIF++) {
444 // IF selection.
445 if (iIF < IFsel.nelements()) {
446 cIFs(iIF) = IFsel(iIF);
447 } else {
448 cIFs(iIF) = False;
449 }
450
451 if (!cIFs(iIF)) continue;
452
453
454 // Channel selection.
455 if (iIF < startChan.nelements()) {
456 cStartChan(iIF) = startChan(iIF);
457
458 if (cStartChan(iIF) <= 0) {
459 cStartChan(iIF) += cNChan(iIF);
460 } else if (cStartChan(iIF) > Int(cNChan(iIF))) {
461 cStartChan(iIF) = cNChan(iIF);
462 }
463 }
464
465 if (iIF < endChan.nelements()) {
466 cEndChan(iIF) = endChan(iIF);
467
468 if (cEndChan(iIF) <= 0) {
469 cEndChan(iIF) += cNChan(iIF);
470 } else if (cEndChan(iIF) > Int(cNChan(iIF))) {
471 cEndChan(iIF) = cNChan(iIF);
472 }
473 }
474
475 if (iIF < refChan.nelements()) {
476 cRefChan(iIF) = refChan(iIF);
477 } else {
478 cRefChan(iIF) = cStartChan(iIF);
479 if (cStartChan(iIF) <= cEndChan(iIF)) {
480 cRefChan(iIF) += (cEndChan(iIF) - cStartChan(iIF) + 1)/2;
481 } else {
482 cRefChan(iIF) -= (cStartChan(iIF) - cEndChan(iIF) + 1)/2;
483 }
484 }
485
486 uInt nChan = abs(cEndChan(iIF) - cStartChan(iIF)) + 1;
487 if (maxNChan < nChan) {
488 maxNChan = nChan;
489 }
490
491 // Inverted Slices are not allowed.
492 Slice outPols;
493 Slice outChans(min(cStartChan(iIF),cEndChan(iIF))-1, nChan);
494 cDataSel(iIF) = Slicer(outPols, outChans);
495 }
496
497 // Get spectral data?
498 cGetSpectra = getSpectra;
499
500 // Get cross-polarization data?
501 cGetXPol = cGetXPol && getXPol;
502
503 // Get feed positions? (Not available.)
504 cGetFeedPos = False;
505
506 return maxNChan;
507}
508
509//---------------------------------------------------- PKSMS2reader::findRange
510
511// Find the range of the data in time and position.
512
513Int PKSMS2reader::findRange(
514 Int &nRow,
515 Int &nSel,
516 Vector<Double> &timeSpan,
517 Matrix<Double> &positions)
518{
519 if (!cMSopen) {
520 return 1;
521 }
522
523 nRow = cNRow;
524
525 // Find the number of rows selected.
526 nSel = 0;
527 Vector<Bool> sel(nRow);
528 for (Int irow = 0; irow < nRow; irow++) {
529 if ((sel(irow) = cBeams(cBeamNoCol(irow)) &&
530 cIFs(cDataDescIdCol(irow)))) {
531 nSel++;
532 }
533 }
534
535 // Find the time range (s).
536 timeSpan.resize(2);
537 timeSpan(0) = cTimeCol(0);
538 timeSpan(1) = cTimeCol(nRow-1);
539
540 // Retrieve positions for selected data.
541 Int isel = 0;
542 positions.resize(2,nSel);
543 for (Int irow = 0; irow < nRow; irow++) {
544 if (sel(irow)) {
545 Matrix<Double> pointingDir = cPointingCol(cFieldIdCol(irow));
546 positions.column(isel++) = pointingDir.column(0);
547 }
548 }
549
550 return 0;
551}
552
553//--------------------------------------------------------- PKSMS2reader::read
554
555// Read the next data record.
556
557Int PKSMS2reader::read(
558 Int &scanNo,
559 Int &cycleNo,
560 Double &mjd,
561 Double &interval,
562 String &fieldName,
563 String &srcName,
564 Vector<Double> &srcDir,
565 Vector<Double> &srcPM,
566 Double &srcVel,
567 String &obsMode,
568 Int &IFno,
569 Double &refFreq,
570 Double &bandwidth,
571 Double &freqInc,
572 Vector<Double> &restFreq,
573 Vector<Float> &tcal,
574 String &tcalTime,
575 Float &azimuth,
576 Float &elevation,
577 Float &parAngle,
578 Float &focusAxi,
579 Float &focusTan,
580 Float &focusRot,
581 Float &temperature,
582 Float &pressure,
583 Float &humidity,
584 Float &windSpeed,
585 Float &windAz,
586 Int &refBeam,
587 Int &beamNo,
588 Vector<Double> &direction,
589 Vector<Double> &scanRate,
590 Vector<Float> &tsys,
591 Vector<Float> &sigma,
592 Vector<Float> &calFctr,
593 Matrix<Float> &baseLin,
594 Matrix<Float> &baseSub,
595 Matrix<Float> &spectra,
596 Matrix<uChar> &flagged,
597 Complex &xCalFctr,
598 Vector<Complex> &xPol)
599{
600 if (!cMSopen) {
601 return 1;
602 }
603
604 // Check for EOF.
605 if (cIdx >= cNRow) {
606 return -1;
607 }
608
609 // Find the next selected beam and IF.
610 Int ibeam;
611 Int iIF;
612 Int iDataDesc;
613
614 while (True) {
615 ibeam = cBeamNoCol(cIdx);
616 iDataDesc = cDataDescIdCol(cIdx);
617 iIF =cSpWinIdCol(iDataDesc);
618 if (cBeams(ibeam) && cIFs(iIF)) {
619 break;
620 }
621
622 // Check for EOF.
623 if (++cIdx >= cNRow) {
624 return -1;
625 }
626 }
627 // Renumerate scan no. Here still is 1-based
628 //scanNo = cScanNoCol(cIdx) - cScanNoCol(0) + 1;
629 scanNo = cScanNoCol(cIdx);
630
631 if (scanNo != cScanNo) {
632 // Start of new scan.
633 cScanNo = scanNo;
634 cCycleNo = 1;
635 cTime = cTimeCol(cIdx);
636 }
637
638 Double time = cTimeCol(cIdx);
639 mjd = time/86400.0;
640 interval = cIntervalCol(cIdx);
641
642 // Reconstruct the integration cycle number; due to small latencies the
643 // integration time is usually slightly less than the time between cycles,
644 // resetting cTime will prevent the difference from accumulating.
645 cCycleNo += nint((time - cTime)/interval);
646 cycleNo = cCycleNo;
647 cTime = time;
648
649 Int fieldId = cFieldIdCol(cIdx);
650 fieldName = cFieldNameCol(fieldId);
651
652 Int srcId = cSrcIdCol(fieldId);
653 //For source with multiple spectral window setting, this is not
654 // correct. Source name of srcId may not be at 'srcId'th row of SrcNameCol
655 //srcName = cSrcNameCol(srcId);
656 for (uInt irow = 0; irow < cSrcId2Col.nrow(); irow++) {
657 if (cSrcId2Col(irow) == srcId) {
658 srcName = cSrcNameCol(irow);
659 }
660 }
661
662 srcDir = cSrcDirCol(srcId);
663 srcPM = cSrcPMCol(srcId);
664
665 // Systemic velocity.
666 if (!cHaveSrcVel || cATF) {
667 srcVel = 0.0f;
668 } else {
669 srcVel = cSrcVelCol(srcId)(IPosition(1,0));
670 }
671
672 ROMSAntennaColumns antennaCols(cPKSMS.antenna());
673 String telescope = antennaCols.name()(0);
674 Bool cGBT = telescope.contains("GBT");
675 // Observation type.
676 // check if State Table exist
677 //Bool cHaveStateTab=Table::isReadable(cPKSMS.stateTableName());
678 Int stateId = 0;
679 Int StateNRow = 0;
680 StateNRow=cObsModeCol.nrow();
681 if (Table::isReadable(cPKSMS.stateTableName())) {
682 obsMode = " ";
683 if (StateNRow > 0) {
684 stateId = cStateIdCol(cIdx);
685 if (stateId == -1) {
686 //obsMode = " ";
687 } else {
688 obsMode = cObsModeCol(stateId);
689 Bool sigState =cSigStateCol(stateId);
690 Bool refState =cRefStateCol(stateId);
691 //DEBUG
692 //cerr <<"stateid="<<stateId<<" obsmode="<<obsMode<<endl;
693 if (cGBT) {
694 // split the obsMode string and append a proper label
695 // (these are GBT specific)
696 int epos = obsMode.find_first_of(':');
697 int nextpos = obsMode.find_first_of(':',epos+1);
698 string obsMode1 = obsMode.substr(0,epos);
699 string obsMode2 = obsMode.substr(epos+1,nextpos-epos-1);
700
701 //cerr <<"obsMode2= "<<obsMode2<<endl;
702 if (!srcName.contains("_ps")
703 &&!srcName.contains("_psr")
704 &&!srcName.contains("_nod")
705 &&!srcName.contains("_fs")
706 &&!srcName.contains("_fsr")) {
707 // if Nod mode observation , append '_nod'
708 if (obsMode1 == "Nod") {
709 srcName.append("_nod");
710 } else if (obsMode1 == "OffOn") {
711 // for GBT position switch observations (OffOn or OnOff)
712 if (obsMode2 == "PSWITCHON") srcName.append("_ps");
713 if (obsMode2 == "PSWITCHOFF") srcName.append("_psr");
714 } else {
715 if (obsMode2 == "FSWITCH") {
716 // for GBT frequency switch mode
717 if (sigState) srcName.append("_fs");
718 if (refState) srcName.append("_fsr");
719 }
720 }
721 }
722 }
723 }
724 }
725 }
726 // CAL state
727 // this should be apply just for GBT data?
728 Double Cal;
729 if (stateId==-1 || StateNRow==0) {
730 Cal = 0;
731 } else {
732 Cal = cCalCol(stateId);
733 }
734 if (cGBT) {
735 if (Cal > 0 && !srcName.contains("_calon")) {
736 srcName.append("_calon");
737 }
738 }
739
740 IFno = iIF + 1;
741 Int nChan = abs(cEndChan(iIF) - cStartChan(iIF)) + 1;
742
743 // Minimal handling on continuum data.
744 Vector<Double> chanFreq = cChanFreqCol(iIF);
745 if (nChan == 1) {
746 freqInc = chanFreq(0);
747 refFreq = chanFreq(0);
748 restFreq = 0.0f;
749 } else {
750
751 if (cStartChan(iIF) <= cEndChan(iIF)) {
752 freqInc = chanFreq(1) - chanFreq(0);
753 } else {
754 freqInc = chanFreq(0) - chanFreq(1);
755 }
756 refFreq = chanFreq(cRefChan(iIF)-1);
757 Bool HaveSrcRestFreq= cSrcRestFrqCol.isDefined(srcId);
758 if (HaveSrcRestFreq) {
759 //restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0));
760 restFreq = cSrcRestFrqCol(srcId);
761 } else {
762 restFreq = 0.0f;
763 }
764 }
765 bandwidth = abs(freqInc * nChan);
766
767 tcal.resize(cNPol(iIF));
768 tcal = 0.0f;
769 tcalTime = "";
770 //azimuth = 0.0f;
771 //elevation = 0.0f;
772 parAngle = 0.0f;
773 focusAxi = 0.0f;
774 focusTan = 0.0f;
775 focusRot = 0.0f;
776
777 // Find the appropriate entry in the WEATHER subtable.
778 //Bool cHaveStateTab=Table::isReadable(cPKSMS.stateTableName());
779 Bool cHaveWeatherTab = Table::isReadable(cPKSMS.weatherTableName());
780 Int weatherIdx=-1;
781 if (cHaveWeatherTab) {
782 Vector<Double> wTimes = cWeatherTimeCol.getColumn();
783 for (weatherIdx = wTimes.nelements()-1; weatherIdx >= 0; weatherIdx--) {
784 if (cWeatherTimeCol(weatherIdx) <= time) {
785 break;
786 }
787 }
788 }
789 if (weatherIdx < 0 || !cHaveWeatherTab) {
790 // No appropriate WEATHER entry.
791 pressure = 0.0f;
792 humidity = 0.0f;
793 temperature = 0.0f;
794 } else {
795 pressure = cPressureCol(weatherIdx);
796 humidity = cHumidityCol(weatherIdx);
797 temperature = cTemperatureCol(weatherIdx);
798 }
799
800 windSpeed = 0.0f;
801 windAz = 0.0f;
802
803 refBeam = 0;
804 beamNo = ibeam + 1;
805
806 //Matrix<Double> pointingDir = cPointingCol(fieldId);
807 //pointingDir = cPointingCol(fieldId);
808 //direction = pointingDir.column(0);
809 //uInt ncols = pointingDir.ncolumn();
810 //if (ncols == 1) {
811 // scanRate = 0.0f;
812 //} else {
813 // scanRate = pointingDir.column(1);
814 //}
815
816 // Get direction from FIELD table
817 // here, assume direction to be the field direction not pointing
818 Matrix<Double> delayDir = cFieldDelayDirCol(fieldId);
819 direction = delayDir.column(0);
820 uInt ncols = delayDir.ncolumn();
821 if (ncols == 1) {
822 scanRate = 0.0f;
823 } else {
824 scanRate = delayDir.column(1);
825 }
826
827 // caluculate azimuth and elevation
828 // first, get the reference frame
829 MVPosition mvpos(antennaCols.position()(0));
830 MPosition mp(mvpos);
831 Quantum<Double> qt(time,"s");
832 MVEpoch mvt(qt);
833 MEpoch me(mvt);
834 MeasFrame frame(mp, me);
835 //
836 ROMSFieldColumns fldCols(cPKSMS.field());
837 Vector<MDirection> vmd(1);
838 MDirection md;
839 fldCols.delayDirMeasCol().get(fieldId,vmd);
840 md = vmd[0];
841 //Vector<Double> dircheck = md.getAngle("rad").getValue();
842 //cerr<<"dircheck="<<dircheck<<endl;
843
844 Vector<Double> azel =
845 MDirection::Convert(md, MDirection::Ref(MDirection::AZEL,
846 frame)
847 )().getAngle("rad").getValue();
848 //cerr<<"azel="<<azel<<endl;
849 azimuth = azel[0];
850 elevation = azel[1];
851
852 // Get Tsys assuming that entries in the SYSCAL table match the main table.
853 if (cHaveTsys) {
854 Int nTsysColRow = cTsysCol.nrow();
855 if (nTsysColRow != cNRow) {
856 cHaveTsys=0;
857 }
858 }
859 if (cHaveTsys) {
860 cTsysCol.get(cIdx, tsys, True);
861 } else {
862 Int numReceptor;
863 cNumReceptorCol.get(0, numReceptor);
864 tsys.resize(numReceptor);
865 tsys = 1.0f;
866 }
867 cSigmaCol.get(cIdx, sigma, True);
868
869 //get Tcal if available
870 if (cHaveTcal) {
871 Int nTcalColRow = cTcalCol.nrow();
872 uInt nBeam = cBeams.nelements();
873 uInt nIF = cIFs.nelements();
874 uInt nrws = nBeam * nIF;
875 if (nTcalColRow > 0) {
876 // find tcal match with the data with the data time stamp
877 Double mjds = mjd*(24*3600);
878 Double dtcalTime;
879 if ( mjd > lastmjd || cIdx==0 ) {
880 //Table tmptab = cSysCalTab(near(cSysCalTab.col("TIME"),mjds));
881 tmptab = cSysCalTab(near(cSysCalTab.col("TIME"),mjds), nrws);
882 //DEBUG
883 //if (cIdx == 0) {
884 // cerr<<"inital table retrieved"<<endl;
885 //}
886
887 }
888
889 if (nBeam == 1) {
890 tmptab2 = tmptab( tmptab.col("SPECTRAL_WINDOW_ID") == iIF, 1);
891 } else {
892 tmptab2 = tmptab( tmptab.col("SPECTRAL_WINDOW_ID") == iIF &&
893 tmptab.col("FEED_ID") == ibeam , 1);
894 }
895 //cerr<<"first subtab rows="<<tmptab.nrow()<<endl;
896 int syscalrow = tmptab2.nrow();
897 ROArrayColumn<Float> tcalCol(tmptab2, "TCAL");
898 ROScalarColumn<Double> tcalTimeCol(tmptab2, "TIME");
899 if (syscalrow==0) {
900 cerr<<"Cannot find any matching Tcal at/near the data timestamp."
901 << " Set Tcal=0.0"<<endl;
902 } else {
903 tcalCol.get(0, tcal);
904 tcalTimeCol.get(0,dtcalTime);
905 tcalTime = MVTime(dtcalTime/(24*3600)).string(MVTime::YMD);
906 //DEBUG
907 //cerr<<"cIdx:"<<cIdx<<" tcal="<<tcal<<" tcalTime="<<tcalTime<<endl;
908 tmptab.markForDelete();
909 tmptab2.markForDelete();
910 }
911 }
912 lastmjd = mjd;
913 }
914
915 // Calibration factors (if available).
916 calFctr.resize(cNPol(iIF));
917 if (cHaveCalFctr) {
918 cCalFctrCol.get(cIdx, calFctr);
919 } else {
920 calFctr = 0.0f;
921 }
922
923 // Baseline parameters (if available).
924 if (cHaveBaseLin) {
925 baseLin.resize(2,cNPol(iIF));
926 cBaseLinCol.get(cIdx, baseLin);
927
928 baseSub.resize(9,cNPol(iIF));
929 cBaseSubCol.get(cIdx, baseSub);
930
931 } else {
932 baseLin.resize(0,0);
933 baseSub.resize(0,0);
934 }
935 // Get spectral data.
936 if (cGetSpectra) {
937 Matrix<Float> tmpData;
938 Matrix<Bool> tmpFlag;
939 if (cHaveDataCol) {
940 Matrix<Complex> tmpCmplxData;
941 Matrix<Float> tmpReData;
942 Matrix<Float> tmpImData;
943 //cerr<<"reading spectra..."<<endl;
944 //# TODO - should have a flag to user to select DATA or CORRECTED_DATA
945 //# currently just automatically determined, --- read CORRECTED one
946 //# if the column exist.
947 if (cHaveCorrectedDataCol) {
948 cCorrectedDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True);
949 } else {
950 cDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True);
951 }
952 tmpReData = real(tmpCmplxData);
953 tmpImData = imag(tmpCmplxData);
954 tmpData = sqrt(tmpReData*tmpReData + tmpImData*tmpImData);
955 } else {
956 cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
957 }
958 cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
959
960 // Transpose spectra.
961 Int nPol = tmpData.nrow();
962 spectra.resize(nChan, nPol);
963 flagged.resize(nChan, nPol);
964 if (cEndChan(iIF) >= cStartChan(iIF)) {
965 // Simple transposition.
966 for (Int ipol = 0; ipol < nPol; ipol++) {
967 for (Int ichan = 0; ichan < nChan; ichan++) {
968 spectra(ichan,ipol) = tmpData(ipol,ichan);
969 flagged(ichan,ipol) = tmpFlag(ipol,ichan);
970 }
971 }
972
973 } else {
974 // Transpose with inversion.
975 Int jchan = nChan - 1;
976 for (Int ipol = 0; ipol < nPol; ipol++) {
977 for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
978 spectra(ichan,ipol) = tmpData(ipol,jchan);
979 flagged(ichan,ipol) = tmpFlag(ipol,jchan);
980 }
981 }
982 }
983 }
984
985 // Get cross-polarization data.
986 if (cGetXPol) {
987 //cerr<<"cGetXPol="<<cGetXPol<<endl;
988 //cerr<<"cHaveXCalFctr="<<cHaveXCalFctr<<endl;
989
990 if (cHaveXCalFctr) {
991 cXCalFctrCol.get(cIdx, xCalFctr);
992 } else {
993 xCalFctr = Complex(0.0f, 0.0f);
994 }
995
996 if(!cATF) {
997 cDataCol.get(cIdx, xPol, True);
998
999 if (cEndChan(iIF) < cStartChan(iIF)) {
1000 Complex ctmp;
1001 Int jchan = nChan - 1;
1002 for (Int ichan = 0; ichan < nChan/2; ichan++, jchan--) {
1003 ctmp = xPol(ichan);
1004 xPol(ichan) = xPol(jchan);
1005 xPol(jchan) = ctmp;
1006 }
1007 }
1008 }
1009 }
1010 /**
1011 cerr<<"scanNo="<<scanNo<<endl;
1012 cerr<<"cycleNo="<<cycleNo<<endl;
1013 cerr<<"mjd="<<mjd<<endl;
1014 cerr<<"interval="<<interval<<endl;
1015 cerr<<"fieldName="<<fieldName<<endl;
1016 cerr<<"srcNmae="<<srcName<<endl;
1017 cerr<<"srcDir="<<srcDir<<endl;
1018 cerr<<"srcPM="<<srcPM<<endl;
1019 cerr<<"srcVel="<<srcVel<<endl;
1020 cerr<<"obsMode="<<obsMode<<endl;
1021 cerr<<"IFno="<<IFno<<endl;
1022 cerr<<"refFreq="<<refFreq<<endl;
1023 cerr<<"tcal="<<tcal<<endl;
1024 cerr<<"direction="<<direction<<endl;
1025 cerr<<"scanRate="<<scanRate<<endl;
1026 cerr<<"tsys="<<tsys<<endl;
1027 cerr<<"sigma="<<sigma<<endl;
1028 cerr<<"calFctr="<<calFctr<<endl;
1029 cerr<<"baseLin="<<baseLin<<endl;
1030 cerr<<"baseSub="<<baseSub<<endl;
1031 cerr<<"spectra="<<spectra.shape()<<endl;
1032 cerr<<"flagged="<<flagged.shape()<<endl;
1033 cerr<<"xCalFctr="<<xCalFctr<<endl;
1034 cerr<<"xPol="<<xPol<<endl;
1035 **/
1036 cIdx++;
1037
1038 return 0;
1039}
1040
1041//--------------------------------------------------------- PKSMS2reader::read
1042
1043// Read the next data record, just the basics.
1044
1045Int PKSMS2reader::read(
1046 Int &IFno,
1047 Vector<Float> &tsys,
1048 Vector<Float> &calFctr,
1049 Matrix<Float> &baseLin,
1050 Matrix<Float> &baseSub,
1051 Matrix<Float> &spectra,
1052 Matrix<uChar> &flagged)
1053{
1054 if (!cMSopen) {
1055 return 1;
1056 }
1057
1058 // Check for EOF.
1059 if (cIdx >= cNRow) {
1060 return -1;
1061 }
1062
1063 // Find the next selected beam and IF.
1064 Int ibeam;
1065 Int iIF;
1066 Int iDataDesc;
1067 while (True) {
1068 ibeam = cBeamNoCol(cIdx);
1069 //iIF = cDataDescIdCol(cIdx);
1070 iDataDesc = cDataDescIdCol(cIdx);
1071 iIF = cSpWinIdCol(iDataDesc);
1072 if (cBeams(ibeam) && cIFs(iIF)) {
1073 break;
1074 }
1075
1076 // Check for EOF.
1077 if (++cIdx >= cNRow) {
1078 return -1;
1079 }
1080 }
1081
1082 IFno = iIF + 1;
1083 // Get Tsys assuming that entries in the SYSCAL table match the main table.
1084 cTsysCol.get(cIdx, tsys, True);
1085
1086 // Calibration factors (if available).
1087 if (cHaveCalFctr) {
1088 cCalFctrCol.get(cIdx, calFctr, True);
1089 } else {
1090 calFctr.resize(cNPol(iIF));
1091 calFctr = 0.0f;
1092 }
1093
1094 // Baseline parameters (if available).
1095 if (cHaveBaseLin) {
1096 baseLin.resize(2,cNPol(iIF));
1097 cBaseLinCol.get(cIdx, baseLin);
1098
1099 baseSub.resize(9,cNPol(iIF));
1100 cBaseSubCol.get(cIdx, baseSub);
1101
1102 } else {
1103 baseLin.resize(0,0);
1104 baseSub.resize(0,0);
1105 }
1106
1107 if (cGetSpectra) {
1108 // Get spectral data.
1109 Matrix<Float> tmpData;
1110 Matrix<Bool> tmpFlag;
1111 if (cHaveDataCol) {
1112 Matrix<Complex> tmpCmplxData;
1113 cDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True);
1114 tmpData = real(tmpCmplxData);
1115 } else {
1116 cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
1117 }
1118 cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
1119
1120 // Transpose spectra.
1121 Int nChan = tmpData.ncolumn();
1122 Int nPol = tmpData.nrow();
1123 spectra.resize(nChan, nPol);
1124 flagged.resize(nChan, nPol);
1125 if (cEndChan(iIF) >= cStartChan(iIF)) {
1126 // Simple transposition.
1127 for (Int ipol = 0; ipol < nPol; ipol++) {
1128 for (Int ichan = 0; ichan < nChan; ichan++) {
1129 spectra(ichan,ipol) = tmpData(ipol,ichan);
1130 flagged(ichan,ipol) = tmpFlag(ipol,ichan);
1131 }
1132 }
1133
1134 } else {
1135 // Transpose with inversion.
1136 Int jchan = nChan - 1;
1137 for (Int ipol = 0; ipol < nPol; ipol++) {
1138 for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
1139 spectra(ichan,ipol) = tmpData(ipol,jchan);
1140 flagged(ichan,ipol) = tmpFlag(ipol,jchan);
1141 }
1142 }
1143 }
1144 }
1145
1146 cIdx++;
1147
1148 return 0;
1149}
1150
1151//-------------------------------------------------------- PKSMS2reader::close
1152
1153// Close the MS.
1154
1155void PKSMS2reader::close()
1156{
1157 cPKSMS = MeasurementSet();
1158 cMSopen = False;
1159}
Note: See TracBrowser for help on using the repository browser.