source: branches/plotter2/external-alma/atnf/PKSIO/PKSFITSreader.cc@ 3037

Last change on this file since 3037 was 1868, checked in by Takeshi Nakazato, 14 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): atnf

Description: Describe your changes here...

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