source: branches/alma/external/atnf/PKSIO/PKSFITSreader.cc@ 1393

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

Changes to handle GBT MS data
File size: 14.8 KB
Line 
1//#---------------------------------------------------------------------------
2//# PKSFITSreader.cc: Class to read Parkes multibeam data from a FITS file.
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/02, Mark Calabretta, ATNF
31//#---------------------------------------------------------------------------
32
33#include <atnf/PKSIO/MBFITSreader.h>
34#include <atnf/PKSIO/SDFITSreader.h>
35#include <atnf/PKSIO/PKSFITSreader.h>
36#include <atnf/PKSIO/PKSMBrecord.h>
37
38#include <casa/Arrays/Array.h>
39#include <casa/BasicMath/Math.h>
40#include <casa/Quanta/MVTime.h>
41
42#include <casa/stdio.h>
43
44//----------------------------------------------- PKSFITSreader::PKSFITSreader
45
46// Constructor sets the method of position interpolation.
47
48PKSFITSreader::PKSFITSreader(
49 const String fitsType,
50 const Int retry,
51 const Bool interpolate)
52{
53 cMBrec.setNIFs(1);
54
55 if (fitsType == "SDFITS") {
56 cReader = new SDFITSreader();
57 } else {
58 cReader = new MBFITSreader(retry, interpolate ? 1 : 0);
59 }
60}
61
62//---------------------------------------------- PKSFITSreader::~PKSFITSreader
63
64// Destructor.
65
66PKSFITSreader::~PKSFITSreader()
67{
68 close();
69 delete cReader;
70}
71
72//-------------------------------------------------------- PKSFITSreader::open
73
74// Open the FITS file for reading.
75
76Int PKSFITSreader::open(
77 const String fitsName,
78 Vector<Bool> &beams,
79 Vector<Bool> &IFs,
80 Vector<uInt> &nChan,
81 Vector<uInt> &nPol,
82 Vector<Bool> &haveXPol,
83 Bool &haveBase,
84 Bool &haveSpectra)
85{
86 int extraSysCal, haveBase_, *haveXPol_, haveSpectra_, nBeam, *nChan_,
87 nIF, *nPol_, status;
88 if ((status = cReader->open((char *)fitsName.chars(), nBeam, cBeams, nIF,
89 cIFs, nChan_, nPol_, haveXPol_, haveBase_,
90 haveSpectra_, extraSysCal))) {
91 return status;
92 }
93
94 // Beams present in data.
95 beams.resize(nBeam);
96 for (Int ibeam = 0; ibeam < nBeam; ibeam++) {
97 beams(ibeam) = cBeams[ibeam];
98 }
99
100 // IFs, channels, and polarizations present in data.
101 IFs.resize(nIF);
102 nChan.resize(nIF);
103 nPol.resize(nIF);
104 haveXPol.resize(nIF);
105
106 for (Int iIF = 0; iIF < nIF; iIF++) {
107 IFs(iIF) = cIFs[iIF];
108 nChan(iIF) = nChan_[iIF];
109 nPol(iIF) = nPol_[iIF];
110
111 // Cross-polarization data present?
112 haveXPol(iIF) = haveXPol_[iIF];
113 }
114
115 cNBeam = nBeam;
116 cNIF = nIF;
117 cNChan.assign(nChan);
118 cNPol.assign(nPol);
119 cHaveXPol.assign(haveXPol);
120
121 // Baseline parameters present?
122 haveBase = haveBase_;
123
124 // Spectral data present?
125 haveSpectra = haveSpectra_;
126
127 return 0;
128}
129
130//--------------------------------------------------- PKSFITSreader::getHeader
131
132// Get parameters describing the data.
133
134Int PKSFITSreader::getHeader(
135 String &observer,
136 String &project,
137 String &antName,
138 Vector<Double> &antPosition,
139 String &obsType,
140 Float &equinox,
141 String &dopplerFrame,
142 Double &mjd,
143 Double &refFreq,
144 Double &bandwidth,
145 String &fluxunit)
146{
147 char datobs[32], dopplerFrame_[32], observer_[32], obsType_[32],
148 project_[32], radecsys[32], telescope[32];
149 float equinox_;
150 double antPos[3], utc;
151
152 if (cReader->getHeader(observer_, project_, telescope, antPos, obsType_,
153 equinox_, radecsys, dopplerFrame_, datobs, utc,
154 refFreq, bandwidth)) {
155 return 1;
156 }
157
158 fluxunit = "";
159 observer = trim(observer_);
160 project = trim(project_);
161 antName = trim(telescope);
162 antPosition.resize(3);
163 antPosition(0) = antPos[0];
164 antPosition(1) = antPos[1];
165 antPosition(2) = antPos[2];
166 obsType = trim(obsType_);
167 equinox = equinox_;
168 dopplerFrame = trim(dopplerFrame_);
169
170 // Get UTC in MJD form.
171 Int day, month, year;
172 sscanf(datobs, "%4d-%2d-%2d", &year, &month, &day);
173 MVTime date(year, month, Double(day));
174 mjd = date.day() + utc/86400.0;
175
176 return 0;
177}
178
179//------------------------------------------------- PKSFITSreader::getFreqInfo
180
181// Get frequency parameters for each IF.
182
183Int PKSFITSreader::getFreqInfo(
184 Vector<Double> &startFreq,
185 Vector<Double> &endFreq)
186{
187 int nIF;
188 double *startfreq, *endfreq;
189
190 Int status;
191 if (!(status = cReader->getFreqInfo(nIF, startfreq, endfreq))) {
192 startFreq.takeStorage(IPosition(1,nIF), startfreq, TAKE_OVER);
193 endFreq.takeStorage(IPosition(1,nIF), endfreq, TAKE_OVER);
194 }
195
196 return status;
197}
198
199//------------------------------------------------------ PKSFITSreader::select
200
201// Set data selection by beam number and channel.
202
203uInt PKSFITSreader::select(
204 const Vector<Bool> beamSel,
205 const Vector<Bool> IFsel,
206 const Vector<Int> startChan,
207 const Vector<Int> endChan,
208 const Vector<Int> refChan,
209 const Bool getSpectra,
210 const Bool getXPol,
211 const Bool getFeedPos)
212{
213 // Apply beam selection.
214 uInt nBeamSel = beamSel.nelements();
215 for (uInt ibeam = 0; ibeam < cNBeam; ibeam++) {
216 // This modifies FITSreader::cBeams[].
217 if (ibeam < nBeamSel) {
218 cBeams[ibeam] = cBeams[ibeam] && beamSel(ibeam);
219 } else {
220 cBeams[ibeam] = 0;
221 }
222 }
223
224 // Apply IF selection.
225 int *end = new int[cNIF];
226 int *ref = new int[cNIF];
227 int *start = new int[cNIF];
228
229 for (uInt iIF = 0; iIF < cNIF; iIF++) {
230 // This modifies FITSreader::cIFs[].
231 if (iIF < IFsel.nelements()) {
232 cIFs[iIF] = cIFs[iIF] && IFsel(iIF);
233 } else {
234 cIFs[iIF] = 0;
235 }
236
237 if (cIFs[iIF]) {
238 if (iIF < startChan.nelements()) {
239 start[iIF] = startChan(iIF);
240
241 if (start[iIF] <= 0) {
242 start[iIF] += cNChan(iIF);
243 } else if (start[iIF] > Int(cNChan(iIF))) {
244 start[iIF] = cNChan(iIF);
245 }
246
247 } else {
248 start[iIF] = 1;
249 }
250
251 if (iIF < endChan.nelements()) {
252 end[iIF] = endChan(iIF);
253
254 if (end[iIF] <= 0) {
255 end[iIF] += cNChan(iIF);
256 } else if (end[iIF] > Int(cNChan(iIF))) {
257 end[iIF] = cNChan(iIF);
258 }
259
260 } else {
261 end[iIF] = cNChan(iIF);
262 }
263
264 if (iIF < refChan.nelements()) {
265 ref[iIF] = refChan(iIF);
266 } else {
267 if (start[iIF] <= end[iIF]) {
268 ref[iIF] = start[iIF] + (end[iIF] - start[iIF] + 1)/2;
269 } else {
270 ref[iIF] = start[iIF] - (start[iIF] - end[iIF] + 1)/2;
271 }
272 }
273 }
274 }
275
276 cGetSpectra = getSpectra;
277 cGetXPol = getXPol;
278 cGetFeedPos = getFeedPos;
279
280 uInt maxNChan = cReader->select(start, end, ref, cGetSpectra, cGetXPol,
281 cGetFeedPos);
282
283 delete [] end;
284 delete [] ref;
285 delete [] start;
286
287 return maxNChan;
288}
289
290//--------------------------------------------------- PKSFITSreader::findRange
291
292// Read the TIME column.
293
294Int PKSFITSreader::findRange(
295 Int &nRow,
296 Int &nSel,
297 Vector<Double> &timeSpan,
298 Matrix<Double> &positions)
299{
300 char dateSpan[2][32];
301 double utcSpan[2];
302 double* posns;
303
304 Int status;
305 if (!(status = cReader->findRange(nRow, nSel, dateSpan, utcSpan, posns))) {
306 timeSpan.resize(2);
307
308 Int day, month, year;
309 sscanf(dateSpan[0], "%4d-%2d-%2d", &year, &month, &day);
310 timeSpan(0) = MVTime(year, month, Double(day)).second() + utcSpan[0];
311 sscanf(dateSpan[1], "%4d-%2d-%2d", &year, &month, &day);
312 timeSpan(1) = MVTime(year, month, Double(day)).second() + utcSpan[1];
313
314 positions.takeStorage(IPosition(2,2,nSel), posns, TAKE_OVER);
315 }
316
317 return status;
318}
319
320//-------------------------------------------------------- PKSFITSreader::read
321
322// Read the next data record.
323
324Int PKSFITSreader::read(
325 Int &scanNo,
326 Int &cycleNo,
327 Double &mjd,
328 Double &interval,
329 String &fieldName,
330 String &srcName,
331 Vector<Double> &srcDir,
332 Vector<Double> &srcPM,
333 Double &srcVel,
334 String &obsType,
335 Int &IFno,
336 Double &refFreq,
337 Double &bandwidth,
338 Double &freqInc,
339 Double &restFreq,
340 Vector<Float> &tcal,
341 String &tcalTime,
342 Float &azimuth,
343 Float &elevation,
344 Float &parAngle,
345 Float &focusAxi,
346 Float &focusTan,
347 Float &focusRot,
348 Float &temperature,
349 Float &pressure,
350 Float &humidity,
351 Float &windSpeed,
352 Float &windAz,
353 Int &refBeam,
354 Int &beamNo,
355 Vector<Double> &direction,
356 Vector<Double> &scanRate,
357 Vector<Float> &tsys,
358 Vector<Float> &sigma,
359 Vector<Float> &calFctr,
360 Matrix<Float> &baseLin,
361 Matrix<Float> &baseSub,
362 Matrix<Float> &spectra,
363 Matrix<uChar> &flagged,
364 Complex &xCalFctr,
365 Vector<Complex> &xPol)
366{
367 Int status;
368
369 if ((status = cReader->read(cMBrec))) {
370 if (status != -1) {
371 status = 1;
372 }
373
374 return status;
375 }
376
377
378 uInt nChan = cMBrec.nChan[0];
379 uInt nPol = cMBrec.nPol[0];
380
381 scanNo = cMBrec.scanNo;
382 cycleNo = cMBrec.cycleNo;
383
384 // Extract MJD.
385 Int day, month, year;
386 sscanf(cMBrec.datobs, "%4d-%2d-%2d", &year, &month, &day);
387 mjd = MVTime(year, month, Double(day)).day() + cMBrec.utc/86400.0;
388
389 interval = cMBrec.exposure;
390
391 fieldName = trim(cMBrec.srcName);
392 srcName = fieldName;
393 srcDir(0) = cMBrec.srcRA;
394 srcDir(1) = cMBrec.srcDec;
395 srcPM(0) = 0.0;
396 srcPM(1) = 0.0;
397 srcVel = 0.0;
398 obsType = trim(cMBrec.obsType);
399
400 IFno = cMBrec.IFno[0];
401 Double chanWidth = fabs(cMBrec.fqDelt[0]);
402 refFreq = cMBrec.fqRefVal[0];
403 bandwidth = chanWidth * nChan;
404 freqInc = cMBrec.fqDelt[0];
405 restFreq = cMBrec.restFreq;
406
407 tcal.resize(nPol);
408 for (uInt ipol = 0; ipol < nPol; ipol++) {
409 tcal(ipol) = cMBrec.tcal[0][ipol];
410 }
411 tcalTime = trim(cMBrec.tcalTime);
412 azimuth = cMBrec.azimuth;
413 elevation = cMBrec.elevation;
414 parAngle = cMBrec.parAngle;
415 focusAxi = cMBrec.focusAxi;
416 focusTan = cMBrec.focusTan;
417 focusRot = cMBrec.focusRot;
418
419 temperature = cMBrec.temp;
420 pressure = cMBrec.pressure;
421 humidity = cMBrec.humidity;
422 windSpeed = cMBrec.windSpeed;
423 windAz = cMBrec.windAz;
424
425 refBeam = cMBrec.refBeam;
426 beamNo = cMBrec.beamNo;
427
428 direction(0) = cMBrec.ra;
429 direction(1) = cMBrec.dec;
430 scanRate(0) = cMBrec.raRate;
431 scanRate(1) = cMBrec.decRate;
432
433 tsys.resize(nPol);
434 sigma.resize(nPol);
435 calFctr.resize(nPol);
436 for (uInt ipol = 0; ipol < nPol; ipol++) {
437 tsys(ipol) = cMBrec.tsys[0][ipol];
438 sigma(ipol) = tsys(ipol) / 0.81 / sqrt(interval * chanWidth);
439 calFctr(ipol) = cMBrec.calfctr[0][ipol];
440 }
441
442 if (cMBrec.haveBase) {
443 baseLin.resize(2,nPol);
444 baseSub.resize(9,nPol);
445
446 for (uInt ipol = 0; ipol < nPol; ipol++) {
447 baseLin(0,ipol) = cMBrec.baseLin[0][ipol][0];
448 baseLin(1,ipol) = cMBrec.baseLin[0][ipol][1];
449
450 for (uInt j = 0; j < 9; j++) {
451 baseSub(j,ipol) = cMBrec.baseSub[0][ipol][j];
452 }
453 }
454
455 } else {
456 baseLin.resize(0,0);
457 baseSub.resize(0,0);
458 }
459
460 if (cGetSpectra && cMBrec.haveSpectra) {
461 spectra.resize(nChan,nPol);
462 spectra.takeStorage(IPosition(2,nChan,nPol), cMBrec.spectra[0], SHARE);
463
464 flagged.resize(nChan,nPol);
465 flagged.takeStorage(IPosition(2,nChan,nPol), cMBrec.flagged[0], SHARE);
466
467 } else {
468 spectra.resize(0,0);
469 flagged.resize(0,0);
470 }
471
472 if (cGetXPol) {
473 xCalFctr = Complex(cMBrec.xcalfctr[0][0], cMBrec.xcalfctr[0][1]);
474 xPol.resize(nChan);
475 xPol.takeStorage(IPosition(1,nChan), (Complex *)cMBrec.xpol[0], SHARE);
476 }
477
478 return 0;
479}
480
481//-------------------------------------------------------- PKSFITSreader::read
482
483// Read the next data record, just the basics.
484
485Int PKSFITSreader::read(
486 Int &IFno,
487 Vector<Float> &tsys,
488 Vector<Float> &calFctr,
489 Matrix<Float> &baseLin,
490 Matrix<Float> &baseSub,
491 Matrix<Float> &spectra,
492 Matrix<uChar> &flagged)
493{
494 Int status;
495
496 if ((status = cReader->read(cMBrec))) {
497 if (status != -1) {
498 status = 1;
499 }
500
501 return status;
502 }
503
504 IFno = cMBrec.IFno[0];
505
506 uInt nChan = cMBrec.nChan[0];
507 uInt nPol = cMBrec.nPol[0];
508
509 tsys.resize(nPol);
510 calFctr.resize(nPol);
511 for (uInt ipol = 0; ipol < nPol; ipol++) {
512 tsys(ipol) = cMBrec.tsys[0][ipol];
513 calFctr(ipol) = cMBrec.calfctr[0][ipol];
514 }
515
516 if (cMBrec.haveBase) {
517 baseLin.resize(2,nPol);
518 baseSub.resize(9,nPol);
519
520 for (uInt ipol = 0; ipol < nPol; ipol++) {
521 baseLin(0,ipol) = cMBrec.baseLin[0][ipol][0];
522 baseLin(1,ipol) = cMBrec.baseLin[0][ipol][1];
523
524 for (uInt j = 0; j < 9; j++) {
525 baseSub(j,ipol) = cMBrec.baseSub[0][ipol][j];
526 }
527 }
528
529 } else {
530 baseLin.resize(0,0);
531 baseSub.resize(0,0);
532 }
533
534 if (cGetSpectra && cMBrec.haveSpectra) {
535 spectra.resize(nChan,nPol);
536 spectra.takeStorage(IPosition(2,nChan,nPol), cMBrec.spectra[0], SHARE);
537
538 flagged.resize(nChan,nPol);
539 flagged.takeStorage(IPosition(2,nChan,nPol), cMBrec.flagged[0], SHARE);
540
541 } else {
542 spectra.resize(0,0);
543 flagged.resize(0,0);
544 }
545
546 return 0;
547}
548
549//------------------------------------------------------- PKSFITSreader::close
550
551// Close the FITS file.
552
553void PKSFITSreader::close()
554{
555 cReader->close();
556}
557
558//-------------------------------------------------------- PKSFITSreader::trim
559
560// Trim trailing blanks from a null-terminated character string.
561
562char* PKSFITSreader::trim(char *string)
563{
564 int j = 0, k = 0;
565 while (string[j] != '\0') {
566 if (string[j++] != ' ') {
567 k = j;
568 }
569 }
570
571 string[k] = '\0';
572
573 return string;
574}
Note: See TracBrowser for help on using the repository browser.