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

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

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