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

Last change on this file since 1400 was 1399, checked in by Malte Marquarding, 17 years ago

Mark C added brightness unit to getHeader()
File size: 14.8 KB
Line 
1//# PKSFITSreader.cc: Class to read Parkes multibeam data from a FITS file.
2//#---------------------------------------------------------------------------
3//# Copyright (C) 2000-2007
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.13 2007/11/12 03:37:56 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 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
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(
324 Int &scanNo,
325 Int &cycleNo,
326 Double &mjd,
327 Double &interval,
328 String &fieldName,
329 String &srcName,
330 Vector<Double> &srcDir,
331 Vector<Double> &srcPM,
332 Double &srcVel,
333 String &obsType,
334 Int &IFno,
335 Double &refFreq,
336 Double &bandwidth,
337 Double &freqInc,
338 Double &restFreq,
339 Vector<Float> &tcal,
340 String &tcalTime,
341 Float &azimuth,
342 Float &elevation,
343 Float &parAngle,
344 Float &focusAxi,
345 Float &focusTan,
346 Float &focusRot,
347 Float &temperature,
348 Float &pressure,
349 Float &humidity,
350 Float &windSpeed,
351 Float &windAz,
352 Int &refBeam,
353 Int &beamNo,
354 Vector<Double> &direction,
355 Vector<Double> &scanRate,
356 Vector<Float> &tsys,
357 Vector<Float> &sigma,
358 Vector<Float> &calFctr,
359 Matrix<Float> &baseLin,
360 Matrix<Float> &baseSub,
361 Matrix<Float> &spectra,
362 Matrix<uChar> &flagged,
363 Complex &xCalFctr,
364 Vector<Complex> &xPol)
365{
366 Int status;
367
368 if ((status = cReader->read(cMBrec))) {
369 if (status != -1) {
370 status = 1;
371 }
372
373 return status;
374 }
375
376
377 uInt nChan = cMBrec.nChan[0];
378 uInt nPol = cMBrec.nPol[0];
379
380 scanNo = cMBrec.scanNo;
381 cycleNo = cMBrec.cycleNo;
382
383 // Extract MJD.
384 Int day, month, year;
385 sscanf(cMBrec.datobs, "%4d-%2d-%2d", &year, &month, &day);
386 mjd = MVTime(year, month, Double(day)).day() + cMBrec.utc/86400.0;
387
388 interval = cMBrec.exposure;
389
390 fieldName = trim(cMBrec.srcName);
391 srcName = fieldName;
392 srcDir(0) = cMBrec.srcRA;
393 srcDir(1) = cMBrec.srcDec;
394 srcPM(0) = 0.0;
395 srcPM(1) = 0.0;
396 srcVel = 0.0;
397 obsType = trim(cMBrec.obsType);
398
399 IFno = cMBrec.IFno[0];
400 Double chanWidth = fabs(cMBrec.fqDelt[0]);
401 refFreq = cMBrec.fqRefVal[0];
402 bandwidth = chanWidth * nChan;
403 freqInc = cMBrec.fqDelt[0];
404 restFreq = cMBrec.restFreq;
405
406 tcal.resize(nPol);
407 for (uInt ipol = 0; ipol < nPol; ipol++) {
408 tcal(ipol) = cMBrec.tcal[0][ipol];
409 }
410 tcalTime = trim(cMBrec.tcalTime);
411 azimuth = cMBrec.azimuth;
412 elevation = cMBrec.elevation;
413 parAngle = cMBrec.parAngle;
414 focusAxi = cMBrec.focusAxi;
415 focusTan = cMBrec.focusTan;
416 focusRot = cMBrec.focusRot;
417
418 temperature = cMBrec.temp;
419 pressure = cMBrec.pressure;
420 humidity = cMBrec.humidity;
421 windSpeed = cMBrec.windSpeed;
422 windAz = cMBrec.windAz;
423
424 refBeam = cMBrec.refBeam;
425 beamNo = cMBrec.beamNo;
426
427 direction(0) = cMBrec.ra;
428 direction(1) = cMBrec.dec;
429 scanRate(0) = cMBrec.raRate;
430 scanRate(1) = cMBrec.decRate;
431
432 tsys.resize(nPol);
433 sigma.resize(nPol);
434 calFctr.resize(nPol);
435 for (uInt ipol = 0; ipol < nPol; ipol++) {
436 tsys(ipol) = cMBrec.tsys[0][ipol];
437 sigma(ipol) = tsys(ipol) / 0.81 / sqrt(interval * chanWidth);
438 calFctr(ipol) = cMBrec.calfctr[0][ipol];
439 }
440
441 if (cMBrec.haveBase) {
442 baseLin.resize(2,nPol);
443 baseSub.resize(9,nPol);
444
445 for (uInt ipol = 0; ipol < nPol; ipol++) {
446 baseLin(0,ipol) = cMBrec.baseLin[0][ipol][0];
447 baseLin(1,ipol) = cMBrec.baseLin[0][ipol][1];
448
449 for (uInt j = 0; j < 9; j++) {
450 baseSub(j,ipol) = cMBrec.baseSub[0][ipol][j];
451 }
452 }
453
454 } else {
455 baseLin.resize(0,0);
456 baseSub.resize(0,0);
457 }
458
459 if (cGetSpectra && cMBrec.haveSpectra) {
460 spectra.resize(nChan,nPol);
461 spectra.takeStorage(IPosition(2,nChan,nPol), cMBrec.spectra[0], SHARE);
462
463 flagged.resize(nChan,nPol);
464 flagged.takeStorage(IPosition(2,nChan,nPol), cMBrec.flagged[0], SHARE);
465
466 } else {
467 spectra.resize(0,0);
468 flagged.resize(0,0);
469 }
470
471 if (cGetXPol) {
472 xCalFctr = Complex(cMBrec.xcalfctr[0][0], cMBrec.xcalfctr[0][1]);
473 xPol.resize(nChan);
474 xPol.takeStorage(IPosition(1,nChan), (Complex *)cMBrec.xpol[0], SHARE);
475 }
476
477 return 0;
478}
479
480//-------------------------------------------------------- PKSFITSreader::read
481
482// Read the next data record, just the basics.
483
484Int PKSFITSreader::read(
485 Int &IFno,
486 Vector<Float> &tsys,
487 Vector<Float> &calFctr,
488 Matrix<Float> &baseLin,
489 Matrix<Float> &baseSub,
490 Matrix<Float> &spectra,
491 Matrix<uChar> &flagged)
492{
493 Int status;
494
495 if ((status = cReader->read(cMBrec))) {
496 if (status != -1) {
497 status = 1;
498 }
499
500 return status;
501 }
502
503 IFno = cMBrec.IFno[0];
504
505 uInt nChan = cMBrec.nChan[0];
506 uInt nPol = cMBrec.nPol[0];
507
508 tsys.resize(nPol);
509 calFctr.resize(nPol);
510 for (uInt ipol = 0; ipol < nPol; ipol++) {
511 tsys(ipol) = cMBrec.tsys[0][ipol];
512 calFctr(ipol) = cMBrec.calfctr[0][ipol];
513 }
514
515 if (cMBrec.haveBase) {
516 baseLin.resize(2,nPol);
517 baseSub.resize(9,nPol);
518
519 for (uInt ipol = 0; ipol < nPol; ipol++) {
520 baseLin(0,ipol) = cMBrec.baseLin[0][ipol][0];
521 baseLin(1,ipol) = cMBrec.baseLin[0][ipol][1];
522
523 for (uInt j = 0; j < 9; j++) {
524 baseSub(j,ipol) = cMBrec.baseSub[0][ipol][j];
525 }
526 }
527
528 } else {
529 baseLin.resize(0,0);
530 baseSub.resize(0,0);
531 }
532
533 if (cGetSpectra && cMBrec.haveSpectra) {
534 spectra.resize(nChan,nPol);
535 spectra.takeStorage(IPosition(2,nChan,nPol), cMBrec.spectra[0], SHARE);
536
537 flagged.resize(nChan,nPol);
538 flagged.takeStorage(IPosition(2,nChan,nPol), cMBrec.flagged[0], SHARE);
539
540 } else {
541 spectra.resize(0,0);
542 flagged.resize(0,0);
543 }
544
545 return 0;
546}
547
548//------------------------------------------------------- PKSFITSreader::close
549
550// Close the FITS file.
551
552void PKSFITSreader::close()
553{
554 cReader->close();
555}
556
557//-------------------------------------------------------- PKSFITSreader::trim
558
559// Trim trailing blanks from a null-terminated character string.
560
561char* PKSFITSreader::trim(char *string)
562{
563 int j = 0, k = 0;
564 while (string[j] != '\0') {
565 if (string[j++] != ' ') {
566 k = j;
567 }
568 }
569
570 string[k] = '\0';
571
572 return string;
573}
Note: See TracBrowser for help on using the repository browser.