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

Last change on this file since 1510 was 1452, checked in by Malte Marquarding, 16 years ago

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