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

Last change on this file since 1386 was 1325, checked in by mar637, 18 years ago

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