source: branches/alma/external/atnf/PKSIO/PKSFITSreader.cc@ 1539

Last change on this file since 1539 was 1453, checked in by TakTsutsumi, 16 years ago

New Development: No

JIRA Issue: No

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: many

Test Programs: sd.scantable(), sd.scantable.save()

Put in Release Notes: N/A

Description: copied from current casapy code tree


File size: 14.8 KB
RevLine 
[1325]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//#
[1393]28//# $Id$
[1325]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,
[1393]144 Double &bandwidth,
145 String &fluxunit)
[1325]146{
147 char datobs[32], dopplerFrame_[32], observer_[32], obsType_[32],
148 project_[32], radecsys[32], telescope[32];
149 float equinox_;
150 double antPos[3], utc;
151
152 if (cReader->getHeader(observer_, project_, telescope, antPos, obsType_,
153 equinox_, radecsys, dopplerFrame_, datobs, utc,
154 refFreq, bandwidth)) {
155 return 1;
156 }
157
[1393]158 fluxunit = "";
[1325]159 observer = trim(observer_);
160 project = trim(project_);
161 antName = trim(telescope);
162 antPosition.resize(3);
163 antPosition(0) = antPos[0];
164 antPosition(1) = antPos[1];
165 antPosition(2) = antPos[2];
166 obsType = trim(obsType_);
167 equinox = equinox_;
168 dopplerFrame = trim(dopplerFrame_);
169
170 // Get UTC in MJD form.
171 Int day, month, year;
172 sscanf(datobs, "%4d-%2d-%2d", &year, &month, &day);
173 MVTime date(year, month, Double(day));
174 mjd = date.day() + utc/86400.0;
175
176 return 0;
177}
178
179//------------------------------------------------- PKSFITSreader::getFreqInfo
180
181// Get frequency parameters for each IF.
182
183Int PKSFITSreader::getFreqInfo(
184 Vector<Double> &startFreq,
185 Vector<Double> &endFreq)
186{
187 int nIF;
188 double *startfreq, *endfreq;
189
190 Int status;
191 if (!(status = cReader->getFreqInfo(nIF, startfreq, endfreq))) {
192 startFreq.takeStorage(IPosition(1,nIF), startfreq, TAKE_OVER);
193 endFreq.takeStorage(IPosition(1,nIF), endfreq, TAKE_OVER);
194 }
195
196 return status;
197}
198
199//------------------------------------------------------ PKSFITSreader::select
200
201// Set data selection by beam number and channel.
202
203uInt PKSFITSreader::select(
204 const Vector<Bool> beamSel,
205 const Vector<Bool> IFsel,
206 const Vector<Int> startChan,
207 const Vector<Int> endChan,
208 const Vector<Int> refChan,
209 const Bool getSpectra,
210 const Bool getXPol,
211 const Bool getFeedPos)
212{
213 // Apply beam selection.
214 uInt nBeamSel = beamSel.nelements();
215 for (uInt ibeam = 0; ibeam < cNBeam; ibeam++) {
216 // This modifies FITSreader::cBeams[].
217 if (ibeam < nBeamSel) {
218 cBeams[ibeam] = cBeams[ibeam] && beamSel(ibeam);
219 } else {
220 cBeams[ibeam] = 0;
221 }
222 }
223
224 // Apply IF selection.
225 int *end = new int[cNIF];
226 int *ref = new int[cNIF];
227 int *start = new int[cNIF];
228
229 for (uInt iIF = 0; iIF < cNIF; iIF++) {
230 // This modifies FITSreader::cIFs[].
231 if (iIF < IFsel.nelements()) {
232 cIFs[iIF] = cIFs[iIF] && IFsel(iIF);
233 } else {
234 cIFs[iIF] = 0;
235 }
236
237 if (cIFs[iIF]) {
238 if (iIF < startChan.nelements()) {
239 start[iIF] = startChan(iIF);
240
241 if (start[iIF] <= 0) {
242 start[iIF] += cNChan(iIF);
243 } else if (start[iIF] > Int(cNChan(iIF))) {
244 start[iIF] = cNChan(iIF);
245 }
246
247 } else {
248 start[iIF] = 1;
249 }
250
251 if (iIF < endChan.nelements()) {
252 end[iIF] = endChan(iIF);
253
254 if (end[iIF] <= 0) {
255 end[iIF] += cNChan(iIF);
256 } else if (end[iIF] > Int(cNChan(iIF))) {
257 end[iIF] = cNChan(iIF);
258 }
259
260 } else {
261 end[iIF] = cNChan(iIF);
262 }
263
264 if (iIF < refChan.nelements()) {
265 ref[iIF] = refChan(iIF);
266 } else {
267 if (start[iIF] <= end[iIF]) {
268 ref[iIF] = start[iIF] + (end[iIF] - start[iIF] + 1)/2;
269 } else {
270 ref[iIF] = start[iIF] - (start[iIF] - end[iIF] + 1)/2;
271 }
272 }
273 }
274 }
275
276 cGetSpectra = getSpectra;
277 cGetXPol = getXPol;
278 cGetFeedPos = getFeedPos;
279
280 uInt maxNChan = cReader->select(start, end, ref, cGetSpectra, cGetXPol,
281 cGetFeedPos);
282
283 delete [] end;
284 delete [] ref;
285 delete [] start;
286
287 return maxNChan;
288}
289
290//--------------------------------------------------- PKSFITSreader::findRange
291
292// Read the TIME column.
293
294Int PKSFITSreader::findRange(
295 Int &nRow,
296 Int &nSel,
297 Vector<Double> &timeSpan,
298 Matrix<Double> &positions)
299{
300 char dateSpan[2][32];
301 double utcSpan[2];
302 double* posns;
303
304 Int status;
305 if (!(status = cReader->findRange(nRow, nSel, dateSpan, utcSpan, posns))) {
306 timeSpan.resize(2);
307
308 Int day, month, year;
309 sscanf(dateSpan[0], "%4d-%2d-%2d", &year, &month, &day);
310 timeSpan(0) = MVTime(year, month, Double(day)).second() + utcSpan[0];
311 sscanf(dateSpan[1], "%4d-%2d-%2d", &year, &month, &day);
312 timeSpan(1) = MVTime(year, month, Double(day)).second() + utcSpan[1];
313
314 positions.takeStorage(IPosition(2,2,nSel), posns, TAKE_OVER);
315 }
316
317 return status;
318}
319
320//-------------------------------------------------------- PKSFITSreader::read
321
322// Read the next data record.
323
324Int PKSFITSreader::read(
325 Int &scanNo,
326 Int &cycleNo,
327 Double &mjd,
328 Double &interval,
329 String &fieldName,
330 String &srcName,
331 Vector<Double> &srcDir,
332 Vector<Double> &srcPM,
333 Double &srcVel,
334 String &obsType,
335 Int &IFno,
336 Double &refFreq,
337 Double &bandwidth,
338 Double &freqInc,
[1453]339 Vector<Double> &restFreq,
[1325]340 Vector<Float> &tcal,
341 String &tcalTime,
342 Float &azimuth,
343 Float &elevation,
344 Float &parAngle,
345 Float &focusAxi,
346 Float &focusTan,
347 Float &focusRot,
348 Float &temperature,
349 Float &pressure,
350 Float &humidity,
351 Float &windSpeed,
352 Float &windAz,
353 Int &refBeam,
354 Int &beamNo,
355 Vector<Double> &direction,
356 Vector<Double> &scanRate,
357 Vector<Float> &tsys,
358 Vector<Float> &sigma,
359 Vector<Float> &calFctr,
360 Matrix<Float> &baseLin,
361 Matrix<Float> &baseSub,
362 Matrix<Float> &spectra,
363 Matrix<uChar> &flagged,
364 Complex &xCalFctr,
365 Vector<Complex> &xPol)
366{
367 Int status;
368
369 if ((status = cReader->read(cMBrec))) {
370 if (status != -1) {
371 status = 1;
372 }
373
374 return status;
375 }
376
377
378 uInt nChan = cMBrec.nChan[0];
379 uInt nPol = cMBrec.nPol[0];
380
381 scanNo = cMBrec.scanNo;
382 cycleNo = cMBrec.cycleNo;
383
384 // Extract MJD.
385 Int day, month, year;
386 sscanf(cMBrec.datobs, "%4d-%2d-%2d", &year, &month, &day);
387 mjd = MVTime(year, month, Double(day)).day() + cMBrec.utc/86400.0;
388
389 interval = cMBrec.exposure;
390
391 fieldName = trim(cMBrec.srcName);
392 srcName = fieldName;
393 srcDir(0) = cMBrec.srcRA;
394 srcDir(1) = cMBrec.srcDec;
395 srcPM(0) = 0.0;
396 srcPM(1) = 0.0;
397 srcVel = 0.0;
398 obsType = trim(cMBrec.obsType);
399
400 IFno = cMBrec.IFno[0];
401 Double chanWidth = fabs(cMBrec.fqDelt[0]);
402 refFreq = cMBrec.fqRefVal[0];
403 bandwidth = chanWidth * nChan;
404 freqInc = cMBrec.fqDelt[0];
[1453]405 //restFreq = cMBrec.restFreq;
406 restFreq.resize(1);
407 restFreq(0) = cMBrec.restFreq;
[1325]408
409 tcal.resize(nPol);
410 for (uInt ipol = 0; ipol < nPol; ipol++) {
411 tcal(ipol) = cMBrec.tcal[0][ipol];
412 }
413 tcalTime = trim(cMBrec.tcalTime);
414 azimuth = cMBrec.azimuth;
415 elevation = cMBrec.elevation;
416 parAngle = cMBrec.parAngle;
417 focusAxi = cMBrec.focusAxi;
418 focusTan = cMBrec.focusTan;
419 focusRot = cMBrec.focusRot;
420
421 temperature = cMBrec.temp;
422 pressure = cMBrec.pressure;
423 humidity = cMBrec.humidity;
424 windSpeed = cMBrec.windSpeed;
425 windAz = cMBrec.windAz;
426
427 refBeam = cMBrec.refBeam;
428 beamNo = cMBrec.beamNo;
429
430 direction(0) = cMBrec.ra;
431 direction(1) = cMBrec.dec;
432 scanRate(0) = cMBrec.raRate;
433 scanRate(1) = cMBrec.decRate;
434
435 tsys.resize(nPol);
436 sigma.resize(nPol);
437 calFctr.resize(nPol);
438 for (uInt ipol = 0; ipol < nPol; ipol++) {
439 tsys(ipol) = cMBrec.tsys[0][ipol];
440 sigma(ipol) = tsys(ipol) / 0.81 / sqrt(interval * chanWidth);
441 calFctr(ipol) = cMBrec.calfctr[0][ipol];
442 }
443
444 if (cMBrec.haveBase) {
445 baseLin.resize(2,nPol);
446 baseSub.resize(9,nPol);
447
448 for (uInt ipol = 0; ipol < nPol; ipol++) {
449 baseLin(0,ipol) = cMBrec.baseLin[0][ipol][0];
450 baseLin(1,ipol) = cMBrec.baseLin[0][ipol][1];
451
452 for (uInt j = 0; j < 9; j++) {
453 baseSub(j,ipol) = cMBrec.baseSub[0][ipol][j];
454 }
455 }
456
457 } else {
458 baseLin.resize(0,0);
459 baseSub.resize(0,0);
460 }
461
462 if (cGetSpectra && cMBrec.haveSpectra) {
463 spectra.resize(nChan,nPol);
464 spectra.takeStorage(IPosition(2,nChan,nPol), cMBrec.spectra[0], SHARE);
465
466 flagged.resize(nChan,nPol);
467 flagged.takeStorage(IPosition(2,nChan,nPol), cMBrec.flagged[0], SHARE);
468
469 } else {
470 spectra.resize(0,0);
471 flagged.resize(0,0);
472 }
473
474 if (cGetXPol) {
475 xCalFctr = Complex(cMBrec.xcalfctr[0][0], cMBrec.xcalfctr[0][1]);
476 xPol.resize(nChan);
477 xPol.takeStorage(IPosition(1,nChan), (Complex *)cMBrec.xpol[0], SHARE);
478 }
479
480 return 0;
481}
482
483//-------------------------------------------------------- PKSFITSreader::read
484
485// Read the next data record, just the basics.
486
487Int PKSFITSreader::read(
488 Int &IFno,
489 Vector<Float> &tsys,
490 Vector<Float> &calFctr,
491 Matrix<Float> &baseLin,
492 Matrix<Float> &baseSub,
493 Matrix<Float> &spectra,
494 Matrix<uChar> &flagged)
495{
496 Int status;
497
498 if ((status = cReader->read(cMBrec))) {
499 if (status != -1) {
500 status = 1;
501 }
502
503 return status;
504 }
505
506 IFno = cMBrec.IFno[0];
507
508 uInt nChan = cMBrec.nChan[0];
509 uInt nPol = cMBrec.nPol[0];
510
511 tsys.resize(nPol);
512 calFctr.resize(nPol);
513 for (uInt ipol = 0; ipol < nPol; ipol++) {
514 tsys(ipol) = cMBrec.tsys[0][ipol];
515 calFctr(ipol) = cMBrec.calfctr[0][ipol];
516 }
517
518 if (cMBrec.haveBase) {
519 baseLin.resize(2,nPol);
520 baseSub.resize(9,nPol);
521
522 for (uInt ipol = 0; ipol < nPol; ipol++) {
523 baseLin(0,ipol) = cMBrec.baseLin[0][ipol][0];
524 baseLin(1,ipol) = cMBrec.baseLin[0][ipol][1];
525
526 for (uInt j = 0; j < 9; j++) {
527 baseSub(j,ipol) = cMBrec.baseSub[0][ipol][j];
528 }
529 }
530
531 } else {
532 baseLin.resize(0,0);
533 baseSub.resize(0,0);
534 }
535
536 if (cGetSpectra && cMBrec.haveSpectra) {
537 spectra.resize(nChan,nPol);
538 spectra.takeStorage(IPosition(2,nChan,nPol), cMBrec.spectra[0], SHARE);
539
540 flagged.resize(nChan,nPol);
541 flagged.takeStorage(IPosition(2,nChan,nPol), cMBrec.flagged[0], SHARE);
542
543 } else {
544 spectra.resize(0,0);
545 flagged.resize(0,0);
546 }
547
548 return 0;
549}
550
551//------------------------------------------------------- PKSFITSreader::close
552
553// Close the FITS file.
554
555void PKSFITSreader::close()
556{
557 cReader->close();
558}
559
560//-------------------------------------------------------- PKSFITSreader::trim
561
562// Trim trailing blanks from a null-terminated character string.
563
564char* PKSFITSreader::trim(char *string)
565{
566 int j = 0, k = 0;
567 while (string[j] != '\0') {
568 if (string[j++] != ' ') {
569 k = j;
570 }
571 }
572
573 string[k] = '\0';
574
575 return string;
576}
Note: See TracBrowser for help on using the repository browser.