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

Last change on this file since 1431 was 1393, checked in by TakTsutsumi, 17 years ago

Changes to handle GBT MS data
File size: 28.0 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 String tmpcalName = msName;
95 tmpcalName.append("/SYSCAL");
96 cSysCalTab = Table(tmpcalName);
97
98 cIdx = 0;
99 lastmjd = 0.0;
100 cNRow = cPKSMS.nrow();
101 cMSopen = True;
102
103 // Lock the table for read access.
104 cPKSMS.lock(False);
105
106 // Main MS table and subtable column access.
107 ROMSMainColumns msCols(cPKSMS);
108 ROMSDataDescColumns dataDescCols(cPKSMS.dataDescription());
109 ROMSFeedColumns feedCols(cPKSMS.feed());
110 ROMSFieldColumns fieldCols(cPKSMS.field());
111 ROMSPointingColumns pointingCols(cPKSMS.pointing());
112 ROMSPolarizationColumns polarizationCols(cPKSMS.polarization());
113 ROMSSourceColumns sourceCols(cPKSMS.source());
114 ROMSSpWindowColumns spWinCols(cPKSMS.spectralWindow());
115 ROMSStateColumns stateCols(cPKSMS.state());
116 ROMSSysCalColumns sysCalCols(cPKSMS.sysCal());
117 ROMSWeatherColumns weatherCols(cPKSMS.weather());
118
119 // Column accessors for required columns.
120 cScanNoCol.reference(msCols.scanNumber());
121 cTimeCol.reference(msCols.time());
122 cIntervalCol.reference(msCols.interval());
123
124 cFieldIdCol.reference(msCols.fieldId());
125 cFieldNameCol.reference(fieldCols.name());
126 cFieldDelayDirCol.reference(fieldCols.delayDir());
127
128 cSrcIdCol.reference(fieldCols.sourceId());
129 cSrcId2Col.reference(sourceCols.sourceId());
130 cSrcNameCol.reference(sourceCols.name());
131 cSrcDirCol.reference(sourceCols.direction());
132 cSrcPMCol.reference(sourceCols.properMotion());
133 cSrcRestFrqCol.reference(sourceCols.restFrequency());
134
135 cStateIdCol.reference(msCols.stateId());
136 cObsModeCol.reference(stateCols.obsMode());
137 cCalCol.reference(stateCols.cal());
138 cSigStateCol.reference(stateCols.sig());
139 cRefStateCol.reference(stateCols.ref());
140 cDataDescIdCol.reference(msCols.dataDescId());
141 cSpWinIdCol.reference(dataDescCols.spectralWindowId());
142 cChanFreqCol.reference(spWinCols.chanFreq());
143
144 cWeatherTimeCol.reference(weatherCols.time());
145 cTemperatureCol.reference(weatherCols.temperature());
146 cPressureCol.reference(weatherCols.pressure());
147 cHumidityCol.reference(weatherCols.relHumidity());
148
149 cBeamNoCol.reference(msCols.feed1());
150 cPointingCol.reference(pointingCols.direction());
151 cSigmaCol.reference(msCols.sigma());
152 cNumReceptorCol.reference(feedCols.numReceptors());
153
154 // Optional columns.
155 if ((cHaveSrcVel = cPKSMS.source().tableDesc().isColumn("SYSVEL"))) {
156 cSrcVelCol.attach(cPKSMS.source(), "SYSVEL");
157 }
158
159 if ((cHaveTsys = cPKSMS.sysCal().tableDesc().isColumn("TSYS"))) {
160 cTsysCol.attach(cPKSMS.sysCal(), "TSYS");
161 }
162
163 if ((cHaveTcal = cPKSMS.sysCal().tableDesc().isColumn("TCAL"))) {
164 cTcalCol.attach(cPKSMS.sysCal(), "TCAL");
165 }
166
167 if ((cHaveCalFctr = cPKSMS.tableDesc().isColumn("CALFCTR"))) {
168 cCalFctrCol.attach(cPKSMS, "CALFCTR");
169 }
170
171 if ((cHaveBaseLin = cPKSMS.tableDesc().isColumn("BASELIN"))) {
172 cBaseLinCol.attach(cPKSMS, "BASELIN");
173 cBaseSubCol.attach(cPKSMS, "BASESUB");
174 }
175
176 // Spectral data should always be present.
177 haveSpectra = True;
178 cFloatDataCol.reference(msCols.floatData());
179 cFlagCol.reference(msCols.flag());
180
181
182 if ((cGetXPol = cPKSMS.isColumn(MSMainEnums::DATA))) {
183 if ((cHaveXCalFctr = cPKSMS.tableDesc().isColumn("XCALFCTR"))) {
184 cXCalFctrCol.attach(cPKSMS, "XCALFCTR");
185 }
186
187 cDataCol.reference(msCols.data());
188 }
189
190 // Find which beams are present in the data.
191 Vector<Int> beamNos = cBeamNoCol.getColumn();
192 Int maxBeamNo = max(beamNos) + 1;
193 beams.resize(maxBeamNo);
194
195 beams = False;
196 for (uInt irow = 0; irow < beamNos.nelements(); irow++) {
197 beams(beamNos(irow)) = True;
198 }
199
200
201 // Number of IFs.
202 //uInt nIF = dataDescCols.nrow();
203 uInt nIF =spWinCols.nrow();
204 IFs.resize(nIF);
205 IFs = True;
206
207 // Number of polarizations and channels in each IF.
208 ROScalarColumn<Int> spWinIdCol(dataDescCols.spectralWindowId());
209 ROScalarColumn<Int> numChanCol(spWinCols.numChan());
210
211 ROScalarColumn<Int> polIdCol(dataDescCols.polarizationId());
212 ROScalarColumn<Int> numPolCol(polarizationCols.numCorr());
213
214 nChan.resize(nIF);
215 nPol.resize(nIF);
216 for (uInt iIF = 0; iIF < nIF; iIF++) {
217 nChan(iIF) = numChanCol(spWinIdCol(iIF));
218 nPol(iIF) = numPolCol(polIdCol(iIF));
219 }
220
221 // Cross-polarization data present?
222 haveXPol.resize(nIF);
223 haveXPol = False;
224
225 if (cGetXPol) {
226 for (Int irow = 0; irow < cNRow; irow++) {
227 if (cDataCol.isDefined(irow)) {
228 Int iIF = cDataDescIdCol(irow);
229 haveXPol(iIF) = True;
230 }
231 }
232 }
233
234
235 // Initialize member data.
236 cBeams.assign(beams);
237 cIFs.assign(IFs);
238 cNChan.assign(nChan);
239 cNPol.assign(nPol);
240 cHaveXPol.assign(haveXPol);
241
242
243 // Default channel range selection.
244 cStartChan.resize(nIF);
245 cEndChan.resize(nIF);
246 cRefChan.resize(nIF);
247
248 for (uInt iIF = 0; iIF < nIF; iIF++) {
249 cStartChan(iIF) = 1;
250 cEndChan(iIF) = cNChan(iIF);
251 cRefChan(iIF) = cNChan(iIF)/2 + 1;
252 }
253
254 Slice all;
255 cDataSel.resize(nIF);
256 cDataSel = Slicer(all, all);
257
258 cScanNo = 0;
259 cCycleNo = 1;
260 cTime = cTimeCol(0);
261
262 return 0;
263}
264
265//---------------------------------------------------- PKSMS2reader::getHeader
266
267// Get parameters describing the data.
268
269Int PKSMS2reader::getHeader(
270 String &observer,
271 String &project,
272 String &antName,
273 Vector<Double> &antPosition,
274 String &obsMode,
275 Float &equinox,
276 String &dopplerFrame,
277 Double &mjd,
278 Double &refFreq,
279 Double &bandwidth,
280 String &fluxunit)
281{
282 if (!cMSopen) {
283 return 1;
284 }
285
286 // Observer and project.
287 ROMSObservationColumns observationCols(cPKSMS.observation());
288 observer = observationCols.observer()(0);
289 project = observationCols.project()(0);
290
291 // Antenna name and ITRF coordinates.
292 ROMSAntennaColumns antennaCols(cPKSMS.antenna());
293 antName = antennaCols.name()(0);
294 antPosition = antennaCols.position()(0);
295
296 // Observation type.
297 if (cObsModeCol.nrow()) {
298 obsMode = cObsModeCol(0);
299 if (obsMode == "\0") obsMode = "RF";
300 } else {
301 obsMode = "RF";
302 }
303
304 fluxunit = "";
305 const TableRecord& keywordSet
306 = cFloatDataCol.columnDesc().keywordSet();
307 if(keywordSet.isDefined("UNIT")) {
308 fluxunit = keywordSet.asString("UNIT");
309 }
310
311 // Coordinate equinox.
312 ROMSPointingColumns pointingCols(cPKSMS.pointing());
313 String dirref = pointingCols.direction().keywordSet().asRecord("MEASINFO").
314 asString("Ref");
315 sscanf(dirref.chars()+1, "%f", &equinox);
316
317 // Frequency/velocity reference frame.
318 ROMSSpWindowColumns spWinCols(cPKSMS.spectralWindow());
319 dopplerFrame = MFrequency::showType(spWinCols.measFreqRef()(0));
320
321 // Translate to FITS standard names.
322 if (dopplerFrame == "TOPO") {
323 dopplerFrame = "TOPOCENT";
324 } else if (dopplerFrame == "GEO") {
325 dopplerFrame = "GEOCENTR";
326 } else if (dopplerFrame == "BARY") {
327 dopplerFrame = "BARYCENT";
328 } else if (dopplerFrame == "GALACTO") {
329 dopplerFrame = "GALACTOC";
330 } else if (dopplerFrame == "LGROUP") {
331 dopplerFrame = "LOCALGRP";
332 } else if (dopplerFrame == "CMB") {
333 dopplerFrame = "CMBDIPOL";
334 } else if (dopplerFrame == "REST") {
335 dopplerFrame = "SOURCE";
336 }
337
338 // MJD at start of observation.
339 mjd = cTimeCol(0)/86400.0;
340
341 // Reference frequency and bandwidth.
342 refFreq = spWinCols.refFrequency()(0);
343 bandwidth = spWinCols.totalBandwidth()(0);
344
345 return 0;
346}
347
348//-------------------------------------------------- PKSMS2reader::getFreqInfo
349
350// Get frequency parameters for each IF.
351
352Int PKSMS2reader::getFreqInfo(
353 Vector<Double> &startFreq,
354 Vector<Double> &endFreq)
355{
356 uInt nIF = cIFs.nelements();
357 startFreq.resize(nIF);
358 endFreq.resize(nIF);
359
360 for (uInt iIF = 0; iIF < nIF; iIF++) {
361 Vector<Double> chanFreq = cChanFreqCol(iIF);
362
363 Int nChan = chanFreq.nelements();
364 startFreq(iIF) = chanFreq(0);
365 endFreq(iIF) = chanFreq(nChan-1);
366 }
367
368 return 0;
369}
370
371//------------------------------------------------------- PKSMS2reader::select
372
373// Set data selection by beam number and channel.
374
375uInt PKSMS2reader::select(
376 const Vector<Bool> beamSel,
377 const Vector<Bool> IFsel,
378 const Vector<Int> startChan,
379 const Vector<Int> endChan,
380 const Vector<Int> refChan,
381 const Bool getSpectra,
382 const Bool getXPol,
383 const Bool getFeedPos)
384{
385 if (!cMSopen) {
386 return 1;
387 }
388
389 // Beam selection.
390 uInt nBeam = cBeams.nelements();
391 uInt nBeamSel = beamSel.nelements();
392 for (uInt ibeam = 0; ibeam < nBeam; ibeam++) {
393 if (ibeam < nBeamSel) {
394 cBeams(ibeam) = beamSel(ibeam);
395 } else {
396 cBeams(ibeam) = False;
397 }
398 }
399
400 uInt nIF = cIFs.nelements();
401 uInt maxNChan = 0;
402 for (uInt iIF = 0; iIF < nIF; iIF++) {
403 // IF selection.
404 if (iIF < IFsel.nelements()) {
405 cIFs(iIF) = IFsel(iIF);
406 } else {
407 cIFs(iIF) = False;
408 }
409
410 if (!cIFs(iIF)) continue;
411
412
413 // Channel selection.
414 if (iIF < startChan.nelements()) {
415 cStartChan(iIF) = startChan(iIF);
416
417 if (cStartChan(iIF) <= 0) {
418 cStartChan(iIF) += cNChan(iIF);
419 } else if (cStartChan(iIF) > Int(cNChan(iIF))) {
420 cStartChan(iIF) = cNChan(iIF);
421 }
422 }
423
424 if (iIF < endChan.nelements()) {
425 cEndChan(iIF) = endChan(iIF);
426
427 if (cEndChan(iIF) <= 0) {
428 cEndChan(iIF) += cNChan(iIF);
429 } else if (cEndChan(iIF) > Int(cNChan(iIF))) {
430 cEndChan(iIF) = cNChan(iIF);
431 }
432 }
433
434 if (iIF < refChan.nelements()) {
435 cRefChan(iIF) = refChan(iIF);
436 } else {
437 cRefChan(iIF) = cStartChan(iIF);
438 if (cStartChan(iIF) <= cEndChan(iIF)) {
439 cRefChan(iIF) += (cEndChan(iIF) - cStartChan(iIF) + 1)/2;
440 } else {
441 cRefChan(iIF) -= (cStartChan(iIF) - cEndChan(iIF) + 1)/2;
442 }
443 }
444
445 uInt nChan = abs(cEndChan(iIF) - cStartChan(iIF)) + 1;
446 if (maxNChan < nChan) {
447 maxNChan = nChan;
448 }
449
450 // Inverted Slices are not allowed.
451 Slice outPols;
452 Slice outChans(min(cStartChan(iIF),cEndChan(iIF))-1, nChan);
453 cDataSel(iIF) = Slicer(outPols, outChans);
454 }
455
456 // Get spectral data?
457 cGetSpectra = getSpectra;
458
459 // Get cross-polarization data?
460 cGetXPol = cGetXPol && getXPol;
461
462 // Get feed positions? (Not available.)
463 cGetFeedPos = False;
464
465 return maxNChan;
466}
467
468//---------------------------------------------------- PKSMS2reader::findRange
469
470// Find the range of the data in time and position.
471
472Int PKSMS2reader::findRange(
473 Int &nRow,
474 Int &nSel,
475 Vector<Double> &timeSpan,
476 Matrix<Double> &positions)
477{
478 if (!cMSopen) {
479 return 1;
480 }
481
482 nRow = cNRow;
483
484 // Find the number of rows selected.
485 nSel = 0;
486 Vector<Bool> sel(nRow);
487 for (Int irow = 0; irow < nRow; irow++) {
488 if ((sel(irow) = cBeams(cBeamNoCol(irow)) &&
489 cIFs(cDataDescIdCol(irow)))) {
490 nSel++;
491 }
492 }
493
494 // Find the time range (s).
495 timeSpan.resize(2);
496 timeSpan(0) = cTimeCol(0);
497 timeSpan(1) = cTimeCol(nRow-1);
498
499 // Retrieve positions for selected data.
500 Int isel = 0;
501 positions.resize(2,nSel);
502 for (Int irow = 0; irow < nRow; irow++) {
503 if (sel(irow)) {
504 Matrix<Double> pointingDir = cPointingCol(cFieldIdCol(irow));
505 positions.column(isel++) = pointingDir.column(0);
506 }
507 }
508
509 return 0;
510}
511
512//--------------------------------------------------------- PKSMS2reader::read
513
514// Read the next data record.
515
516Int PKSMS2reader::read(
517 Int &scanNo,
518 Int &cycleNo,
519 Double &mjd,
520 Double &interval,
521 String &fieldName,
522 String &srcName,
523 Vector<Double> &srcDir,
524 Vector<Double> &srcPM,
525 Double &srcVel,
526 String &obsMode,
527 Int &IFno,
528 Double &refFreq,
529 Double &bandwidth,
530 Double &freqInc,
531 Double &restFreq,
532 Vector<Float> &tcal,
533 String &tcalTime,
534 Float &azimuth,
535 Float &elevation,
536 Float &parAngle,
537 Float &focusAxi,
538 Float &focusTan,
539 Float &focusRot,
540 Float &temperature,
541 Float &pressure,
542 Float &humidity,
543 Float &windSpeed,
544 Float &windAz,
545 Int &refBeam,
546 Int &beamNo,
547 Vector<Double> &direction,
548 Vector<Double> &scanRate,
549 Vector<Float> &tsys,
550 Vector<Float> &sigma,
551 Vector<Float> &calFctr,
552 Matrix<Float> &baseLin,
553 Matrix<Float> &baseSub,
554 Matrix<Float> &spectra,
555 Matrix<uChar> &flagged,
556 Complex &xCalFctr,
557 Vector<Complex> &xPol)
558{
559 if (!cMSopen) {
560 return 1;
561 }
562
563 // Check for EOF.
564 if (cIdx >= cNRow) {
565 return -1;
566 }
567
568 // Find the next selected beam and IF.
569 Int ibeam;
570 Int iIF;
571 Int iDataDesc;
572
573 while (True) {
574 ibeam = cBeamNoCol(cIdx);
575 iDataDesc = cDataDescIdCol(cIdx);
576 iIF =cSpWinIdCol(iDataDesc);
577 if (cBeams(ibeam) && cIFs(iIF)) {
578 break;
579 }
580
581 // Check for EOF.
582 if (++cIdx >= cNRow) {
583 return -1;
584 }
585 }
586
587 // Renumerate scan no. Here still is 1-based
588 //scanNo = cScanNoCol(cIdx) - cScanNoCol(0) + 1;
589 scanNo = cScanNoCol(cIdx);
590
591 if (scanNo != cScanNo) {
592 // Start of new scan.
593 cScanNo = scanNo;
594 cCycleNo = 1;
595 cTime = cTimeCol(cIdx);
596 }
597
598 Double time = cTimeCol(cIdx);
599 mjd = time/86400.0;
600 interval = cIntervalCol(cIdx);
601
602 // Reconstruct the integration cycle number; due to small latencies the
603 // integration time is usually slightly less than the time between cycles,
604 // resetting cTime will prevent the difference from accumulating.
605 cCycleNo += nint((time - cTime)/interval);
606 cycleNo = cCycleNo;
607 cTime = time;
608
609 Int fieldId = cFieldIdCol(cIdx);
610 fieldName = cFieldNameCol(fieldId);
611
612 Int srcId = cSrcIdCol(fieldId);
613 //For source with multiple spectral window setting, this is not
614 // correct. Source name of srcId may not be at 'srcId'th row of SrcNameCol
615 //srcName = cSrcNameCol(srcId);
616 for (uInt irow = 0; irow < cSrcId2Col.nrow(); irow++) {
617 if (cSrcId2Col(irow) == srcId) {
618 srcName = cSrcNameCol(irow);
619 }
620 }
621
622 srcDir = cSrcDirCol(srcId);
623 srcPM = cSrcPMCol(srcId);
624
625 // Systemic velocity.
626 if (!cHaveSrcVel) {
627 srcVel = 0.0f;
628 } else {
629 srcVel = cSrcVelCol(srcId)(IPosition(1,0));
630 }
631
632 ROMSAntennaColumns antennaCols(cPKSMS.antenna());
633 String telescope = antennaCols.name()(0);
634 Bool cGBT = telescope.contains("GBT");
635 // Observation type.
636 Int stateId = cStateIdCol(cIdx);
637 if (stateId == -1) {
638 //obsMode = " ";
639 } else {
640 obsMode = cObsModeCol(stateId);
641 Bool sigState =cSigStateCol(stateId);
642 Bool refState =cRefStateCol(stateId);
643 //DEBUG
644 //cerr <<"stateid="<<stateId<<" obsmode="<<obsMode<<endl;
645 if (cGBT) {
646 // split the obsMode string and append a proper label
647 // (these are GBT specific)
648 int epos = obsMode.find_first_of(':');
649 int nextpos = obsMode.find_first_of(':',epos+1);
650 string obsMode1 = obsMode.substr(0,epos);
651 string obsMode2 = obsMode.substr(epos+1,nextpos-epos-1);
652
653 //cerr <<"obsMode2= "<<obsMode2<<endl;
654 if (!srcName.contains("_ps")
655 &&!srcName.contains("_psr")
656 &&!srcName.contains("_nod")
657 &&!srcName.contains("_fs")
658 &&!srcName.contains("_fsr")) {
659 // if Nod mode observation , append '_nod'
660 if (obsMode1 == "Nod") {
661 srcName.append("_nod");
662 } else if (obsMode1 == "OffOn") {
663 // for GBT position switch observations (OffOn or OnOff)
664 if (obsMode2 == "PSWITCHON") srcName.append("_ps");
665 if (obsMode2 == "PSWITCHOFF") srcName.append("_psr");
666 } else {
667 if (obsMode2 == "FSWITCH") {
668 // for GBT frequency switch mode
669 if (sigState) srcName.append("_fs");
670 if (refState) srcName.append("_fsr");
671 }
672 }
673 }
674 }
675 }
676 // CAL state
677 // this should be apply just for GBT data?
678 Double Cal;
679 if (stateId==-1) {
680 Cal = 0;
681 } else {
682 Cal = cCalCol(stateId);
683 }
684 if (cGBT) {
685 if (Cal > 0 && !srcName.contains("_calon")) {
686 srcName.append("_calon");
687 }
688 }
689
690 IFno = iIF + 1;
691 Int nChan = abs(cEndChan(iIF) - cStartChan(iIF)) + 1;
692
693 // Minimal handling on continuum data.
694 Vector<Double> chanFreq = cChanFreqCol(iIF);
695 if (nChan == 1) {
696 cout << "The input is continuum data. "<< endl;
697 freqInc = chanFreq(0);
698 refFreq = chanFreq(0);
699 restFreq = 0.0f;
700 } else {
701 if (cStartChan(iIF) <= cEndChan(iIF)) {
702 freqInc = chanFreq(1) - chanFreq(0);
703 } else {
704 freqInc = chanFreq(0) - chanFreq(1);
705 }
706 refFreq = chanFreq(cRefChan(iIF)-1);
707 restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0));
708 }
709 bandwidth = abs(freqInc * nChan);
710
711 tcal.resize(cNPol(iIF));
712 tcal = 0.0f;
713 tcalTime = "";
714 //azimuth = 0.0f;
715 //elevation = 0.0f;
716 parAngle = 0.0f;
717 focusAxi = 0.0f;
718 focusTan = 0.0f;
719 focusRot = 0.0f;
720
721 // Find the appropriate entry in the WEATHER subtable.
722 Vector<Double> wTimes = cWeatherTimeCol.getColumn();
723 Int weatherIdx;
724 for (weatherIdx = wTimes.nelements()-1; weatherIdx >= 0; weatherIdx--) {
725 if (cWeatherTimeCol(weatherIdx) <= time) {
726 break;
727 }
728 }
729
730 if (weatherIdx < 0) {
731 // No appropriate WEATHER entry.
732 pressure = 0.0f;
733 humidity = 0.0f;
734 temperature = 0.0f;
735 } else {
736 pressure = cPressureCol(weatherIdx);
737 humidity = cHumidityCol(weatherIdx);
738 temperature = cTemperatureCol(weatherIdx);
739 }
740
741 windSpeed = 0.0f;
742 windAz = 0.0f;
743
744 refBeam = 0;
745 beamNo = ibeam + 1;
746
747 //Matrix<Double> pointingDir = cPointingCol(fieldId);
748 //pointingDir = cPointingCol(fieldId);
749 //direction = pointingDir.column(0);
750 //uInt ncols = pointingDir.ncolumn();
751 //if (ncols == 1) {
752 // scanRate = 0.0f;
753 //} else {
754 // scanRate = pointingDir.column(1);
755 //}
756
757 // Get direction from FIELD table
758 // here, assume direction to be the field direction not pointing
759 Matrix<Double> delayDir = cFieldDelayDirCol(fieldId);
760 direction = delayDir.column(0);
761 uInt ncols = delayDir.ncolumn();
762 if (ncols == 1) {
763 scanRate = 0.0f;
764 } else {
765 scanRate = delayDir.column(1);
766 }
767
768 // caluculate azimuth and elevation
769 // first, get the reference frame
770 MVPosition mvpos(antennaCols.position()(0));
771 MPosition mp(mvpos);
772 Quantum<Double> qt(time,"s");
773 MVEpoch mvt(qt);
774 MEpoch me(mvt);
775 MeasFrame frame(mp, me);
776 //
777 ROMSFieldColumns fldCols(cPKSMS.field());
778 Vector<MDirection> vmd(1);
779 MDirection md;
780 fldCols.delayDirMeasCol().get(fieldId,vmd);
781 md = vmd[0];
782 //Vector<Double> dircheck = md.getAngle("rad").getValue();
783 //cerr<<"dircheck="<<dircheck<<endl;
784
785 Vector<Double> azel =
786 MDirection::Convert(md, MDirection::Ref(MDirection::AZEL,
787 frame)
788 )().getAngle("rad").getValue();
789 //cerr<<"azel="<<azel<<endl;
790 azimuth = azel[0];
791 elevation = azel[1];
792
793 // Get Tsys assuming that entries in the SYSCAL table match the main table.
794 if (cHaveTsys) {
795 Int nTsysColRow = cTsysCol.nrow();
796 if (nTsysColRow != cNRow) {
797 cHaveTsys=0;
798 }
799 }
800 if (cHaveTsys) {
801 cTsysCol.get(cIdx, tsys, True);
802 } else {
803 Int numReceptor;
804 cNumReceptorCol.get(0, numReceptor);
805 tsys.resize(numReceptor);
806 tsys = 1.0f;
807 }
808 cSigmaCol.get(cIdx, sigma, True);
809
810 //get Tcal if available
811 if (cHaveTcal) {
812 Int nTcalColRow = cTcalCol.nrow();
813 uInt nBeam = cBeams.nelements();
814 uInt nIF = cIFs.nelements();
815 uInt nrws = nBeam * nIF;
816 if (nTcalColRow > 0) {
817 // find tcal match with the data with the data time stamp
818 Double mjds = mjd*(24*3600);
819 Double dtcalTime;
820 if ( mjd > lastmjd || cIdx==0 ) {
821 //Table tmptab = cSysCalTab(near(cSysCalTab.col("TIME"),mjds));
822 tmptab = cSysCalTab(near(cSysCalTab.col("TIME"),mjds), nrws);
823 //DEBUG
824 //if (cIdx == 0) {
825 // cerr<<"inital table retrieved"<<endl;
826 //}
827
828 }
829
830 if (nBeam == 1) {
831 tmptab2 = tmptab( tmptab.col("SPECTRAL_WINDOW_ID") == iIF, 1);
832 } else {
833 tmptab2 = tmptab( tmptab.col("SPECTRAL_WINDOW_ID") == iIF &&
834 tmptab.col("FEED_ID") == ibeam , 1);
835 }
836 //cerr<<"first subtab rows="<<tmptab.nrow()<<endl;
837 int syscalrow = tmptab2.nrow();
838 ROArrayColumn<Float> tcalCol(tmptab2, "TCAL");
839 ROScalarColumn<Double> tcalTimeCol(tmptab2, "TIME");
840 if (syscalrow==0) {
841 cerr<<"Cannot find any matching Tcal at/near the data timestamp."
842 << " Set Tcal=0.0"<<endl;
843 } else {
844 tcalCol.get(0, tcal);
845 tcalTimeCol.get(0,dtcalTime);
846 tcalTime = MVTime(dtcalTime/(24*3600)).string(MVTime::YMD);
847 //DEBUG
848 //cerr<<"cIdx:"<<cIdx<<" tcal="<<tcal<<" tcalTime="<<tcalTime<<endl;
849 tmptab.markForDelete();
850 tmptab2.markForDelete();
851 }
852 }
853 lastmjd = mjd;
854 }
855
856 // Calibration factors (if available).
857 calFctr.resize(cNPol(iIF));
858 if (cHaveCalFctr) {
859 cCalFctrCol.get(cIdx, calFctr);
860 } else {
861 calFctr = 0.0f;
862 }
863
864 // Baseline parameters (if available).
865 if (cHaveBaseLin) {
866 baseLin.resize(2,cNPol(iIF));
867 cBaseLinCol.get(cIdx, baseLin);
868
869 baseSub.resize(9,cNPol(iIF));
870 cBaseSubCol.get(cIdx, baseSub);
871
872 } else {
873 baseLin.resize(0,0);
874 baseSub.resize(0,0);
875 }
876
877 // Get spectral data.
878 if (cGetSpectra) {
879 Matrix<Float> tmpData;
880 Matrix<Bool> tmpFlag;
881 cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
882 cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
883
884 // Transpose spectra.
885 Int nPol = tmpData.nrow();
886 spectra.resize(nChan, nPol);
887 flagged.resize(nChan, nPol);
888 if (cEndChan(iIF) >= cStartChan(iIF)) {
889 // Simple transposition.
890 for (Int ipol = 0; ipol < nPol; ipol++) {
891 for (Int ichan = 0; ichan < nChan; ichan++) {
892 spectra(ichan,ipol) = tmpData(ipol,ichan);
893 flagged(ichan,ipol) = tmpFlag(ipol,ichan);
894 }
895 }
896
897 } else {
898 // Transpose with inversion.
899 Int jchan = nChan - 1;
900 for (Int ipol = 0; ipol < nPol; ipol++) {
901 for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
902 spectra(ichan,ipol) = tmpData(ipol,jchan);
903 flagged(ichan,ipol) = tmpFlag(ipol,jchan);
904 }
905 }
906 }
907 }
908
909 // Get cross-polarization data.
910 if (cGetXPol) {
911 if (cHaveXCalFctr) {
912 cXCalFctrCol.get(cIdx, xCalFctr);
913 } else {
914 xCalFctr = Complex(0.0f, 0.0f);
915 }
916
917 cDataCol.get(cIdx, xPol, True);
918
919 if (cEndChan(iIF) < cStartChan(iIF)) {
920 Complex ctmp;
921 Int jchan = nChan - 1;
922 for (Int ichan = 0; ichan < nChan/2; ichan++, jchan--) {
923 ctmp = xPol(ichan);
924 xPol(ichan) = xPol(jchan);
925 xPol(jchan) = ctmp;
926 }
927 }
928 }
929
930 cIdx++;
931
932 return 0;
933}
934
935//--------------------------------------------------------- PKSMS2reader::read
936
937// Read the next data record, just the basics.
938
939Int PKSMS2reader::read(
940 Int &IFno,
941 Vector<Float> &tsys,
942 Vector<Float> &calFctr,
943 Matrix<Float> &baseLin,
944 Matrix<Float> &baseSub,
945 Matrix<Float> &spectra,
946 Matrix<uChar> &flagged)
947{
948 if (!cMSopen) {
949 return 1;
950 }
951
952 // Check for EOF.
953 if (cIdx >= cNRow) {
954 return -1;
955 }
956
957 // Find the next selected beam and IF.
958 Int ibeam;
959 Int iIF;
960 Int iDataDesc;
961 while (True) {
962 ibeam = cBeamNoCol(cIdx);
963 //iIF = cDataDescIdCol(cIdx);
964 iDataDesc = cDataDescIdCol(cIdx);
965 iIF = cSpWinIdCol(iDataDesc);
966 if (cBeams(ibeam) && cIFs(iIF)) {
967 break;
968 }
969
970 // Check for EOF.
971 if (++cIdx >= cNRow) {
972 return -1;
973 }
974 }
975
976 IFno = iIF + 1;
977 // Get Tsys assuming that entries in the SYSCAL table match the main table.
978 cTsysCol.get(cIdx, tsys, True);
979
980 // Calibration factors (if available).
981 if (cHaveCalFctr) {
982 cCalFctrCol.get(cIdx, calFctr, True);
983 } else {
984 calFctr.resize(cNPol(iIF));
985 calFctr = 0.0f;
986 }
987
988 // Baseline parameters (if available).
989 if (cHaveBaseLin) {
990 baseLin.resize(2,cNPol(iIF));
991 cBaseLinCol.get(cIdx, baseLin);
992
993 baseSub.resize(9,cNPol(iIF));
994 cBaseSubCol.get(cIdx, baseSub);
995
996 } else {
997 baseLin.resize(0,0);
998 baseSub.resize(0,0);
999 }
1000
1001 if (cGetSpectra) {
1002 // Get spectral data.
1003 Matrix<Float> tmpData;
1004 Matrix<Bool> tmpFlag;
1005 cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
1006 cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
1007
1008 // Transpose spectra.
1009 Int nChan = tmpData.ncolumn();
1010 Int nPol = tmpData.nrow();
1011 spectra.resize(nChan, nPol);
1012 flagged.resize(nChan, nPol);
1013 if (cEndChan(iIF) >= cStartChan(iIF)) {
1014 // Simple transposition.
1015 for (Int ipol = 0; ipol < nPol; ipol++) {
1016 for (Int ichan = 0; ichan < nChan; ichan++) {
1017 spectra(ichan,ipol) = tmpData(ipol,ichan);
1018 flagged(ichan,ipol) = tmpFlag(ipol,ichan);
1019 }
1020 }
1021
1022 } else {
1023 // Transpose with inversion.
1024 Int jchan = nChan - 1;
1025 for (Int ipol = 0; ipol < nPol; ipol++) {
1026 for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
1027 spectra(ichan,ipol) = tmpData(ipol,jchan);
1028 flagged(ichan,ipol) = tmpFlag(ipol,jchan);
1029 }
1030 }
1031 }
1032 }
1033
1034 cIdx++;
1035
1036 return 0;
1037}
1038
1039//-------------------------------------------------------- PKSMS2reader::close
1040
1041// Close the MS.
1042
1043void PKSMS2reader::close()
1044{
1045 cPKSMS = MeasurementSet();
1046 cMSopen = False;
1047}
Note: See TracBrowser for help on using the repository browser.