source: trunk/external/atnf/PKSIO/PKSFITSreader.cc@ 2675

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

Update from livedata CVS repository

File size: 13.6 KB
RevLine 
[1720]1//#---------------------------------------------------------------------------
[1325]2//# PKSFITSreader.cc: Class to read Parkes multibeam data from a FITS file.
3//#---------------------------------------------------------------------------
[1720]4//# livedata - processing pipeline for single-dish, multibeam spectral data.
5//# Copyright (C) 2000-2009, Australia Telescope National Facility, CSIRO
[1325]6//#
[1720]7//# This file is part of livedata.
[1325]8//#
[1720]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.
[1325]13//#
[1720]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.
[1325]18//#
[1720]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/>.
[1325]21//#
[1720]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 $
[1325]32//#---------------------------------------------------------------------------
33//# Original: 2000/08/02, Mark Calabretta, ATNF
34//#---------------------------------------------------------------------------
35
[1452]36#include <atnf/PKSIO/PKSmsg.h>
[1325]37#include <atnf/PKSIO/MBFITSreader.h>
38#include <atnf/PKSIO/SDFITSreader.h>
39#include <atnf/PKSIO/PKSFITSreader.h>
[1452]40#include <atnf/PKSIO/PKSrecord.h>
[1325]41
[1452]42#include <casa/stdio.h>
[1325]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{
[1452]56 cMBrec.setNIFs(1);
[1325]57
58 if (fitsType == "SDFITS") {
59 cReader = new SDFITSreader();
60 } else {
61 cReader = new MBFITSreader(retry, interpolate ? 1 : 0);
62 }
[1452]63
64 // By default, messages are written to stderr.
65 initMsg();
[1325]66}
67
68//---------------------------------------------- PKSFITSreader::~PKSFITSreader
69
70// Destructor.
71
72PKSFITSreader::~PKSFITSreader()
73{
74 close();
75 delete cReader;
76}
77
[1452]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
[1325]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{
[1452]105 clearMsg();
106
[1325]107 int extraSysCal, haveBase_, *haveXPol_, haveSpectra_, nBeam, *nChan_,
108 nIF, *nPol_, status;
[1452]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) {
[1325]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,
[1399]164 String &bunit,
[1325]165 Float &equinox,
166 String &dopplerFrame,
167 Double &mjd,
168 Double &refFreq,
169 Double &bandwidth)
170{
[1399]171 char bunit_[32], datobs[32], dopplerFrame_[32], observer_[32],
172 obsType_[32], project_[32], radecsys[32], telescope[32];
[1452]173 int status;
[1325]174 float equinox_;
175 double antPos[3], utc;
176
[1452]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) {
[1325]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_);
[1399]194 bunit = trim(bunit_);
[1325]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
[1452]218 Int status = cReader->getFreqInfo(nIF, startfreq, endfreq);
219
220 logMsg(cReader->getMsg());
221 cReader->clearMsg();
222 if (!status) {
[1325]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,
[1452]242 const Int coordSys)
[1325]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;
[1452]309 cCoordSys = coordSys;
[1325]310
311 uInt maxNChan = cReader->select(start, end, ref, cGetSpectra, cGetXPol,
[1452]312 cCoordSys);
313 logMsg(cReader->getMsg());
314 cReader->clearMsg();
[1325]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
[1452]337 Int status = cReader->findRange(nRow, nSel, dateSpan, utcSpan, posns);
338 logMsg(cReader->getMsg());
339 cReader->clearMsg();
340
341 if (!status) {
[1325]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
[1452]360Int PKSFITSreader::read(PKSrecord &pksrec)
[1325]361{
[1452]362 Int status = cReader->read(cMBrec);
363 logMsg(cReader->getMsg());
364 cReader->clearMsg();
[1325]365
[1452]366 if (status) {
[1325]367 if (status != -1) {
368 status = 1;
369 }
370
371 return status;
372 }
373
374
[1452]375 uInt nChan = cMBrec.nChan[0];
376 uInt nPol = cMBrec.nPol[0];
[1325]377
[1452]378 pksrec.scanNo = cMBrec.scanNo;
379 pksrec.cycleNo = cMBrec.cycleNo;
[1325]380
381 // Extract MJD.
382 Int day, month, year;
[1452]383 sscanf(cMBrec.datobs, "%4d-%2d-%2d", &year, &month, &day);
384 pksrec.mjd = MVTime(year, month, Double(day)).day() + cMBrec.utc/86400.0;
[1325]385
[1452]386 pksrec.interval = cMBrec.exposure;
[1325]387
[1452]388 pksrec.fieldName = trim(cMBrec.srcName);
389 pksrec.srcName = pksrec.fieldName;
[1325]390
[1452]391 pksrec.srcDir.resize(2);
392 pksrec.srcDir(0) = cMBrec.srcRA;
393 pksrec.srcDir(1) = cMBrec.srcDec;
[1325]394
[1452]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);
[1427]400
[1452]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;
[1427]407
[1452]408 pksrec.tcal.resize(nPol);
[1325]409 for (uInt ipol = 0; ipol < nPol; ipol++) {
[1452]410 pksrec.tcal(ipol) = cMBrec.tcal[0][ipol];
[1325]411 }
[1452]412 pksrec.tcalTime = trim(cMBrec.tcalTime);
413 pksrec.azimuth = cMBrec.azimuth;
414 pksrec.elevation = cMBrec.elevation;
415 pksrec.parAngle = cMBrec.parAngle;
[1325]416
[1452]417 pksrec.focusAxi = cMBrec.focusAxi;
418 pksrec.focusTan = cMBrec.focusTan;
419 pksrec.focusRot = cMBrec.focusRot;
[1325]420
[1452]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;
[1325]426
[1452]427 pksrec.refBeam = cMBrec.refBeam;
428 pksrec.beamNo = cMBrec.beamNo;
[1325]429
[1452]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;
[1427]439
[1452]440 pksrec.tsys.resize(nPol);
441 pksrec.sigma.resize(nPol);
442 pksrec.calFctr.resize(nPol);
[1325]443 for (uInt ipol = 0; ipol < nPol; ipol++) {
[1452]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];
[1325]448 }
449
[1452]450 if (cMBrec.haveBase) {
451 pksrec.baseLin.resize(2,nPol);
[1635]452 pksrec.baseSub.resize(24,nPol);
[1325]453
454 for (uInt ipol = 0; ipol < nPol; ipol++) {
[1452]455 pksrec.baseLin(0,ipol) = cMBrec.baseLin[0][ipol][0];
456 pksrec.baseLin(1,ipol) = cMBrec.baseLin[0][ipol][1];
[1325]457
[1635]458 for (uInt j = 0; j < 24; j++) {
[1452]459 pksrec.baseSub(j,ipol) = cMBrec.baseSub[0][ipol][j];
[1325]460 }
461 }
462
463 } else {
[1452]464 pksrec.baseLin.resize(0,0);
465 pksrec.baseSub.resize(0,0);
[1325]466 }
467
[1452]468 if (cGetSpectra && cMBrec.haveSpectra) {
469 pksrec.spectra.resize(nChan,nPol);
470 pksrec.spectra.takeStorage(IPosition(2,nChan,nPol), cMBrec.spectra[0],
[1427]471 SHARE);
[1325]472
[1452]473 pksrec.flagged.resize(nChan,nPol);
474 pksrec.flagged.takeStorage(IPosition(2,nChan,nPol), cMBrec.flagged[0],
[1427]475 SHARE);
[1325]476
477 } else {
[1452]478 pksrec.spectra.resize(0,0);
479 pksrec.flagged.resize(0,0);
[1325]480 }
481
482 if (cGetXPol) {
[1452]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],
[1427]487 SHARE);
[1325]488 }
489
490 return 0;
491}
492
493//------------------------------------------------------- PKSFITSreader::close
494
495// Close the FITS file.
496
497void PKSFITSreader::close()
498{
499 cReader->close();
[1452]500 logMsg(cReader->getMsg());
501 cReader->clearMsg();
[1325]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.