source: branches/hpc34/external/atnf/PKSIO/PKSFITSreader.cc@ 2874

Last change on this file since 2874 was 1720, checked in by Malte Marquarding, 15 years ago

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