source: branches/mergetest/external/atnf/PKSIO/PKSFITSreader.cc@ 1775

Last change on this file since 1775 was 1757, checked in by Kana Sugimoto, 14 years ago

New Development: Yes

JIRA Issue: Yes (CAS-2211)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: ASAP 3.0.0 interface changes

Test Programs:

Put in Release Notes: Yes

Module(s): all the CASA sd tools and tasks are affected.

Description: Merged ATNF-ASAP 3.0.0 developments to CASA (alma) branch.

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