source: trunk/external-alma/atnf/PKSIO/GBTFITSreader.cc@ 2011

Last change on this file since 2011 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: 84.1 KB
Line 
1//#---------------------------------------------------------------------------
2//# GBTFITSreader.cc: interface class for GBT SDFITS input using CFITSIO.
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: SDFITSreader.cc,v 19.45 2009-09-30 07:23:48 cal103 Exp $
32//#---------------------------------------------------------------------------
33//# The GBTFITSreader class reads single dish FITS files from Green Bank
34//# telescopes. This class is actually a specific version of SDFITSreader.
35//#
36//# Original: 2000/08/09, Mark Calabretta, ATNF
37//# Update: 2010/07/29, Takeshi Nakazato, NAOJ
38//#---------------------------------------------------------------------------
39
40#include <atnf/pks/pks_maths.h>
41#include <atnf/PKSIO/MBrecord.h>
42#include <atnf/PKSIO/GBTFITSreader.h>
43#include <atnf/PKSIO/SrcType.h>
44
45#include <measures/Measures/MeasConvert.h>
46#include <measures/Measures/MCFrequency.h>
47#include <measures/Measures/MFrequency.h>
48#include <measures/Measures/MPosition.h>
49#include <measures/Measures/MEpoch.h>
50#include <measures/Measures/MDirection.h>
51
52#include <casa/Logging/LogIO.h>
53#include <casa/Quanta/MVTime.h>
54#include <casa/math.h>
55#include <casa/stdio.h>
56
57#include <algorithm>
58#include <strings.h>
59#include <cstring>
60
61class FITSparm
62{
63 public:
64 char *name; // Keyword or column name.
65 int type; // Expected keyvalue or column data type.
66 int colnum; // Column number; 0 for keyword; -1 absent.
67 int coltype; // Column data type, as found.
68 long nelem; // Column data repeat count; < 0 for vardim.
69 int tdimcol; // TDIM column number; 0 for keyword; -1 absent.
70 char units[32]; // Units from TUNITn keyword.
71};
72
73// Numerical constants.
74const double PI = 3.141592653589793238462643;
75
76// Factor to convert radians to degrees.
77const double D2R = PI / 180.0;
78
79// Class name
80const string className = "GBTFITSreader" ;
81
82//---------------------------------------------------- GBTFITSreader::(statics)
83
84int GBTFITSreader::sInit = 1;
85int GBTFITSreader::sReset = 0;
86int (*GBTFITSreader::sALFAcalNon)[2] = (int (*)[2])(new float[16]);
87int (*GBTFITSreader::sALFAcalNoff)[2] = (int (*)[2])(new float[16]);
88float (*GBTFITSreader::sALFAcalOn)[2] = (float (*)[2])(new float[16]);
89float (*GBTFITSreader::sALFAcalOff)[2] = (float (*)[2])(new float[16]);
90float (*GBTFITSreader::sALFAcal)[2] = (float (*)[2])(new float[16]);
91
92//------------------------------------------------- GBTFITSreader::GBTFITSreader
93
94GBTFITSreader::GBTFITSreader()
95{
96 // Default constructor.
97 cSDptr = 0x0;
98
99 // Allocate space for data descriptors.
100 cData = new FITSparm[NDATA];
101
102 for (int iData = 0; iData < NDATA; iData++) {
103 cData[iData].colnum = -1;
104 }
105
106 // Initialize pointers.
107 cBeams = 0x0;
108 cIFs = 0x0;
109 cStartChan = 0x0;
110 cEndChan = 0x0;
111 cRefChan = 0x0;
112 cPols = 0x0;
113
114 // other initialization
115 numHDU = 0 ;
116
117}
118
119//------------------------------------------------ GBTFITSreader::~GBTFITSreader
120
121GBTFITSreader::~GBTFITSreader()
122{
123 close();
124
125 delete [] cData;
126}
127
128//--------------------------------------------------------- GBTFITSreader::open
129
130// Open an SDFITS file for reading.
131
132int GBTFITSreader::open(
133 char* sdName,
134 int &nBeam,
135 int* &beams,
136 int &nIF,
137 int* &IFs,
138 int* &nChan,
139 int* &nPol,
140 int* &haveXPol,
141 int &haveBase,
142 int &haveSpectra,
143 int &extraSysCal)
144{
145 const string methodName = "open()" ;
146
147 if (cSDptr) {
148 close();
149 }
150
151 // Open the SDFITS file.
152 cStatus = 0;
153 if (fits_open_file(&cSDptr, sdName, READONLY, &cStatus)) {
154 sprintf(cMsg, "ERROR: Failed to open SDFITS file\n %s", sdName);
155 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, cMsg);
156 return 1;
157 }
158
159 // Get number of HDUs in the file
160 if ( fits_get_num_hdus( cSDptr, &numHDU, &cStatus ) ) {
161 sprintf(cMsg, "ERROR: Failed to get number of HDUs in the file\n %s", sdName );
162 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, cMsg) ;
163 return 1 ;
164 }
165 // ignore Primary HDU
166 numHDU = numHDU - 1 ;
167 sprintf(cMsg, "number of HDUs = %d", numHDU) ;
168 log( LogOrigin( className, methodName, WHERE ), LogIO::DEBUGGING, cMsg ) ;
169
170 // Move to the SDFITS extension.
171 cALFA = cALFA_BD = cALFA_CIMA = 0;
172 char *extName = "SINGLE DISH" ;
173 if (fits_movnam_hdu(cSDptr, BINARY_TBL, "SINGLE DISH", 0, &cStatus)) {
174 // No SDFITS table, look for BDFITS or CIMAFITS.
175 cStatus = 0;
176 if (fits_movnam_hdu(cSDptr, BINARY_TBL, "BDFITS", 0, &cStatus) == 0) {
177 cALFA_BD = 1;
178 extName = "BDFITS" ;
179
180 } else {
181 cStatus = 0;
182 if (fits_movnam_hdu(cSDptr, BINARY_TBL, "CIMAFITS", 0, &cStatus) == 0) {
183 cALFA_CIMA = 1;
184 extName = "CIMAFITS" ;
185
186 // Check for later versions of CIMAFITS.
187 float version;
188 readParm("VERSION", TFLOAT, &version);
189 if (version >= 2.0f) cALFA_CIMA = int(version);
190
191 } else {
192 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Failed to locate SDFITS binary table.");
193 close();
194 return 1;
195 }
196 }
197
198 // Arecibo ALFA data of some kind.
199 cALFA = 1;
200 if (sInit) {
201 for (int iBeam = 0; iBeam < 8; iBeam++) {
202 for (int iPol = 0; iPol < 2; iPol++) {
203 sALFAcalOn[iBeam][iPol] = 0.0f;
204 sALFAcalOff[iBeam][iPol] = 0.0f;
205
206 // Nominal factor to calibrate spectra in Jy.
207 sALFAcal[iBeam][iPol] = 3.0f;
208 }
209 }
210
211 sInit = 0;
212 }
213 }
214
215 // Find the number of rows.
216 cNRowT = new long[numHDU] ;
217 cNRowCum = new long[numHDU] ;
218 hduId = new int[numHDU] ;
219 cNRow = 0 ;
220 fits_movabs_hdu( cSDptr, 1, NULL, &cStatus ) ;
221 for ( int i = 0 ; i < numHDU ; i++ ) {
222 if ( i == 0 ) {
223 fits_movnam_hdu( cSDptr, BINARY_TBL, extName, 0, &cStatus ) ;
224 fits_get_hdu_num( cSDptr, &(hduId[i]) ) ;
225 }
226 else {
227 if ( fits_movrel_hdu( cSDptr, 1, NULL, &cStatus ) ) {
228 log( LogOrigin( className, methodName, WHERE), LogIO::SEVERE, "Failed to move" ) ;
229 }
230 fits_get_hdu_num( cSDptr, &(hduId[i]) ) ;
231 }
232 sprintf(cMsg, "HDU ID for %d: %d", i, hduId[i]) ;
233 log( LogOrigin( className, methodName, WHERE ), LogIO::DEBUGGING, cMsg ) ;
234 fits_get_num_rows( cSDptr, &cNRowT[i], &cStatus ) ;
235 if (!cNRowT[i]) {
236 ostringstream oss ;
237 sprintf(cMsg, "Table %d contains no entries.", i ) ;
238 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, cMsg);
239 close();
240 return 1;
241 }
242 cNRow += cNRowT[i] ;
243 cNRowCum[i] = cNRow ;
244 sprintf(cMsg, "number of rows for HDU %d = %ld", i, cNRowT[i]) ;
245 log( LogOrigin( className, methodName, WHERE ), LogIO::DEBUGGING, cMsg ) ;
246 sprintf(cMsg, "cumulative number of rows for HDU %d = %ld", i, cNRowCum[i]) ;
247 log( LogOrigin( className, methodName, WHERE ), LogIO::DEBUGGING, cMsg ) ;
248 }
249 sprintf(cMsg, "number of rows = %ld", cNRow) ;
250 log( LogOrigin( className, methodName, WHERE ), LogIO::DEBUGGING, cMsg ) ;
251
252 // Move back to first binary table
253 fits_movabs_hdu( cSDptr, hduId[0], NULL, &cStatus ) ;
254
255 //////
256 // For devlopment
257// cNRow = 0 ;
258// for ( int ihdu = 0 ; ihdu < numHDU ; ihdu++ ) {
259// //cNRowT[ihdu] = 10 ;
260// //cNRowT[ihdu] = 100 ;
261// cNRowT[ihdu] = 500 ;
262// cNRow += cNRowT[ihdu] ;
263// cNRowCum[ihdu] = cNRow ;
264// sprintf(cMsg, "number of rows for HDU %d = %ld", ihdu, cNRowT[ihdu]) ;
265// log( LogOrigin( className, methodName, WHERE ), LogIO::DEBUGGING, cMsg ) ;
266// sprintf(cMsg, "cumulative number of rows for HDU %d = %ld", ihdu, cNRowCum[ihdu]) ;
267// log( LogOrigin( className, methodName, WHERE ), LogIO::DEBUGGING, cMsg ) ;
268// }
269// sprintf(cMsg, "number of rows = %ld", cNRow) ;
270// log( LogOrigin( className, methodName, WHERE ), LogIO::DEBUGGING, cMsg ) ;
271 // Original cNRow is cNRowT[0]
272 // cNRow = cNRowT[0] ;
273 // fits_get_num_rows(cSDptr, &cNRow, &cStatus);
274 /////
275
276 if (!cNRow) {
277 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Table contains no entries.");
278 close();
279 return 1;
280 }
281
282 // GBT data.
283 char telescope[32];
284 readParm("TELESCOP", TSTRING, telescope); // Core.
285 cGBT = strncmp(telescope, "GBT", 3) == 0 ||
286 strncmp(telescope, "NRAO_GBT", 8) == 0;
287
288
289 // Check that the DATA array column is present.
290 findData(DATA, "DATA", TFLOAT);
291 haveSpectra = cHaveSpectra = cData[DATA].colnum > 0;
292
293 cNAxisTime = 0;
294 if (cHaveSpectra) {
295 // Find the number of data axes (must be the same for each IF).
296 cNAxes = 5;
297 if (readDim(DATA, 1, &cNAxes, cNAxis)) {
298 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
299 close();
300 return 1;
301 }
302
303 if (cALFA_BD) {
304 // ALFA BDFITS: variable length arrays don't actually vary and there is
305 // no TDIM (or MAXISn) card; use the LAGS_IN value.
306 cNAxes = 5;
307 readParm("LAGS_IN", TLONG, cNAxis);
308 cNAxis[1] = 1;
309 cNAxis[2] = 1;
310 cNAxis[3] = 1;
311 cNAxis[4] = 1;
312 cData[DATA].nelem = cNAxis[0];
313 }
314
315 if (cNAxes < 4) {
316 // Need at least four axes (for now).
317 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "DATA array contains fewer than four axes.");
318 close();
319 return 1;
320 } else if (cNAxes > 5) {
321 // We support up to five axes.
322 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "DATA array contains more than five axes.");
323 close();
324 return 1;
325 }
326
327 findData(FLAGGED, "FLAGGED", TBYTE);
328
329 } else {
330 // DATA column not present, check for a DATAXED keyword.
331 findData(DATAXED, "DATAXED", TSTRING);
332 if (cData[DATAXED].colnum < 0) {
333 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "DATA array column absent from binary table.");
334 close();
335 return 1;
336 }
337
338 // Determine the number of axes and their length.
339 char dataxed[32];
340 readParm("DATAXED", TSTRING, dataxed);
341
342 for (int iaxis = 0; iaxis < 5; iaxis++) cNAxis[iaxis] = 0;
343 sscanf(dataxed, "(%ld,%ld,%ld,%ld,%ld)", cNAxis, cNAxis+1, cNAxis+2,
344 cNAxis+3, cNAxis+4);
345 for (int iaxis = 4; iaxis > -1; iaxis--) {
346 if (cNAxis[iaxis] == 0) cNAxes = iaxis;
347 }
348 }
349
350 char *CTYPE[5] = {"CTYPE1", "CTYPE2", "CTYPE3", "CTYPE4", "CTYPE5"};
351 char *CRPIX[5] = {"CRPIX1", "CRPIX2", "CRPIX3", "CRPIX4", "CRPIX5"};
352 char *CRVAL[5] = {"CRVAL1", "CRVAL2", "CRVAL3", "CRVAL4", "CRVAL5"};
353 char *CDELT[5] = {"CDELT1", "CDELT2", "CDELT3", "CDELT4", "CDELT5"};
354
355 // Find required DATA array axes.
356 char ctype[5][72];
357 for (int iaxis = 0; iaxis < cNAxes; iaxis++) {
358 strcpy(ctype[iaxis], "");
359 readParm(CTYPE[iaxis], TSTRING, ctype[iaxis]); // Core.
360 }
361
362 if (cStatus) {
363 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
364 close();
365 return 1;
366 }
367
368 char *fqCRVAL = 0;
369 char *fqCDELT = 0;
370 char *fqCRPIX = 0;
371 char *raCRVAL = 0;
372 char *decCRVAL = 0;
373 char *timeCRVAL = 0;
374 char *timeCDELT = 0;
375 char *timeCRPIX = 0;
376 char *beamCRVAL = 0;
377 char *polCRVAL = 0;
378
379 cFreqAxis = -1;
380 cStokesAxis = -1;
381 cRaAxis = -1;
382 cDecAxis = -1;
383 cTimeAxis = -1;
384 cBeamAxis = -1;
385
386 for (int iaxis = 0; iaxis < cNAxes; iaxis++) {
387 if (strncmp(ctype[iaxis], "FREQ", 4) == 0) {
388 cFreqAxis = iaxis;
389 fqCRVAL = CRVAL[iaxis];
390 fqCDELT = CDELT[iaxis];
391 fqCRPIX = CRPIX[iaxis];
392
393 } else if (strncmp(ctype[iaxis], "STOKES", 6) == 0) {
394 cStokesAxis = iaxis;
395 polCRVAL = CRVAL[iaxis];
396
397 } else if (strncmp(ctype[iaxis], "RA", 2) == 0) {
398 cRaAxis = iaxis;
399 raCRVAL = CRVAL[iaxis];
400
401 } else if (strncmp(ctype[iaxis], "DEC", 3) == 0) {
402 cDecAxis = iaxis;
403 decCRVAL = CRVAL[iaxis];
404
405 } else if (strcmp(ctype[iaxis], "TIME") == 0) {
406 // TIME (UTC seconds since midnight); axis type, if present, takes
407 // precedence over keyword.
408 cTimeAxis = iaxis;
409 timeCRVAL = CRVAL[iaxis];
410
411 // Check for non-degeneracy.
412 if ((cNAxisTime = cNAxis[iaxis]) > 1) {
413 timeCDELT = CDELT[iaxis];
414 timeCRPIX = CRPIX[iaxis];
415 sprintf(cMsg, "DATA array contains a TIME axis of length %ld.",
416 cNAxisTime);
417 //logMsg(cMsg);
418 log(LogOrigin( className, methodName, WHERE ), LogIO::NORMAL, cMsg);
419 }
420
421 } else if (strcmp(ctype[iaxis], "BEAM") == 0) {
422 // BEAM can be a keyword or axis type.
423 cBeamAxis = iaxis;
424 beamCRVAL = CRVAL[iaxis];
425 }
426 }
427
428 if (cALFA_BD) {
429 // Fixed in ALFA CIMAFITS.
430 cRaAxis = 2;
431 raCRVAL = "CRVAL2A";
432
433 cDecAxis = 3;
434 decCRVAL = "CRVAL3A";
435 }
436
437
438 // Check that required axes are present.
439 if (cFreqAxis < 0 || cStokesAxis < 0 || cRaAxis < 0 || cDecAxis < 0) {
440 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Could not find required DATA array axes.");
441 close();
442 return 1;
443 }
444
445 // Set up machinery for data retrieval.
446 findData(SCAN, "SCAN", TINT); // Shared.
447 findData(CYCLE, "CYCLE", TINT); // Additional.
448 findData(DATE_OBS, "DATE-OBS", TSTRING); // Core.
449
450 if (cTimeAxis >= 0) {
451 // The DATA array has a TIME axis.
452 if (cNAxisTime > 1) {
453 // Non-degenerate.
454 findData(TimeRefVal, timeCRVAL, TDOUBLE); // Time reference value.
455 findData(TimeDelt, timeCDELT, TDOUBLE); // Time increment.
456 findData(TimeRefPix, timeCRPIX, TFLOAT); // Time reference pixel.
457 } else {
458 // Degenerate, treat its like a simple TIME keyword.
459 findData(TIME, timeCRVAL, TDOUBLE);
460 }
461
462 } else {
463 findData(TIME, "TIME", TDOUBLE); // Core.
464 }
465
466 findData(EXPOSURE, "EXPOSURE", TFLOAT); // Core.
467 findData(OBJECT, "OBJECT", TSTRING); // Core.
468 findData(OBJ_RA, "OBJ-RA", TDOUBLE); // Additional.
469 findData(OBJ_DEC, "OBJ-DEC", TDOUBLE); // Additional.
470 findData(RESTFRQ, "RESTFRQ", TDOUBLE); // Additional.
471 findData(OBSMODE, "OBSMODE", TSTRING); // Shared.
472
473 findData(BEAM, "BEAM", TSHORT); // Additional.
474 findData(IF, "IF", TSHORT); // Additional.
475 findData(FqRefVal, fqCRVAL, TDOUBLE); // Frequency reference value.
476 findData(FqDelt, fqCDELT, TDOUBLE); // Frequency increment.
477 findData(FqRefPix, fqCRPIX, TFLOAT); // Frequency reference pixel.
478 findData(RA, raCRVAL, TDOUBLE); // Right ascension.
479 findData(DEC, decCRVAL, TDOUBLE); // Declination.
480 findData(SCANRATE, "SCANRATE", TFLOAT); // Additional.
481
482 findData(TSYS, "TSYS", TFLOAT); // Core.
483 findData(CALFCTR, "CALFCTR", TFLOAT); // Additional.
484 findData(XCALFCTR, "XCALFCTR", TFLOAT); // Additional.
485 findData(BASELIN, "BASELIN", TFLOAT); // Additional.
486 findData(BASESUB, "BASESUB", TFLOAT); // Additional.
487 findData(XPOLDATA, "XPOLDATA", TFLOAT); // Additional.
488
489 findData(REFBEAM, "REFBEAM", TSHORT); // Additional.
490 findData(TCAL, "TCAL", TFLOAT); // Shared.
491 findData(TCALTIME, "TCALTIME", TSTRING); // Additional.
492 findData(AZIMUTH, "AZIMUTH", TFLOAT); // Shared.
493 findData(ELEVATIO, "ELEVATIO", TFLOAT); // Shared.
494 findData(PARANGLE, "PARANGLE", TFLOAT); // Additional.
495 findData(FOCUSAXI, "FOCUSAXI", TFLOAT); // Additional.
496 findData(FOCUSTAN, "FOCUSTAN", TFLOAT); // Additional.
497 findData(FOCUSROT, "FOCUSROT", TFLOAT); // Additional.
498 findData(TAMBIENT, "TAMBIENT", TFLOAT); // Shared.
499 findData(PRESSURE, "PRESSURE", TFLOAT); // Shared.
500 findData(HUMIDITY, "HUMIDITY", TFLOAT); // Shared.
501 findData(WINDSPEE, "WINDSPEE", TFLOAT); // Shared.
502 findData(WINDDIRE, "WINDDIRE", TFLOAT); // Shared.
503
504 findData(STOKES, polCRVAL, TINT);
505 findData(SIG, "SIG", TSTRING);
506 findData(CAL, "CAL", TSTRING);
507
508 findData(RVSYS, "RVSYS", TDOUBLE);
509 findData(VFRAME, "VFRAME", TDOUBLE);
510 findData(VELDEF, "VELDEF", TSTRING);
511
512 findData(TIMESTAMP, "TIMESTAMP", TSTRING);
513 findData(DURATION, "DURATION", TDOUBLE);
514 findData(SAMPLER, "SAMPLER", TSTRING);
515
516 if (cStatus) {
517 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
518 close();
519 return 1;
520 }
521
522
523 // Check for alternative column names.
524 if (cALFA) {
525 // ALFA data.
526 cALFAscan = 0;
527 cScanNo = 0;
528 if (cALFA_CIMA) {
529 findData(SCAN, "SCAN_ID", TINT);
530 if (cALFA_CIMA > 1) {
531 // Note that RECNUM increases by cNAxisTime per row.
532 findData(CYCLE, "RECNUM", TINT);
533 } else {
534 findData(CYCLE, "SUBSCAN", TINT);
535 }
536 } else if (cALFA_BD) {
537 findData(SCAN, "SCAN_NUMBER", TINT);
538 findData(CYCLE, "PATTERN_NUMBER", TINT);
539 }
540 } else {
541 readData(SCAN, 1, &cFirstScanNo);
542 }
543
544 cCycleNo = 0;
545 cLastUTC = 0.0;
546 for ( int i = 0 ; i < 4 ; i++ ) {
547 cGLastUTC[i] = 0.0 ;
548 cGLastScan[i] = -1 ;
549 cGCycleNo[i] = 0 ;
550 }
551
552 // Beam number, 1-relative by default.
553 cBeam_1rel = 1;
554 if (cALFA) {
555 // ALFA INPUT_ID, 0-relative (overrides BEAM column if present).
556 findData(BEAM, "INPUT_ID", TSHORT);
557 cBeam_1rel = 0;
558
559 } else if (cData[BEAM].colnum < 0) {
560 if (beamCRVAL) {
561 // There is a BEAM axis.
562 findData(BEAM, beamCRVAL, TDOUBLE);
563 } else {
564 // ms2sdfits output, 0-relative "feed" number.
565 findData(BEAM, "MAIN_FEED1", TSHORT);
566 cBeam_1rel = 0;
567 }
568 }
569
570 // IF number, 1-relative by default.
571 cIF_1rel = 1;
572 if (cALFA && cData[IF].colnum < 0) {
573 // ALFA data, 0-relative.
574 if (cALFA_CIMA > 1) {
575 findData(IF, "IFN", TSHORT);
576 } else {
577 findData(IF, "IFVAL", TSHORT);
578 }
579 cIF_1rel = 0;
580 }
581
582 // ms2sdfits writes a scalar "TSYS" column that averages the polarizations.
583 int colnum;
584 findCol("SYSCAL_TSYS", &colnum);
585 if (colnum > 0) {
586 // This contains the vector Tsys.
587 findData(TSYS, "SYSCAL_TSYS", TFLOAT);
588 }
589
590 // XPOLDATA?
591
592 if (cData[SCANRATE].colnum < 0) {
593 findData(SCANRATE, "FIELD_POINTING_DIR_RATE", TFLOAT);
594 }
595
596 if (cData[RESTFRQ].colnum < 0) {
597 findData(RESTFRQ, "RESTFREQ", TDOUBLE);
598 if (cData[RESTFRQ].colnum < 0) {
599 findData(RESTFRQ, "SPECTRAL_WINDOW_REST_FREQUENCY", TDOUBLE);
600 }
601 }
602
603 if (cData[OBJ_RA].colnum < 0) {
604 findData(OBJ_RA, "SOURCE_DIRECTION", TDOUBLE);
605 }
606 if (cData[OBJ_DEC].colnum < 0) {
607 findData(OBJ_DEC, "SOURCE_DIRECTION", TDOUBLE);
608 }
609
610 // REFBEAM?
611
612 if (cData[TCAL].colnum < 0) {
613 findData(TCAL, "SYSCAL_TCAL", TFLOAT);
614 } else if (cALFA_BD) {
615 // ALFA BDFITS has a different TCAL with 64 elements - kill it!
616 findData(TCAL, "NO NO NO", TFLOAT);
617 }
618
619 if (cALFA_BD) {
620 // ALFA BDFITS.
621 findData(AZIMUTH, "CRVAL2B", TFLOAT);
622 findData(ELEVATIO, "CRVAL3B", TFLOAT);
623 }
624
625 if (cALFA) {
626 // ALFA data.
627 findData(PARANGLE, "PARA_ANG", TFLOAT);
628 }
629
630 if (cData[TAMBIENT].colnum < 0) {
631 findData(TAMBIENT, "WEATHER_TEMPERATURE", TFLOAT);
632 }
633
634 if (cData[PRESSURE].colnum < 0) {
635 findData(PRESSURE, "WEATHER_PRESSURE", TFLOAT);
636 }
637
638 if (cData[HUMIDITY].colnum < 0) {
639 findData(HUMIDITY, "WEATHER_REL_HUMIDITY", TFLOAT);
640 }
641
642 if (cData[WINDSPEE].colnum < 0) {
643 findData(WINDSPEE, "WEATHER_WIND_SPEED", TFLOAT);
644 }
645
646 if (cData[WINDDIRE].colnum < 0) {
647 findData(WINDDIRE, "WEATHER_WIND_DIRECTION", TFLOAT);
648 }
649
650 // Determine which beams are present in the data.
651 if (cData[BEAM].colnum > 0) {
652 short *beamCol = new short[cNRow];
653 short beamNul = 1;
654 int anynul;
655 if (fits_read_col(cSDptr, TSHORT, cData[BEAM].colnum, 1, 1, cNRow,
656 &beamNul, beamCol, &anynul, &cStatus)) {
657 delete [] beamCol;
658 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
659 close();
660 return 1;
661 }
662
663 // Find the maximum beam number.
664 cNBeam = cBeam_1rel - 1;
665 for (int irow = 0; irow < cNRow; irow++) {
666 if (beamCol[irow] > cNBeam) {
667 cNBeam = beamCol[irow];
668 }
669
670 // Check validity.
671 if (beamCol[irow] < cBeam_1rel) {
672 delete [] beamCol;
673 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "SDFITS file contains invalid beam number.");
674 close();
675 return 1;
676 }
677 }
678
679 if (!cBeam_1rel) cNBeam++;
680
681 // Find all beams present in the data.
682 cBeams = new int[cNBeam];
683 for (int ibeam = 0; ibeam < cNBeam; ibeam++) {
684 cBeams[ibeam] = 0;
685 }
686
687 for (int irow = 0; irow < cNRow; irow++) {
688 cBeams[beamCol[irow] - cBeam_1rel] = 1;
689 }
690
691 delete [] beamCol;
692
693 } else {
694 // No BEAM column.
695 cNBeam = 1;
696 cBeams = new int[1];
697 cBeams[0] = 1;
698 }
699
700 // Passing back the address of the array allows PKSFITSreader::select() to
701 // modify its elements directly.
702 nBeam = cNBeam;
703 beams = cBeams;
704
705
706 // Determine which IFs are present in the data.
707 if (cData[IF].colnum > 0) {
708 short *IFCol = new short[cNRow];
709 short IFNul = 1;
710 int anynul;
711 if (fits_read_col(cSDptr, TSHORT, cData[IF].colnum, 1, 1, cNRow,
712 &IFNul, IFCol, &anynul, &cStatus)) {
713 delete [] IFCol;
714 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
715 close();
716 return 1;
717 }
718
719 // Find the maximum IF number.
720 cNIF = cIF_1rel - 1;
721 for (int irow = 0; irow < cNRow; irow++) {
722 if (IFCol[irow] > cNIF) {
723 cNIF = IFCol[irow];
724 }
725
726 // Check validity.
727 if (IFCol[irow] < cIF_1rel) {
728 delete [] IFCol;
729 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "SDFITS file contains invalid IF number.");
730 close();
731 return 1;
732 }
733 }
734
735 if (!cIF_1rel) cNIF++;
736
737 // Find all IFs present in the data.
738 cIFs = new int[cNIF];
739 cNChan = new int[cNIF];
740 cNPol = new int[cNIF];
741 cHaveXPol = new int[cNIF];
742 cGetXPol = 0;
743
744 for (int iIF = 0; iIF < cNIF; iIF++) {
745 cIFs[iIF] = 0;
746 cNChan[iIF] = 0;
747 cNPol[iIF] = 0;
748 cHaveXPol[iIF] = 0;
749 }
750
751 for (int irow = 0; irow < cNRow; irow++) {
752 int iIF = IFCol[irow] - cIF_1rel;
753 if (cIFs[iIF] == 0) {
754 cIFs[iIF] = 1;
755
756 // Find the axis lengths.
757 if (cHaveSpectra) {
758 if (cData[DATA].nelem < 0) {
759 // Variable dimension array.
760 if (readDim(DATA, irow+1, &cNAxes, cNAxis)) {
761 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
762 close();
763 return 1;
764 }
765 }
766
767 } else {
768 if (cData[DATAXED].colnum > 0) {
769 char dataxed[32];
770 readParm("DATAXED", TSTRING, dataxed);
771
772 sscanf(dataxed, "(%ld,%ld,%ld,%ld,%ld)", cNAxis, cNAxis+1,
773 cNAxis+2, cNAxis+3, cNAxis+4);
774 }
775 }
776
777 // Number of channels and polarizations.
778 cNChan[iIF] = cNAxis[cFreqAxis];
779 cNPol[iIF] = cNAxis[cStokesAxis];
780 cHaveXPol[iIF] = 0;
781
782 // Is cross-polarization data present?
783 if (cData[XPOLDATA].colnum > 0) {
784 // Check that it conforms.
785 int nAxis;
786 long nAxes[2];
787
788 if (readDim(XPOLDATA, irow+1, &nAxis, nAxes)) {
789 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE );
790 close();
791 return 1;
792 }
793
794 // Default is to get it if we have it.
795 if (nAxis == 2 &&
796 nAxes[0] == 2 &&
797 nAxes[1] == cNChan[iIF]) {
798 cGetXPol = cHaveXPol[iIF] = 1;
799 }
800 }
801 }
802 }
803
804 delete [] IFCol;
805
806 } else if ( numHDU > 1 ) {
807 Vector<Int> nIFs = nIFPerHDU() ;
808 cNIF = sum( nIFs ) ;
809 sprintf( cMsg, "cNIF = %d", cNIF ) ;
810 log( LogOrigin( className, methodName, WHERE ), LogIO::NORMAL, cMsg ) ;
811 cIFs = new int[cNIF] ;
812 cNChan = new int[cNIF];
813 cNPol = new int[cNIF];
814 cHaveXPol = new int[cNIF];
815 cGetXPol = 0;
816 for ( int iif = 0 ; iif < numHDU ; iif++ ) {
817 cIFs[iif] = 1 ;
818 cNChan[iif] = 0 ;
819 cNPol[iif] = 1 ;
820 cHaveXPol[iif] = 0 ;
821 }
822 Int cidx = 0 ;
823 for ( int ihdu = 0 ; ihdu < numHDU ; ihdu++ ) {
824 // Move
825 fits_movabs_hdu( cSDptr, hduId[ihdu], NULL, &cStatus ) ;
826
827 if (readDim(DATA, 1, &cNAxes, cNAxis)) {
828 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
829 close();
830 return 1;
831 }
832
833 // Number of channels and polarizations.
834 for ( int irow = 0 ; irow < nIFs[ihdu] ; irow++ ) {
835 cNChan[cidx+irow] = cNAxis[cFreqAxis];
836 cNPol[cidx+irow] = cNAxis[cStokesAxis];
837 cHaveXPol[cidx+irow] = 0;
838
839 // Is cross-polarization data present?
840 if (cData[XPOLDATA].colnum > 0) {
841 // Check that it conforms.
842 int nAxis;
843 long nAxes[2];
844
845 if (readDim(XPOLDATA, 1, &nAxis, nAxes)) {
846 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE );
847 close();
848 return 1;
849 }
850
851 // Default is to get it if we have it.
852 for ( int irow = 0 ; irow < nIFs[ihdu] ; irow++ ) {
853 if (nAxis == 2 &&
854 nAxes[0] == 2 &&
855 nAxes[1] == cNChan[cidx+irow]) {
856 cGetXPol = cHaveXPol[cidx+irow] = 1;
857 }
858 }
859 }
860 }
861 cidx += nIFs[ihdu] ;
862 }
863 } else {
864 // No IF column.
865 cNIF = 1;
866 cIFs = new int[1];
867 cIFs[0] = 1;
868
869 cNChan = new int[1];
870 cNPol = new int[1];
871 cHaveXPol = new int[1];
872 cGetXPol = 0;
873
874 // Number of channels and polarizations.
875 cNChan[0] = cNAxis[cFreqAxis];
876 cNPol[0] = cNAxis[cStokesAxis];
877 cHaveXPol[0] = 0;
878 }
879
880 if (cALFA && cALFA_CIMA < 2) {
881 // Older ALFA data labels each polarization as a separate IF.
882 cNPol[0] = cNIF;
883 cNIF = 1;
884 }
885
886 // For GBT data that stores spectra for each polarization in separate rows
887 if ( cData[STOKES].colnum > 0 ) {
888// int *stokesCol = new int[cNRow];
889 int stokesNul = 1;
890 int anynul;
891// if (fits_read_col(cSDptr, TINT, cData[STOKES].colnum, 1, 1, cNRow,
892// &stokesNul, stokesCol, &anynul, &cStatus)) {
893// delete [] stokesCol;
894// log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
895// close();
896// return 1;
897// }
898 vector<int *> stokesCols( numHDU ) ;
899 for ( int ihdu = 0 ; ihdu < numHDU ; ihdu++ ) {
900 fits_movabs_hdu( cSDptr, hduId[ihdu], NULL, &cStatus ) ;
901 stokesCols[ihdu] = new int[cNRowT[ihdu]] ;
902 fits_movabs_hdu( cSDptr, hduId[ihdu], NULL, &cStatus ) ;
903 if (fits_read_col(cSDptr, TINT, cData[STOKES].colnum, 1, 1, cNRowT[ihdu],
904 &stokesNul, stokesCols[ihdu] , &anynul, &cStatus)) {
905 for ( int jhdu = 0 ; jhdu <= ihdu ; jhdu++ ) {
906 delete [] stokesCols[jhdu] ;
907 }
908 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
909 close();
910 return 1;
911 }
912 }
913 fits_movabs_hdu( cSDptr, hduId[0], NULL, &cStatus ) ;
914
915 vector<int> pols ;
916// pols.push_back( stokesCol[0] ) ;
917// for ( int i = 0 ; i < cNRow ; i++ ) {
918// bool pmatch = false ;
919// for ( uint j = 0 ; j < pols.size() ; j++ ) {
920// if ( stokesCol[i] == pols[j] ) {
921// pmatch = true ;
922// break ;
923// }
924// }
925// if ( !pmatch ) {
926// pols.push_back( stokesCol[i] ) ;
927// }
928// }
929 pols.push_back( stokesCols[0][0] ) ;
930 for ( int k = 0 ; k < numHDU ; k++ ) {
931 for ( int i = 0 ; i < cNRowT[k] ; i++ ) {
932 bool pmatch = false ;
933 for ( uint j = 0 ; j < pols.size() ; j++ ) {
934 if ( stokesCols[k][i] == pols[j] ) {
935 pmatch = true ;
936 break ;
937 }
938 }
939 if ( !pmatch ) {
940 pols.push_back( stokesCols[k][i] ) ;
941 }
942 }
943 }
944
945 cPols = new int[pols.size()] ;
946 for ( uint i = 0 ; i < pols.size() ; i++ ) {
947 cPols[i] = pols[i] ;
948 }
949
950 for ( int i = 0 ; i < cNIF ; i++ ) {
951 cNPol[i] = pols.size() ;
952 }
953
954// delete [] stokesCol ;
955 for ( int k = 0 ; k < numHDU ; k++ ) {
956 delete [] stokesCols[k] ;
957 }
958 }
959
960 // Passing back the address of the array allows PKSFITSreader::select() to
961 // modify its elements directly.
962 nIF = cNIF;
963 IFs = cIFs;
964
965 nChan = cNChan;
966 nPol = cNPol;
967 haveXPol = cHaveXPol;
968
969
970 // Default channel range selection.
971 cStartChan = new int[cNIF];
972 cEndChan = new int[cNIF];
973 cRefChan = new int[cNIF];
974
975 for (int iIF = 0; iIF < cNIF; iIF++) {
976 cStartChan[iIF] = 1;
977 cEndChan[iIF] = cNChan[iIF];
978 cRefChan[iIF] = cNChan[iIF]/2 + 1;
979 }
980
981 // Default is to get it if we have it.
982 cGetSpectra = cHaveSpectra;
983
984
985 // Are baseline parameters present?
986 cHaveBase = 0;
987 if (cData[BASELIN].colnum) {
988 // Check that it conforms.
989 int nAxis, status = 0;
990 long nAxes[2];
991
992 if (fits_read_tdim(cSDptr, cData[BASELIN].colnum, 2, &nAxis, nAxes,
993 &status) == 0) {
994 cHaveBase = (nAxis == 2);
995 }
996 }
997 haveBase = cHaveBase;
998
999
1000 // Is extra system calibration data available?
1001 cExtraSysCal = 0;
1002 for (int iparm = REFBEAM; iparm < NDATA; iparm++) {
1003 if (cData[iparm].colnum >= 0) {
1004 cExtraSysCal = 1;
1005 break;
1006 }
1007 }
1008
1009 extraSysCal = cExtraSysCal;
1010
1011
1012 // Extras for ALFA data.
1013 cALFAacc = 0.0f;
1014 if (cALFA_CIMA > 1) {
1015 // FFTs per second when the Mock correlator operates in RFI blanking mode.
1016 readData("PHFFTACC", TFLOAT, 0, &cALFAacc);
1017 }
1018
1019
1020 cRow = 0;
1021 cTimeIdx = cNAxisTime;
1022
1023 return 0;
1024}
1025
1026//---------------------------------------------------- GBTFITSreader::getHeader
1027
1028// Get parameters describing the data.
1029
1030int GBTFITSreader::getHeader(
1031 char observer[32],
1032 char project[32],
1033 char telescope[32],
1034 double antPos[3],
1035 char obsMode[32],
1036 char bunit[32],
1037 float &equinox,
1038 char radecsys[32],
1039 char dopplerFrame[32],
1040 char datobs[32],
1041 double &utc,
1042 double &refFreq,
1043 double &bandwidth)
1044{
1045 const string methodName = "getHeader()" ;
1046
1047 // Has the file been opened?
1048 if (!cSDptr) {
1049 return 1;
1050 }
1051
1052 // Read parameter values.
1053 readParm("OBSERVER", TSTRING, observer); // Shared.
1054 readParm("PROJID", TSTRING, project); // Shared.
1055 readParm("TELESCOP", TSTRING, telescope); // Core.
1056
1057 antPos[0] = 0.0;
1058 antPos[1] = 0.0;
1059 antPos[2] = 0.0;
1060 if (readParm("ANTENNA_POSITION", TDOUBLE, antPos)) {
1061 readParm("OBSGEO-X", TDOUBLE, antPos); // Additional.
1062 readParm("OBSGEO-Y", TDOUBLE, antPos + 1); // Additional.
1063 readParm("OBSGEO-Z", TDOUBLE, antPos + 2); // Additional.
1064 }
1065
1066 if (antPos[0] == 0.0) {
1067 if (cGBT) {
1068 readParm( "SITELONG", TDOUBLE, antPos ) ;
1069 readParm( "SITELAT", TDOUBLE, antPos+1 ) ;
1070 readParm( "SITEELEV", TDOUBLE, antPos+2 ) ;
1071 Vector<Double> pos( 2 ) ;
1072 pos[0] = Double( antPos[0] ) ;
1073 pos[1] = Double( antPos[1] ) ;
1074 MPosition mp( MVPosition( Quantity( Double(antPos[2]), "m" ),
1075 Quantum< Vector<Double> >( pos, "deg" ) ),
1076 MPosition::WGS84 ) ;
1077 MeasFrame mFrame( mp ) ;
1078 MVPosition mvp ;
1079 mFrame.getITRF( mvp ) ;
1080 Vector<Double> antPosVec = mvp.getValue() ;
1081 antPos[0] = (double)antPosVec[0] ;
1082 antPos[1] = (double)antPosVec[1] ;
1083 antPos[2] = (double)antPosVec[2] ;
1084 }
1085 else if (strncmp(telescope, "ATPKS", 5) == 0) {
1086 // Parkes coordinates.
1087 antPos[0] = -4554232.087;
1088 antPos[1] = 2816759.046;
1089 antPos[2] = -3454035.950;
1090 } else if (strncmp(telescope, "ATMOPRA", 7) == 0) {
1091 // Mopra coordinates.
1092 antPos[0] = -4682768.630;
1093 antPos[1] = 2802619.060;
1094 antPos[2] = -3291759.900;
1095 } else if (strncmp(telescope, "ARECIBO", 7) == 0) {
1096 // Arecibo coordinates.
1097 antPos[0] = 2390486.900;
1098 antPos[1] = -5564731.440;
1099 antPos[2] = 1994720.450;
1100 }
1101 }
1102
1103 readData(OBSMODE, 1, obsMode); // Shared.
1104
1105 // Brightness unit.
1106 if (cData[DATAXED].colnum >= 0) {
1107 strcpy(bunit, "Jy");
1108 } else {
1109 strcpy(bunit, cData[DATA].units);
1110 }
1111
1112 if (strcmp(bunit, "JY") == 0) {
1113 bunit[1] = 'y';
1114 } else if (strcmp(bunit, "JY/BEAM") == 0) {
1115 strcpy(bunit, "Jy/beam");
1116 }
1117
1118 readParm("EQUINOX", TFLOAT, &equinox); // Shared.
1119 if (cStatus == 405) {
1120 // EQUINOX was written as string value in early versions.
1121 cStatus = 0;
1122 char strtmp[32];
1123 readParm("EQUINOX", TSTRING, strtmp);
1124 sscanf(strtmp, "%f", &equinox);
1125 }
1126
1127 if (readParm("RADESYS", TSTRING, radecsys)) { // Additional.
1128 if (readParm("RADECSYS", TSTRING, radecsys)) { // Additional.
1129 strcpy(radecsys, "");
1130 }
1131 }
1132
1133 if (readParm("SPECSYS", TSTRING, dopplerFrame)) { // Additional.
1134 // Fallback value.
1135 strcpy(dopplerFrame, "TOPOCENT");
1136
1137 // Look for VELFRAME, written by earlier versions of Livedata.
1138 //
1139 // Added few more codes currently (as of 2009 Oct) used in the GBT
1140 // SDFITS (based io_sdfits_define.pro of GBTIDL). - TT
1141 if (readParm("VELFRAME", TSTRING, dopplerFrame)) { // Additional.
1142 // No, try digging it out of the CTYPE card (AIPS convention).
1143 char keyw[9], ctype[9];
1144 sprintf(keyw, "CTYPE%ld", cFreqAxis+1);
1145 readParm(keyw, TSTRING, ctype);
1146
1147 if (strncmp(ctype, "FREQ-", 5) == 0) {
1148 strcpy(dopplerFrame, ctype+5);
1149 if (strcmp(dopplerFrame, "LSR") == 0) {
1150 // LSR unqualified usually means LSR (kinematic).
1151 strcpy(dopplerFrame, "LSRK");
1152 } else if (strcmp(dopplerFrame, "LSD") == 0) {
1153 // LSR as a dynamical defintion
1154 strcpy(dopplerFrame, "LSRD");
1155 } else if (strcmp(dopplerFrame, "HEL") == 0) {
1156 // Almost certainly barycentric.
1157 strcpy(dopplerFrame, "BARYCENT");
1158 } else if (strcmp(dopplerFrame, "BAR") == 0) {
1159 // barycentric.
1160 strcpy(dopplerFrame, "BARYCENT");
1161 } else if (strcmp(dopplerFrame, "OBS") == 0) {
1162 strcpy(dopplerFrame, "TOPO");
1163 } else if (strcmp(dopplerFrame, "GEO") == 0) {
1164 // geocentric
1165 strcpy(dopplerFrame, "GEO");
1166 } else if (strcmp(dopplerFrame, "GAL") == 0) {
1167 // galactic
1168 strcpy(dopplerFrame, "GAL");
1169 } else if (strcmp(dopplerFrame, "LGR") == 0) {
1170 // Local group
1171 strcpy(dopplerFrame, "LGROUP");
1172 } else if (strcmp(dopplerFrame, "CMB") == 0) {
1173 // Cosimic Microwave Backgroup
1174 strcpy(dopplerFrame, "CMB");
1175 }
1176 } else {
1177 strcpy(dopplerFrame, "");
1178 }
1179 }
1180
1181 // always LSRK if GBT data
1182 if (cGBT) {
1183 strcpy(dopplerFrame, "LSRK");
1184 }
1185
1186
1187 // Translate to FITS standard names.
1188 if (strncmp(dopplerFrame, "TOP", 3) == 0) {
1189 strcpy(dopplerFrame, "TOPOCENT");
1190 } else if (strncmp(dopplerFrame, "GEO", 3) == 0) {
1191 strcpy(dopplerFrame, "GEOCENTR");
1192 } else if (strncmp(dopplerFrame, "HEL", 3) == 0) {
1193 strcpy(dopplerFrame, "HELIOCEN");
1194 } else if (strncmp(dopplerFrame, "BARY", 4) == 0) {
1195 strcpy(dopplerFrame, "BARYCENT");
1196 } else if (strncmp(dopplerFrame, "GAL", 3) == 0) {
1197 strcpy(dopplerFrame, "GALACTOC");
1198 } else if (strncmp(dopplerFrame, "LGROUP", 6) == 0) {
1199 strcpy(dopplerFrame, "LOCALGRP");
1200 } else if (strncmp(dopplerFrame, "CMB", 3) == 0) {
1201 strcpy(dopplerFrame, "CMBDIPOL");
1202 }
1203 }
1204
1205 if (cStatus) {
1206 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1207 return 1;
1208 }
1209
1210 // Get parameters from first row of table.
1211 readTime(1, 1, datobs, utc);
1212 readData(FqRefVal, 1, &refFreq);
1213 readParm("BANDWID", TDOUBLE, &bandwidth); // Core.
1214
1215 if (cStatus) {
1216 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1217 return 1;
1218 }
1219
1220 return 0;
1221}
1222
1223//-------------------------------------------------- GBTFITSreader::getFreqInfo
1224
1225// Get frequency parameters for each IF.
1226
1227int GBTFITSreader::getFreqInfo(
1228 int &nIF,
1229 double* &startFreq,
1230 double* &endFreq)
1231{
1232 const string methodName = "getFreqInfo()" ;
1233
1234 float fqRefPix;
1235 double fqDelt, fqRefVal;
1236
1237 nIF = cNIF;
1238 startFreq = new double[nIF];
1239 endFreq = new double[nIF];
1240
1241 if (cData[IF].colnum > 0) {
1242 short *IFCol = new short[cNRow];
1243 short IFNul = 1;
1244 int anynul;
1245 if (fits_read_col(cSDptr, TSHORT, cData[IF].colnum, 1, 1, cNRow,
1246 &IFNul, IFCol, &anynul, &cStatus)) {
1247 delete [] IFCol;
1248 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1249 close();
1250 return 1;
1251 }
1252
1253 for (int iIF = 0; iIF < nIF; iIF++) {
1254 if (cIFs[iIF]) {
1255 // Find the first occurrence of this IF in the table.
1256 int IFno = iIF + cIF_1rel;
1257 for (int irow = 0; irow < cNRow;) {
1258 if (IFCol[irow++] == IFno) {
1259 readData(FqRefPix, irow, &fqRefPix);
1260 readData(FqRefVal, irow, &fqRefVal);
1261 readData(FqDelt, irow, &fqDelt);
1262
1263 if (cALFA_BD) {
1264 unsigned char invert;
1265 readData("UPPERSB", TBYTE, irow, &invert);
1266
1267 if (invert) {
1268 fqDelt = -fqDelt;
1269 }
1270 }
1271
1272 startFreq[iIF] = fqRefVal + ( 1 - fqRefPix) * fqDelt;
1273 endFreq[iIF] = fqRefVal + (cNChan[iIF] - fqRefPix) * fqDelt;
1274
1275 break;
1276 }
1277 }
1278
1279 } else {
1280 startFreq[iIF] = 0.0;
1281 endFreq[iIF] = 0.0;
1282 }
1283 }
1284
1285 delete [] IFCol;
1286
1287 } else {
1288 // No IF column, read the first table entry.
1289 readData(FqRefPix, 1, &fqRefPix);
1290 readData(FqRefVal, 1, &fqRefVal);
1291 readData(FqDelt, 1, &fqDelt);
1292
1293 startFreq[0] = fqRefVal + ( 1 - fqRefPix) * fqDelt;
1294 endFreq[0] = fqRefVal + (cNChan[0] - fqRefPix) * fqDelt;
1295 }
1296
1297 return cStatus;
1298}
1299
1300//---------------------------------------------------- GBTFITSreader::findRange
1301
1302// Find the range of the data in time and position.
1303
1304int GBTFITSreader::findRange(
1305 int &nRow,
1306 int &nSel,
1307 char dateSpan[2][32],
1308 double utcSpan[2],
1309 double* &positions)
1310{
1311 const string methodName = "findRange()" ;
1312
1313 // Has the file been opened?
1314 if (!cSDptr) {
1315 return 1;
1316 }
1317
1318 nRow = cNRow;
1319
1320 // Find the number of rows selected.
1321 short *sel = new short[cNRow];
1322 for (int irow = 0; irow < cNRow; irow++) {
1323 sel[irow] = 1;
1324 }
1325
1326 int anynul;
1327 if (cData[BEAM].colnum > 0) {
1328 short *beamCol = new short[cNRow];
1329 short beamNul = 1;
1330 if (fits_read_col(cSDptr, TSHORT, cData[BEAM].colnum, 1, 1, cNRow,
1331 &beamNul, beamCol, &anynul, &cStatus)) {
1332 delete [] beamCol;
1333 delete [] sel;
1334 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1335 return 1;
1336 }
1337
1338 for (int irow = 0; irow < cNRow; irow++) {
1339 if (!cBeams[beamCol[irow]-cBeam_1rel]) {
1340 sel[irow] = 0;
1341 }
1342 }
1343
1344 delete [] beamCol;
1345 }
1346
1347 if (cData[IF].colnum > 0) {
1348 short *IFCol = new short[cNRow];
1349 short IFNul = 1;
1350 if (fits_read_col(cSDptr, TSHORT, cData[IF].colnum, 1, 1, cNRow,
1351 &IFNul, IFCol, &anynul, &cStatus)) {
1352 delete [] IFCol;
1353 delete [] sel;
1354 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1355 return 1;
1356 }
1357
1358 for (int irow = 0; irow < cNRow; irow++) {
1359 if (!cIFs[IFCol[irow]-cIF_1rel]) {
1360 sel[irow] = 0;
1361 }
1362 }
1363
1364 delete [] IFCol;
1365 }
1366
1367 nSel = 0;
1368 for (int irow = 0; irow < cNRow; irow++) {
1369 nSel += sel[irow];
1370 }
1371
1372
1373 // Find the time range assuming the data is in chronological order.
1374 readTime(1, 1, dateSpan[0], utcSpan[0]);
1375 readTime(cNRow, cNAxisTime, dateSpan[1], utcSpan[1]);
1376
1377
1378 // Retrieve positions for selected data.
1379 int isel = 0;
1380 positions = new double[2*nSel];
1381
1382 if (cCoordSys == 1) {
1383 // Horizontal (Az,El).
1384 if (cData[AZIMUTH].colnum < 0 ||
1385 cData[ELEVATIO].colnum < 0) {
1386 log(LogOrigin( className, methodName, WHERE ), LogIO::WARN, "Azimuth/elevation information absent.");
1387 cStatus = -1;
1388
1389 } else {
1390 float *az = new float[cNRow];
1391 float *el = new float[cNRow];
1392 readCol(AZIMUTH, az);
1393 readCol(ELEVATIO, el);
1394
1395 if (!cStatus) {
1396 for (int irow = 0; irow < cNRow; irow++) {
1397 if (sel[irow]) {
1398 positions[isel++] = az[irow] * D2R;
1399 positions[isel++] = el[irow] * D2R;
1400 }
1401 }
1402 }
1403
1404 delete [] az;
1405 delete [] el;
1406 }
1407
1408 } else if (cCoordSys == 3) {
1409 // ZPA-EL.
1410 if (cData[BEAM].colnum < 0 ||
1411 cData[FOCUSROT].colnum < 0 ||
1412 cData[ELEVATIO].colnum < 0) {
1413 log(LogOrigin( className, methodName, WHERE ), LogIO::WARN, "ZPA/elevation information absent.");
1414 cStatus = -1;
1415
1416 } else {
1417 short *beam = new short[cNRow];
1418 float *rot = new float[cNRow];
1419 float *el = new float[cNRow];
1420 readCol(BEAM, beam);
1421 readCol(FOCUSROT, rot);
1422 readCol(ELEVATIO, el);
1423
1424 if (!cStatus) {
1425 for (int irow = 0; irow < cNRow; irow++) {
1426 if (sel[irow]) {
1427 Int beamNo = beam[irow];
1428 Double zpa = rot[irow];
1429 if (beamNo > 1) {
1430 // Beam geometry for the Parkes multibeam.
1431 if (beamNo < 8) {
1432 zpa += -60.0 + 60.0*(beamNo-2);
1433 } else {
1434 zpa += -90.0 + 60.0*(beamNo-8);
1435 }
1436
1437 if (zpa < -180.0) {
1438 zpa += 360.0;
1439 } else if (zpa > 180.0) {
1440 zpa -= 360.0;
1441 }
1442 }
1443
1444 positions[isel++] = zpa * D2R;
1445 positions[isel++] = el[irow] * D2R;
1446 }
1447 }
1448 }
1449
1450 delete [] beam;
1451 delete [] rot;
1452 delete [] el;
1453 }
1454
1455 } else {
1456 double *ra = new double[cNRow];
1457 double *dec = new double[cNRow];
1458 readCol(RA, ra);
1459 readCol(DEC, dec);
1460
1461 if (cStatus) {
1462 delete [] ra;
1463 delete [] dec;
1464 goto cleanup;
1465 }
1466
1467 if (cALFA_BD) {
1468 for (int irow = 0; irow < cNRow; irow++) {
1469 // Convert hours to degrees.
1470 ra[irow] *= 15.0;
1471 }
1472 }
1473
1474 if (cCoordSys == 0) {
1475 // Equatorial (RA,Dec).
1476 for (int irow = 0; irow < cNRow; irow++) {
1477 if (sel[irow]) {
1478 positions[isel++] = ra[irow] * D2R;
1479 positions[isel++] = dec[irow] * D2R;
1480 }
1481 }
1482
1483 } else if (cCoordSys == 2) {
1484 // Feed-plane.
1485 if (cData[OBJ_RA].colnum < 0 ||
1486 cData[OBJ_DEC].colnum < 0 ||
1487 cData[PARANGLE].colnum < 0 ||
1488 cData[FOCUSROT].colnum < 0) {
1489 log( LogOrigin( className, methodName, WHERE ), LogIO::WARN,
1490 "Insufficient information to compute feed-plane\n"
1491 " coordinates.");
1492 cStatus = -1;
1493
1494 } else {
1495 double *srcRA = new double[cNRow];
1496 double *srcDec = new double[cNRow];
1497 float *par = new float[cNRow];
1498 float *rot = new float[cNRow];
1499
1500 readCol(OBJ_RA, srcRA);
1501 readCol(OBJ_DEC, srcDec);
1502 readCol(PARANGLE, par);
1503 readCol(FOCUSROT, rot);
1504
1505 if (!cStatus) {
1506 for (int irow = 0; irow < cNRow; irow++) {
1507 if (sel[irow]) {
1508 // Convert to feed-plane coordinates.
1509 Double dist, pa;
1510 distPA(ra[irow]*D2R, dec[irow]*D2R, srcRA[irow]*D2R,
1511 srcDec[irow]*D2R, dist, pa);
1512
1513 Double spin = (par[irow] + rot[irow])*D2R - pa;
1514 if (spin > 2.0*PI) spin -= 2.0*PI;
1515 Double squint = PI/2.0 - dist;
1516
1517 positions[isel++] = spin;
1518 positions[isel++] = squint;
1519 }
1520 }
1521 }
1522
1523 delete [] srcRA;
1524 delete [] srcDec;
1525 delete [] par;
1526 delete [] rot;
1527 }
1528 }
1529
1530 delete [] ra;
1531 delete [] dec;
1532 }
1533
1534cleanup:
1535 delete [] sel;
1536
1537 if (cStatus) {
1538 nSel = 0;
1539 delete [] positions;
1540 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1541 cStatus = 0;
1542 return 1;
1543 }
1544
1545 return 0;
1546}
1547
1548
1549//--------------------------------------------------------- GBTFITSreader::read
1550
1551// Read the next data record.
1552
1553int GBTFITSreader::read(
1554 MBrecord &mbrec)
1555{
1556 const string methodName = "read()" ;
1557
1558 // for multiple binary table
1559 long cRowR = -1 ;
1560
1561 // Has the file been opened?
1562 if (!cSDptr) {
1563 return 1;
1564 }
1565 // Find the next selected beam and IF.
1566 short iBeam = 0, iIF = 0;
1567 int iPol = -1 ;
1568 while (1) {
1569 if (++cTimeIdx > cNAxisTime) {
1570 if (++cRow > cNRow) break;
1571 cTimeIdx = 1;
1572 }
1573
1574 //cRowR = cRow ;
1575 cRowR = cRowRelative( cRow ) ;
1576
1577 if (cData[BEAM].colnum > 0) {
1578 readData(BEAM, cRowR, &iBeam);
1579
1580 // Convert to 0-relative.
1581 if (cBeam_1rel) iBeam--;
1582 }
1583
1584
1585 if (cBeams[iBeam]) {
1586 if (cData[IF].colnum > 0) {
1587 readData(IF, cRowR, &iIF);
1588
1589 // Convert to 0-relative.
1590 if (cIF_1rel) iIF--;
1591 }
1592
1593 if (cIFs[iIF]) {
1594 if (cALFA) {
1595 // ALFA data, check for calibration data.
1596 char chars[32];
1597 readData(OBSMODE, cRowR, chars);
1598 if (strcmp(chars, "DROP") == 0) {
1599 // Completely flagged integration.
1600 continue;
1601
1602 } else if (strcmp(chars, "CAL") == 0) {
1603 sReset = 1;
1604 if (cALFA_CIMA > 1) {
1605 for (short iPol = 0; iPol < cNPol[iIF]; iPol++) {
1606 alfaCal(iBeam, iIF, iPol);
1607 }
1608 continue;
1609 } else {
1610 // iIF is really the polarization in older ALFA data.
1611 alfaCal(iBeam, 0, iIF);
1612 continue;
1613 }
1614
1615 } else {
1616 // Reset for the next CAL record.
1617 if (sReset) {
1618 for (short iPol = 0; iPol < cNPol[iIF]; iPol++) {
1619 sALFAcalNon[iBeam][iPol] = 0;
1620 sALFAcalNoff[iBeam][iPol] = 0;
1621 sALFAcalOn[iBeam][iPol] = 0.0f;
1622 sALFAcalOff[iBeam][iPol] = 0.0f;
1623 }
1624 sReset = 0;
1625
1626 sprintf(cMsg, "ALFA cal factors for beam %d: %.3e, %.3e",
1627 iBeam+1, sALFAcal[iBeam][0], sALFAcal[iBeam][1]);
1628 log(LogOrigin( className, methodName, WHERE ), LogIO::NORMAL, cMsg);
1629 //logMsg(cMsg);
1630 }
1631 }
1632 }
1633
1634 // for GBT SDFITS
1635 if (cData[STOKES].colnum > 0 ) {
1636 readData(STOKES, cRowR, &iPol ) ;
1637 for ( int i = 0 ; i < cNPol[iIF] ; i++ ) {
1638 if ( cPols[i] == iPol ) {
1639 iPol = i ;
1640 break ;
1641 }
1642 }
1643 }
1644 break;
1645 }
1646 }
1647 }
1648
1649 // EOF?
1650 if (cRow > cNRow) {
1651 return -1;
1652 }
1653
1654 int idx_GBT = 0 ;
1655 if ( numHDU > 1 ) {
1656 char sampler[10];
1657 readData(SAMPLER,cRowR,sampler);
1658 readData(RESTFRQ, cRowR, &mbrec.restFreq);
1659 int iiIF = 0 ;
1660 int hduid ;
1661 fits_get_hdu_num( cSDptr, &hduid ) ;
1662 while ( iiIF < numHDU && hduId[iiIF] != hduid ) {
1663 //iiIF++ ;
1664 idx_GBT += cNumKeys[iiIF++] ;
1665 }
1666 char sarf[48] ;
1667 sprintf( sarf, "%d/%s/%lf", iiIF, sampler, mbrec.restFreq ) ;
1668 String sr = String( sarf ) ;
1669 while ( idx_GBT < (int)cRefValKey.size() && sr != cRefValKey[idx_GBT] )
1670 idx_GBT++ ;
1671 iIF = (short)cIFNO[idx_GBT] ;
1672 sprintf( cMsg, "iIF = %d (key=%s), cRow = %ld" , iIF, sarf, cRow ) ;
1673 log( LogOrigin( className, methodName, WHERE ), LogIO::DEBUGGING, cMsg ) ;
1674 }
1675
1676 if (cALFA) {
1677 int scanNo;
1678 readData(SCAN, cRowR, &scanNo);
1679 if (scanNo != cALFAscan) {
1680 cScanNo++;
1681 cALFAscan = scanNo;
1682 }
1683 mbrec.scanNo = cScanNo;
1684
1685 } else {
1686 readData(SCAN, cRowR, &mbrec.scanNo);
1687
1688 // Ensure that scan number is 1-relative.
1689 mbrec.scanNo -= (cFirstScanNo - 1);
1690 }
1691
1692 // Times.
1693 char datobs[32];
1694 readTime(cRowR, cTimeIdx, datobs, mbrec.utc);
1695 strcpy(mbrec.datobs, datobs);
1696
1697 if (cData[CYCLE].colnum > 0) {
1698 readData(CYCLE, cRowR, &mbrec.cycleNo);
1699 mbrec.cycleNo += cTimeIdx - 1;
1700 if (cALFA_BD) mbrec.cycleNo++;
1701 } else {
1702 // Cycle number not recorded, must do our own bookkeeping.
1703 if (mbrec.utc != cLastUTC) {
1704 mbrec.cycleNo = ++cCycleNo;
1705 cLastUTC = mbrec.utc;
1706 }
1707 }
1708
1709 if ( iPol != -1 ) {
1710 if ( mbrec.scanNo != cGLastScan[iPol] ) {
1711 cGLastScan[iPol] = mbrec.scanNo ;
1712 cGCycleNo[iPol] = 0 ;
1713 mbrec.cycleNo = ++cGCycleNo[iPol] ;
1714 }
1715 else {
1716 mbrec.cycleNo = ++cGCycleNo[iPol] ;
1717 }
1718 }
1719
1720 readData(EXPOSURE, cRowR, &mbrec.exposure);
1721
1722 // Source identification.
1723 readData(OBJECT, cRowR, mbrec.srcName);
1724
1725 if ( iPol != -1 ) {
1726 char obsmode[32] ;
1727 readData( OBSMODE, cRowR, obsmode ) ;
1728 char sig[1] ;
1729 char cal[1] ;
1730 readData( SIG, cRowR, sig ) ;
1731 readData( CAL, cRowR, cal ) ;
1732 if ( strstr( obsmode, "PSWITCH" ) != NULL ) {
1733 // position switch
1734 //strcat( mbrec.srcName, "_p" ) ;
1735 if ( strstr( obsmode, "PSWITCHON" ) != NULL ) {
1736 //strcat( mbrec.srcName, "s" ) ;
1737 mbrec.srcType = SrcType::PSON ;
1738 }
1739 else if ( strstr( obsmode, "PSWITCHOFF" ) != NULL ) {
1740 //strcat( mbrec.srcName, "r" ) ;
1741 mbrec.srcType = SrcType::PSOFF ;
1742 }
1743 }
1744 else if ( strstr( obsmode, "Nod" ) != NULL ) {
1745 // nod
1746 //strcat( mbrec.srcName, "_n" ) ;
1747 mbrec.srcType = SrcType::NOD ;
1748 if ( sig[0] == 'T' ) {
1749 //strcat( mbrec.srcName, "s" ) ;
1750 }
1751 else {
1752 //strcat( mbrec.srcName, "r" ) ;
1753 }
1754 }
1755 else if ( strstr( obsmode, "FSWITCH" ) != NULL ) {
1756 // frequency switch
1757 //strcat( mbrec.srcName, "_f" ) ;
1758 if ( sig[0] == 'T' ) {
1759 //strcat( mbrec.srcName, "s" ) ;
1760 mbrec.srcType = SrcType::FSON ;
1761 }
1762 else {
1763 //strcat( mbrec.srcName, "r" ) ;
1764 mbrec.srcType = SrcType::FSOFF ;
1765 }
1766 }
1767 if ( cal[0] == 'T' ) {
1768 //strcat( mbrec.srcName, "c" ) ;
1769 if ( mbrec.srcType == SrcType::PSON )
1770 mbrec.srcType = SrcType::PONCAL ;
1771 else if ( mbrec.srcType == SrcType::PSOFF )
1772 mbrec.srcType = SrcType::POFFCAL ;
1773 else if ( mbrec.srcType == SrcType::NOD )
1774 mbrec.srcType = SrcType::NODCAL ;
1775 else if ( mbrec.srcType == SrcType::FSON )
1776 mbrec.srcType = SrcType::FONCAL ;
1777 else if ( mbrec.srcType == SrcType::FSOFF )
1778 mbrec.srcType = SrcType::FOFFCAL ;
1779 }
1780 else {
1781 //strcat( mbrec.srcName, "o" ) ;
1782 }
1783 }
1784
1785 readData(OBJ_RA, cRowR, &mbrec.srcRA);
1786 if (strcmp(cData[OBJ_RA].name, "OBJ-RA") == 0) {
1787 mbrec.srcRA *= D2R;
1788 }
1789
1790 if (strcmp(cData[OBJ_DEC].name, "OBJ-DEC") == 0) {
1791 readData(OBJ_DEC, cRowR, &mbrec.srcDec);
1792 mbrec.srcDec *= D2R;
1793 }
1794
1795 // Line rest frequency (Hz).
1796 readData(RESTFRQ, cRowR, &mbrec.restFreq);
1797 if (mbrec.restFreq == 0.0 && cALFA_BD) {
1798 mbrec.restFreq = 1420.40575e6;
1799 }
1800
1801 // Observation mode.
1802 readData(OBSMODE, cRowR, mbrec.obsType);
1803
1804 // Beam-dependent parameters.
1805 mbrec.beamNo = iBeam + 1;
1806
1807 readData(RA, cRowR, &mbrec.ra);
1808 readData(DEC, cRowR, &mbrec.dec);
1809 mbrec.ra *= D2R;
1810 mbrec.dec *= D2R;
1811
1812 if (cALFA_BD) mbrec.ra *= 15.0;
1813
1814 float scanrate[2];
1815 readData(SCANRATE, cRowR, &scanrate);
1816 if (strcmp(cData[SCANRATE].name, "SCANRATE") == 0) {
1817 mbrec.raRate = scanrate[0] * D2R;
1818 mbrec.decRate = scanrate[1] * D2R;
1819 }
1820 mbrec.paRate = 0.0f;
1821
1822 // IF-dependent parameters.
1823 int startChan = cStartChan[iIF];
1824 int endChan = cEndChan[iIF];
1825 int refChan = cRefChan[iIF];
1826
1827 // Allocate data storage.
1828 int nChan = abs(endChan - startChan) + 1;
1829 int nPol = cNPol[iIF];
1830
1831 if ( cData[STOKES].colnum > 0 )
1832 nPol = 1 ;
1833
1834 if (cGetSpectra || cGetXPol) {
1835 int nxpol = cGetXPol ? 2*nChan : 0;
1836 mbrec.allocate(0, nChan*nPol, nxpol);
1837 }
1838
1839 mbrec.nIF = 1;
1840 mbrec.IFno[0] = iIF + 1;
1841 mbrec.nChan[0] = nChan;
1842 mbrec.nPol[0] = nPol;
1843 mbrec.polNo = iPol ;
1844
1845 readData(FqRefPix, cRowR, mbrec.fqRefPix);
1846 readData(FqRefVal, cRowR, mbrec.fqRefVal);
1847 readData(FqDelt, cRowR, mbrec.fqDelt);
1848
1849 // Convert frequency reference to LSRK if GBT TOPO
1850 if ( cGBT ) {
1851 char sampler[10];
1852 readData(SAMPLER,cRowR,sampler);
1853 //Double fqrefval = getRefValLSR( cRow, idx_GBT, sampler, mbrec.restFreq ) ;
1854 sprintf( cMsg, "Row %ld (%s): reference frequency %lf", cRow, sampler, cRefVal[idx_GBT] ) ;
1855 log( LogOrigin( className, methodName, WHERE ), LogIO::DEBUGGING, cMsg ) ;
1856
1857 if ( idx_GBT == (int)cRefValKey.size() ) {
1858 log( LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Failed to get reference frequency value" ) ;
1859 mbrec.fqRefVal[0] = 0.0 ;
1860 }
1861 else {
1862 mbrec.fqRefVal[0] = (double)cRefVal[idx_GBT] ;
1863 }
1864 }
1865
1866 if (cALFA_BD) {
1867 unsigned char invert;
1868 int anynul, colnum;
1869 findCol("UPPERSB", &colnum);
1870 fits_read_col(cSDptr, TBYTE, colnum, cRow, 1, 1, 0, &invert, &anynul,
1871 &cStatus);
1872
1873 if (invert) {
1874 mbrec.fqDelt[0] = -mbrec.fqDelt[0];
1875 }
1876 }
1877
1878 if (cStatus) {
1879 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1880 return 1;
1881 }
1882
1883 // Adjust for channel selection.
1884 if (mbrec.fqRefPix[0] != refChan) {
1885 mbrec.fqRefVal[0] += (refChan - mbrec.fqRefPix[0]) * mbrec.fqDelt[0];
1886 mbrec.fqRefPix[0] = refChan;
1887 }
1888
1889 if (endChan < startChan) {
1890 mbrec.fqDelt[0] = -mbrec.fqDelt[0];
1891 }
1892
1893 // The data may only have a scalar Tsys value.
1894 mbrec.tsys[0][0] = 0.0f;
1895 mbrec.tsys[0][1] = 0.0f;
1896 if (cData[TSYS].nelem >= nPol) {
1897 readData(TSYS, cRowR, mbrec.tsys[0]);
1898 }
1899
1900 for (int j = 0; j < 2; j++) {
1901 mbrec.calfctr[0][j] = 0.0f;
1902 }
1903 if (cData[CALFCTR].colnum > 0) {
1904 readData(CALFCTR, cRowR, mbrec.calfctr);
1905 }
1906
1907 if (cHaveBase) {
1908 mbrec.haveBase = 1;
1909 readData(BASELIN, cRowR, mbrec.baseLin);
1910 readData(BASESUB, cRowR, mbrec.baseSub);
1911 } else {
1912 mbrec.haveBase = 0;
1913 }
1914
1915 if (cStatus) {
1916 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1917 return 1;
1918 }
1919
1920 // Read data, sectioning and transposing it in the process.
1921 long *blc = new long[cNAxes+1];
1922 long *trc = new long[cNAxes+1];
1923 long *inc = new long[cNAxes+1];
1924 for (int iaxis = 0; iaxis <= cNAxes; iaxis++) {
1925 blc[iaxis] = 1;
1926 trc[iaxis] = 1;
1927 inc[iaxis] = 1;
1928 }
1929
1930 blc[cFreqAxis] = std::min(startChan, endChan);
1931 trc[cFreqAxis] = std::max(startChan, endChan);
1932 if (cTimeAxis >= 0) {
1933 blc[cTimeAxis] = cTimeIdx;
1934 trc[cTimeAxis] = cTimeIdx;
1935 }
1936 blc[cNAxes] = cRowR;
1937 trc[cNAxes] = cRowR;
1938
1939 mbrec.haveSpectra = cGetSpectra;
1940 if (cGetSpectra) {
1941 int anynul;
1942
1943 for (int iPol = 0; iPol < nPol; iPol++) {
1944 blc[cStokesAxis] = iPol+1;
1945 trc[cStokesAxis] = iPol+1;
1946
1947 if (cALFA && cALFA_CIMA < 2) {
1948 // ALFA data: polarizations are stored in successive rows.
1949 blc[cStokesAxis] = 1;
1950 trc[cStokesAxis] = 1;
1951
1952 if (iPol) {
1953 if (++cRow > cNRow) {
1954 return -1;
1955 }
1956
1957 blc[cNAxes] = cRowR;
1958 trc[cNAxes] = cRowR;
1959 }
1960
1961 } else if (cData[DATA].nelem < 0) {
1962 // Variable dimension array; get axis lengths.
1963 int naxes = 5, status;
1964
1965 if ((status = readDim(DATA, cRowR, &naxes, cNAxis))) {
1966 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1967
1968 } else if ((status = (naxes != cNAxes))) {
1969 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "DATA array dimensions changed.");
1970 }
1971
1972 if (status) {
1973 delete [] blc;
1974 delete [] trc;
1975 delete [] inc;
1976 return 1;
1977 }
1978 }
1979
1980 if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAxes, cNAxis,
1981 blc, trc, inc, 0, mbrec.spectra[0] + iPol*nChan, &anynul,
1982 &cStatus)) {
1983 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1984 delete [] blc;
1985 delete [] trc;
1986 delete [] inc;
1987 return 1;
1988 }
1989
1990 if (endChan < startChan) {
1991 // Reverse the spectrum.
1992 float *iptr = mbrec.spectra[0] + iPol*nChan;
1993 float *jptr = iptr + nChan - 1;
1994 float *mid = iptr + nChan/2;
1995 while (iptr < mid) {
1996 float tmp = *iptr;
1997 *(iptr++) = *jptr;
1998 *(jptr--) = tmp;
1999 }
2000 }
2001
2002 if (cALFA) {
2003 // ALFA data, rescale the spectrum.
2004 float el, zd;
2005 readData(ELEVATIO, cRowR, &el);
2006 zd = 90.0f - el;
2007
2008 float factor = sALFAcal[iBeam][iPol] / alfaGain(zd);
2009
2010 if (cALFA_CIMA > 1) {
2011 // Rescale according to the number of unblanked accumulations.
2012 int colnum, naccum;
2013 findCol("STAT", &colnum);
2014 fits_read_col(cSDptr, TINT, colnum, cRowR, 10*(cTimeIdx-1)+2, 1, 0,
2015 &naccum, &anynul, &cStatus);
2016 factor *= cALFAacc / naccum;
2017 }
2018
2019 float *chan = mbrec.spectra[0] + iPol*nChan;
2020 float *chanN = chan + nChan;
2021 while (chan < chanN) {
2022 // Approximate conversion to Jy.
2023 *(chan++) *= factor;
2024 }
2025 }
2026
2027 if (mbrec.tsys[0][iPol] == 0.0) {
2028 // Compute Tsys as the average across the spectrum.
2029 float *chan = mbrec.spectra[0] + iPol*nChan;
2030 float *chanN = chan + nChan;
2031 float *tsys = mbrec.tsys[0] + iPol;
2032 while (chan < chanN) {
2033 *tsys += *(chan++);
2034 }
2035
2036 *tsys /= nChan;
2037 }
2038
2039 // Read data flags.
2040 if (cData[FLAGGED].colnum > 0) {
2041 if (fits_read_subset_byt(cSDptr, cData[FLAGGED].colnum, cNAxes,
2042 cNAxis, blc, trc, inc, 0, mbrec.flagged[0] + iPol*nChan, &anynul,
2043 &cStatus)) {
2044 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
2045 delete [] blc;
2046 delete [] trc;
2047 delete [] inc;
2048 return 1;
2049 }
2050
2051 if (endChan < startChan) {
2052 // Reverse the flag vector.
2053 unsigned char *iptr = mbrec.flagged[0] + iPol*nChan;
2054 unsigned char *jptr = iptr + nChan - 1;
2055 for (int ichan = 0; ichan < nChan/2; ichan++) {
2056 unsigned char tmp = *iptr;
2057 *(iptr++) = *jptr;
2058 *(jptr--) = tmp;
2059 }
2060 }
2061
2062 } else {
2063 // All channels are unflagged by default.
2064 unsigned char *iptr = mbrec.flagged[0] + iPol*nChan;
2065 for (int ichan = 0; ichan < nChan; ichan++) {
2066 *(iptr++) = 0;
2067 }
2068 }
2069 }
2070 }
2071
2072
2073 // Read cross-polarization data.
2074 if (cGetXPol) {
2075 int anynul;
2076 for (int j = 0; j < 2; j++) {
2077 mbrec.xcalfctr[0][j] = 0.0f;
2078 }
2079 if (cData[XCALFCTR].colnum > 0) {
2080 readData(XCALFCTR, cRowR, mbrec.xcalfctr);
2081 }
2082
2083 blc[0] = 1;
2084 trc[0] = 2;
2085 blc[1] = std::min(startChan, endChan);
2086 trc[1] = std::max(startChan, endChan);
2087 blc[2] = cRowR;
2088 trc[2] = cRowR;
2089
2090 int nAxis = 2;
2091 long nAxes[] = {2, nChan};
2092
2093 if (fits_read_subset_flt(cSDptr, cData[XPOLDATA].colnum, nAxis, nAxes,
2094 blc, trc, inc, 0, mbrec.xpol[0], &anynul, &cStatus)) {
2095 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
2096 delete [] blc;
2097 delete [] trc;
2098 delete [] inc;
2099 return 1;
2100 }
2101
2102 if (endChan < startChan) {
2103 // Invert the cross-polarization spectrum.
2104 float *iptr = mbrec.xpol[0];
2105 float *jptr = iptr + nChan - 2;
2106 for (int ichan = 0; ichan < nChan/2; ichan++) {
2107 float tmp = *iptr;
2108 *iptr = *jptr;
2109 *jptr = tmp;
2110
2111 tmp = *(iptr+1);
2112 *(iptr+1) = *(jptr+1);
2113 *(jptr+1) = tmp;
2114
2115 iptr += 2;
2116 jptr -= 2;
2117 }
2118 }
2119 }
2120
2121 delete [] blc;
2122 delete [] trc;
2123 delete [] inc;
2124
2125 if (cStatus) {
2126 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
2127 return 1;
2128 }
2129
2130 mbrec.extraSysCal = cExtraSysCal;
2131 readData(REFBEAM, cRowR, &mbrec.refBeam);
2132 readData(TCAL, cRowR, &mbrec.tcal[0]);
2133 readData(TCALTIME, cRowR, mbrec.tcalTime);
2134
2135 readData(AZIMUTH, cRowR, &mbrec.azimuth);
2136 readData(ELEVATIO, cRowR, &mbrec.elevation);
2137 readData(PARANGLE, cRowR, &mbrec.parAngle);
2138
2139 readData(FOCUSAXI, cRowR, &mbrec.focusAxi);
2140 readData(FOCUSTAN, cRowR, &mbrec.focusTan);
2141 readData(FOCUSROT, cRowR, &mbrec.focusRot);
2142
2143 readData(TAMBIENT, cRowR, &mbrec.temp);
2144 readData(PRESSURE, cRowR, &mbrec.pressure);
2145 readData(HUMIDITY, cRowR, &mbrec.humidity);
2146 readData(WINDSPEE, cRowR, &mbrec.windSpeed);
2147 readData(WINDDIRE, cRowR, &mbrec.windAz);
2148
2149 if (cALFA_BD) {
2150 // ALFA BDFITS stores zenith angle rather than elevation.
2151 mbrec.elevation = 90.0 - mbrec.elevation;
2152 }
2153
2154 mbrec.azimuth *= D2R;
2155 mbrec.elevation *= D2R;
2156 mbrec.parAngle *= D2R;
2157 mbrec.focusRot *= D2R;
2158 mbrec.windAz *= D2R;
2159
2160 // For GBT data, source velocity can be evaluated
2161 if ( cData[RVSYS].colnum > 0 && cData[VFRAME].colnum > 0 ) {
2162 double vframe;
2163 readData(VFRAME, cRowR, &vframe);
2164 double rvsys;
2165 readData(RVSYS, cRowR, &rvsys);
2166 //mbrec.srcVelocity = rvsys - vframe ;
2167 mbrec.srcVelocity = rvsys ;
2168 }
2169
2170 if (cStatus) {
2171 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
2172 return 1;
2173 }
2174
2175 return 0;
2176}
2177
2178//-------------------------------------------------------- GBTFITSreader::close
2179
2180// Close the SDFITS file.
2181
2182void GBTFITSreader::close()
2183{
2184 if (cSDptr) {
2185 int status = 0;
2186 fits_close_file(cSDptr, &status);
2187 cSDptr = 0x0;
2188
2189 if (cBeams) delete [] cBeams;
2190 if (cIFs) delete [] cIFs;
2191 if (cStartChan) delete [] cStartChan;
2192 if (cEndChan) delete [] cEndChan;
2193 if (cRefChan) delete [] cRefChan;
2194 if (cNRowT) delete [] cNRowT;
2195 if (cNRowCum) delete [] cNRowCum;
2196 if (hduId) delete [] hduId;
2197 }
2198}
2199
2200//------------------------------------------------------- GBTFITSreader::log
2201
2202// Log a message. If the current CFITSIO status value is non-zero, also log
2203// the corresponding error message and the CFITSIO message stack.
2204
2205void GBTFITSreader::log(LogOrigin origin, LogIO::Command cmd, const char *msg)
2206{
2207 LogIO os( origin ) ;
2208
2209 os << cmd ;
2210
2211 os << msg << endl ;
2212
2213 if (cStatus > 0) {
2214 fits_get_errstatus(cStatus, cMsg);
2215 os << cMsg << endl ;
2216
2217 while (fits_read_errmsg(cMsg)) {
2218 os << cMsg << endl ;
2219 }
2220 }
2221 os << LogIO::POST ;
2222}
2223
2224//----------------------------------------------------- GBTFITSreader::findData
2225
2226// Locate a data item in the SDFITS file.
2227
2228void GBTFITSreader::findData(
2229 int iData,
2230 char *name,
2231 int type)
2232{
2233 cData[iData].name = name;
2234 cData[iData].type = type;
2235
2236 int colnum;
2237 findCol(name, &colnum);
2238 cData[iData].colnum = colnum;
2239
2240 // Determine the number of data elements.
2241 if (colnum > 0) {
2242 int coltype;
2243 long nelem, width;
2244 fits_get_coltype(cSDptr, colnum, &coltype, &nelem, &width, &cStatus);
2245 fits_get_bcolparms(cSDptr, colnum, 0x0, cData[iData].units, 0x0, 0x0, 0x0,
2246 0x0, 0x0, 0x0, &cStatus);
2247
2248 // Look for a TDIMnnn keyword or column.
2249 char tdim[8];
2250 sprintf(tdim, "TDIM%d", colnum);
2251 findCol(tdim, &cData[iData].tdimcol);
2252
2253 if (coltype < 0) {
2254 // CFITSIO returns coltype < 0 for variable length arrays.
2255 cData[iData].coltype = -coltype;
2256 cData[iData].nelem = -nelem;
2257
2258 } else {
2259 cData[iData].coltype = coltype;
2260
2261 // Is there a TDIMnnn column?
2262 if (cData[iData].tdimcol > 0) {
2263 // Yes, dimensions of the fixed-length array could still vary.
2264 cData[iData].nelem = -nelem;
2265 } else {
2266 cData[iData].nelem = nelem;
2267 }
2268 }
2269
2270 } else if (colnum == 0) {
2271 // Keyword.
2272 cData[iData].coltype = 0;
2273 cData[iData].nelem = 1;
2274 cData[iData].tdimcol = -1;
2275 }
2276}
2277
2278//------------------------------------------------------ GBTFITSreader::findCol
2279
2280// Locate a parameter in the SDFITS file.
2281
2282void GBTFITSreader::findCol(
2283 char *name,
2284 int *colnum)
2285{
2286 *colnum = 0;
2287 int status = 0;
2288 fits_get_colnum(cSDptr, CASESEN, name, colnum, &status);
2289
2290 if (status) {
2291 // Not a real column - maybe it's virtual.
2292 char card[81];
2293
2294 status = 0;
2295 fits_read_card(cSDptr, name, card, &status);
2296 if (status) {
2297 // Not virtual either.
2298 *colnum = -1;
2299 }
2300
2301 // Clear error messages.
2302 fits_clear_errmsg();
2303 }
2304}
2305
2306//------------------------------------------------------ GBTFITSreader::readDim
2307
2308// Determine the dimensions of an array in the SDFITS file.
2309
2310int GBTFITSreader::readDim(
2311 int iData,
2312 long iRow,
2313 int *naxes,
2314 long naxis[])
2315{
2316 int colnum = cData[iData].colnum;
2317 if (colnum <= 0) {
2318 return 1;
2319 }
2320
2321 int maxdim = *naxes;
2322 if (cData[iData].tdimcol < 0) {
2323 // No TDIMnnn column for this array.
2324 if (cData[iData].nelem < 0) {
2325 // Variable length array; read the array descriptor.
2326 *naxes = 1;
2327 long dummy;
2328 if (fits_read_descript(cSDptr, colnum, iRow, naxis, &dummy, &cStatus)) {
2329 return 1;
2330 }
2331
2332 } else {
2333 // Read the repeat count from TFORMnnn.
2334 if (fits_read_tdim(cSDptr, colnum, maxdim, naxes, naxis, &cStatus)) {
2335 return 1;
2336 }
2337 }
2338
2339 } else {
2340 // Read the TDIMnnn value from the header or table.
2341 char tdim[8], tdimval[64];
2342 sprintf(tdim, "TDIM%d", colnum);
2343 readData(tdim, TSTRING, iRow, tdimval);
2344
2345 // fits_decode_tdim() checks that the TDIMnnn value is within the length
2346 // of the array in the specified column number but unfortunately doesn't
2347 // recognize variable-length arrays. Hence we must decode it here.
2348 char *tp = tdimval;
2349 if (*tp != '(') return 1;
2350
2351 tp++;
2352 *naxes = 0;
2353 for (size_t j = 1; j < strlen(tdimval); j++) {
2354 if (tdimval[j] == ',' || tdimval[j] == ')') {
2355 sscanf(tp, "%ld", naxis + (*naxes)++);
2356 if (tdimval[j] == ')') break;
2357 tp = tdimval + j + 1;
2358 }
2359 }
2360 }
2361
2362 return 0;
2363}
2364
2365//----------------------------------------------------- GBTFITSreader::readParm
2366
2367// Read a parameter value from the SDFITS file.
2368
2369int GBTFITSreader::readParm(
2370 char *name,
2371 int type,
2372 void *value)
2373{
2374 return readData(name, type, 1, value);
2375}
2376
2377//----------------------------------------------------- GBTFITSreader::readData
2378
2379// Read a data value from the SDFITS file.
2380
2381int GBTFITSreader::readData(
2382 char *name,
2383 int type,
2384 long iRow,
2385 void *value)
2386{
2387 int colnum;
2388 findCol(name, &colnum);
2389
2390 if (colnum > 0 && iRow > 0) {
2391 // Read the first value from the specified row of the table.
2392 int coltype;
2393 long nelem, width;
2394 fits_get_coltype(cSDptr, colnum, &coltype, &nelem, &width, &cStatus);
2395
2396 int anynul;
2397 if (type == TSTRING) {
2398 if (nelem) {
2399 fits_read_col(cSDptr, type, colnum, iRow, 1, 1, 0, &value, &anynul,
2400 &cStatus);
2401 } else {
2402 strcpy((char *)value, "");
2403 }
2404
2405 } else {
2406 if (nelem) {
2407 fits_read_col(cSDptr, type, colnum, iRow, 1, 1, 0, value, &anynul,
2408 &cStatus);
2409 } else {
2410 if (type == TSHORT) {
2411 *((short *)value) = 0;
2412 } else if (type == TINT) {
2413 *((int *)value) = 0;
2414 } else if (type == TFLOAT) {
2415 *((float *)value) = 0.0f;
2416 } else if (type == TDOUBLE) {
2417 *((double *)value) = 0.0;
2418 }
2419 }
2420 }
2421
2422 } else if (colnum == 0) {
2423 // Read keyword value.
2424 fits_read_key(cSDptr, type, name, value, 0, &cStatus);
2425
2426 } else {
2427 // Not present.
2428 if (type == TSTRING) {
2429 strcpy((char *)value, "");
2430 } else if (type == TSHORT) {
2431 *((short *)value) = 0;
2432 } else if (type == TINT) {
2433 *((int *)value) = 0;
2434 } else if (type == TFLOAT) {
2435 *((float *)value) = 0.0f;
2436 } else if (type == TDOUBLE) {
2437 *((double *)value) = 0.0;
2438 }
2439 }
2440
2441 return colnum < 0;
2442}
2443
2444//----------------------------------------------------- GBTFITSreader::readData
2445
2446// Read data from the SDFITS file.
2447
2448int GBTFITSreader::readData(
2449 int iData,
2450 long iRow,
2451 void *value)
2452{
2453 int type = cData[iData].type;
2454 int colnum = cData[iData].colnum;
2455
2456 if (colnum > 0 && iRow > 0) {
2457 // Read the required number of values from the specified row of the table.
2458 long nelem = cData[iData].nelem;
2459 int anynul;
2460 if (type == TSTRING) {
2461 if (nelem) {
2462 fits_read_col(cSDptr, type, colnum, iRow, 1, 1, 0, &value, &anynul,
2463 &cStatus);
2464 } else {
2465 strcpy((char *)value, "");
2466 }
2467
2468 } else {
2469 if (nelem) {
2470 fits_read_col(cSDptr, type, colnum, iRow, 1, abs(nelem), 0, value,
2471 &anynul, &cStatus);
2472 } else {
2473 if (type == TSHORT) {
2474 *((short *)value) = 0;
2475 } else if (type == TINT) {
2476 *((int *)value) = 0;
2477 } else if (type == TFLOAT) {
2478 *((float *)value) = 0.0f;
2479 } else if (type == TDOUBLE) {
2480 *((double *)value) = 0.0;
2481 }
2482 }
2483 }
2484
2485 } else if (colnum == 0) {
2486 // Read keyword value.
2487 char *name = cData[iData].name;
2488 fits_read_key(cSDptr, type, name, value, 0, &cStatus);
2489
2490 } else {
2491 // Not present.
2492 if (type == TSTRING) {
2493 strcpy((char *)value, "");
2494 } else if (type == TSHORT) {
2495 *((short *)value) = 0;
2496 } else if (type == TINT) {
2497 *((int *)value) = 0;
2498 } else if (type == TFLOAT) {
2499 *((float *)value) = 0.0f;
2500 } else if (type == TDOUBLE) {
2501 *((double *)value) = 0.0;
2502 }
2503 }
2504
2505 return colnum < 0;
2506}
2507
2508//------------------------------------------------------ GBTFITSreader::readCol
2509
2510// Read a scalar column from the SDFITS file.
2511
2512int GBTFITSreader::readCol(
2513 int iData,
2514 void *value)
2515{
2516 int type = cData[iData].type;
2517
2518 if (cData[iData].colnum > 0) {
2519 // Table column.
2520 int anynul;
2521 int offset = 0 ;
2522 for ( int ihdu = 0 ; ihdu < numHDU ; ihdu++ ) {
2523 fits_movabs_hdu(cSDptr, hduId[ihdu], NULL, &cStatus) ;
2524 if ( type == TSHORT ) {
2525 fits_read_col(cSDptr, type, cData[iData].colnum, 1, 1, cNRowT[ihdu], 0,
2526 &((short *)value)[offset], &anynul, &cStatus);
2527 }
2528 else if ( type == TINT ) {
2529 fits_read_col(cSDptr, type, cData[iData].colnum, 1, 1, cNRowT[ihdu], 0,
2530 &((int *)value)[offset], &anynul, &cStatus);
2531 }
2532 else if ( type == TFLOAT ) {
2533 fits_read_col(cSDptr, type, cData[iData].colnum, 1, 1, cNRowT[ihdu], 0,
2534 &((float *)value)[offset], &anynul, &cStatus);
2535 }
2536 else if ( type == TDOUBLE ) {
2537 fits_read_col(cSDptr, type, cData[iData].colnum, 1, 1, cNRowT[ihdu], 0,
2538 &((double *)value)[offset], &anynul, &cStatus);
2539 }
2540 offset = cNRowCum[ihdu] ;
2541 }
2542 } else {
2543 // Header keyword.
2544 readData(iData, 0, value);
2545 for (int irow = 1; irow < cNRow; irow++) {
2546 if (type == TSHORT) {
2547 ((short *)value)[irow] = *((short *)value);
2548 } else if (type == TINT) {
2549 ((int *)value)[irow] = *((int *)value);
2550 } else if (type == TFLOAT) {
2551 ((float *)value)[irow] = *((float *)value);
2552 } else if (type == TDOUBLE) {
2553 ((double *)value)[irow] = *((double *)value);
2554 }
2555 }
2556 }
2557
2558 return cData[iData].colnum < 0;
2559}
2560
2561//----------------------------------------------------- GBTFITSreader::readTime
2562
2563// Read the time from the SDFITS file.
2564
2565int GBTFITSreader::readTime(
2566 long iRow,
2567 int iPix,
2568 char *datobs,
2569 double &utc)
2570{
2571 readData(DATE_OBS, iRow, datobs);
2572 if (cData[TIME].colnum >= 0) {
2573 readData(TIME, iRow, &utc);
2574 } else if (cGBT) {
2575 Int yy, mm ;
2576 Double dd, hour, min, sec ;
2577 sscanf( datobs, "%d-%d-%lfT%lf:%lf:%lf", &yy, &mm, &dd, &hour, &min, &sec ) ;
2578 dd = dd + ( hour * 3600.0 + min * 60.0 + sec ) / 86400.0 ;
2579 MVTime mvt( yy, mm, dd ) ;
2580 dd = mvt.day() ;
2581 utc = fmod( dd, 1.0 ) * 86400.0 ;
2582 } else if (cNAxisTime > 1) {
2583 double timeDelt, timeRefPix, timeRefVal;
2584 readData(TimeRefVal, iRow, &timeRefVal);
2585 readData(TimeDelt, iRow, &timeDelt);
2586 readData(TimeRefPix, iRow, &timeRefPix);
2587 utc = timeRefVal + (iPix - timeRefPix) * timeDelt;
2588 }
2589
2590 if (cALFA_BD) utc *= 3600.0;
2591
2592 // Check DATE-OBS format.
2593 if (datobs[2] == '/') {
2594 // Translate an old-format DATE-OBS.
2595 datobs[9] = datobs[1];
2596 datobs[8] = datobs[0];
2597 datobs[2] = datobs[6];
2598 datobs[5] = datobs[3];
2599 datobs[3] = datobs[7];
2600 datobs[6] = datobs[4];
2601 datobs[7] = '-';
2602 datobs[4] = '-';
2603 datobs[1] = '9';
2604 datobs[0] = '1';
2605
2606 } else if (datobs[10] == 'T' && cData[TIME].colnum < 0) {
2607 // Dig UTC out of a new-format DATE-OBS.
2608 int hh, mm;
2609 float ss;
2610 sscanf(datobs+11, "%d:%d:%f", &hh, &mm, &ss);
2611 utc = (hh*60 + mm)*60 + ss;
2612 }
2613
2614 datobs[10] = '\0';
2615
2616 return 0;
2617}
2618
2619//------------------------------------------------------ GBTFITSreader::alfaCal
2620
2621// Process ALFA calibration data.
2622
2623int GBTFITSreader::alfaCal(
2624 short iBeam,
2625 short iIF,
2626 short iPol)
2627{
2628 const string methodName = "alfaCal()" ;
2629
2630 int calOn;
2631 char chars[32];
2632 if (cALFA_BD) {
2633 readData("OBS_NAME", TSTRING, cRow, chars);
2634 } else {
2635 readData("SCANTYPE", TSTRING, cRow, chars);
2636 }
2637
2638 if (strcmp(chars, "ON") == 0) {
2639 calOn = 1;
2640 } else if (strcmp(chars, "OFF") == 0) {
2641 calOn = 0;
2642 } else {
2643 return 1;
2644 }
2645
2646 // Read cal data.
2647 long *blc = new long[cNAxes+1];
2648 long *trc = new long[cNAxes+1];
2649 long *inc = new long[cNAxes+1];
2650 for (int iaxis = 0; iaxis <= cNAxes; iaxis++) {
2651 blc[iaxis] = 1;
2652 trc[iaxis] = 1;
2653 inc[iaxis] = 1;
2654 }
2655
2656 // User channel selection.
2657 int startChan = cStartChan[iIF];
2658 int endChan = cEndChan[iIF];
2659
2660 blc[cFreqAxis] = std::min(startChan, endChan);
2661 trc[cFreqAxis] = std::max(startChan, endChan);
2662 if (cALFA_CIMA > 1) {
2663 // CIMAFITS 2.x has a legitimate STOKES axis...
2664 blc[cStokesAxis] = iPol+1;
2665 trc[cStokesAxis] = iPol+1;
2666 } else {
2667 // ...older ALFA data does not.
2668 blc[cStokesAxis] = 1;
2669 trc[cStokesAxis] = 1;
2670 }
2671 if (cTimeAxis >= 0) {
2672 blc[cTimeAxis] = cTimeIdx;
2673 trc[cTimeAxis] = cTimeIdx;
2674 }
2675 blc[cNAxes] = cRow;
2676 trc[cNAxes] = cRow;
2677
2678 float spectrum[endChan];
2679 int anynul;
2680 if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAxes, cNAxis,
2681 blc, trc, inc, 0, spectrum, &anynul, &cStatus)) {
2682 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
2683 delete [] blc;
2684 delete [] trc;
2685 delete [] inc;
2686 return 1;
2687 }
2688
2689 // Factor to rescale according to the number of unblanked accumulations.
2690 float factor = 1.0f;
2691 if (cALFA_CIMA > 1) {
2692 int colnum, naccum;
2693 findCol("STAT", &colnum);
2694 fits_read_col(cSDptr, TINT, colnum, cRow, 2, 1, 0, &naccum, &anynul,
2695 &cStatus);
2696 factor = cALFAacc / naccum;
2697 }
2698
2699 // Average the spectrum.
2700 float mean = 1e9f;
2701 for (int k = 0; k < 2; k++) {
2702 float discrim = 2.0f * mean;
2703
2704 int nChan = 0;
2705 float sum = 0.0f;
2706
2707 float *chanN = spectrum + abs(endChan - startChan) + 1;
2708 for (float *chan = spectrum; chan < chanN; chan++) {
2709 // Simple discriminant that eliminates strong radar interference.
2710 if (*chan < discrim) {
2711 nChan++;
2712 sum += *chan * factor;
2713 }
2714 }
2715
2716 mean = sum / nChan;
2717 }
2718
2719 if (calOn) {
2720 sALFAcalOn[iBeam][iPol] *= sALFAcalNon[iBeam][iPol];
2721 sALFAcalOn[iBeam][iPol] += mean;
2722 sALFAcalOn[iBeam][iPol] /= ++sALFAcalNon[iBeam][iPol];
2723 } else {
2724 sALFAcalOff[iBeam][iPol] *= sALFAcalNoff[iBeam][iPol];
2725 sALFAcalOff[iBeam][iPol] += mean;
2726 sALFAcalOff[iBeam][iPol] /= ++sALFAcalNoff[iBeam][iPol];
2727 }
2728
2729 if (sALFAcalNon[iBeam][iPol] && sALFAcalNoff[iBeam][iPol]) {
2730 // Tcal should come from the TCAL table, it varies weakly with beam,
2731 // polarization, and frequency. However, TCAL is not written properly.
2732 float Tcal = 12.0f;
2733 sALFAcal[iBeam][iPol] = Tcal / (sALFAcalOn[iBeam][iPol] -
2734 sALFAcalOff[iBeam][iPol]);
2735
2736 // Scale from K to Jy; the gain also varies weakly with beam,
2737 // polarization, frequency, and zenith angle.
2738 float fluxCal = 10.0f;
2739 sALFAcal[iBeam][iPol] /= fluxCal;
2740 }
2741
2742 return 0;
2743}
2744
2745//----------------------------------------------------- GBTFITSreader::alfaGain
2746
2747// ALFA gain factor.
2748
2749float GBTFITSreader::alfaGain(
2750 float zd)
2751{
2752 // Gain vs zenith distance table from Robert Minchin, 2008/12/08.
2753 const int nZD = 37;
2754 const float zdLim[] = {1.5f, 19.5f};
2755 const float zdInc = (nZD - 1) / (zdLim[1] - zdLim[0]);
2756 float zdGain[] = { 1.00723708,
2757 1.16644573, 1.15003645, 1.07117307, 1.02532673,
2758 1.01788402, 1.01369524, 1.00000000, 0.989855111,
2759 0.990888834, 0.993996620, 0.989964068, 0.982213855,
2760 0.978662670, 0.979349494, 0.978478372, 0.974631131,
2761 0.972126007, 0.972835243, 0.972742677, 0.968671739,
2762 0.963891327, 0.963452935, 0.966831207, 0.969585896,
2763 0.970700860, 0.972644389, 0.973754644, 0.967344403,
2764 0.952168941, 0.937160134, 0.927843094, 0.914048433,
2765 0.886700928, 0.864701211, 0.869126320, 0.854309499};
2766
2767 float gain;
2768 // Do table lookup by linear interpolation.
2769 float lambda = zdInc * (zd - zdLim[0]);
2770 int j = int(lambda);
2771 if (j < 0) {
2772 gain = zdGain[0];
2773 } else if (j >= nZD-1) {
2774 gain = zdGain[nZD-1];
2775 } else {
2776 gain = zdGain[j] + (lambda - j) * (zdGain[j+1] - zdGain[j]);
2777 }
2778
2779 return gain;
2780}
2781
2782//----------------------------------------------------- GBTFITSreader::cRowRelative
2783
2784// get relative row index in appropriate binary table
2785
2786long GBTFITSreader::cRowRelative( long iRow )
2787{
2788 // iRow is 1-relative
2789 const string methodName = "cRowRelative" ;
2790 if ( numHDU == 0 ) {
2791 log( LogOrigin( className, methodName, WHERE ), LogIO::SEVERE,
2792 "You should call open() first" ) ;
2793 return 0L ;
2794 }
2795 else if ( numHDU == 1 ) {
2796 if ( iRow > cNRow )
2797 return 0L ;
2798 else
2799 return iRow ;
2800 }
2801
2802 long jRow = 0L ;
2803 sprintf( cMsg, "iRow = %ld", iRow );
2804 log(LogOrigin( className, methodName, WHERE ), LogIO::DEBUGGING, cMsg) ;
2805 if ( iRow <= cNRowCum[0] ) {
2806 jRow = iRow ;
2807 sprintf( cMsg, "Moving %d(%d)", 0, hduId[0] ) ;
2808 log(LogOrigin( className, methodName, WHERE ), LogIO::DEBUGGING, cMsg) ;
2809 if ( fits_movabs_hdu( cSDptr, hduId[0], NULL, &cStatus ) ) {
2810 log( LogOrigin( className, methodName, WHERE ), LogIO::SEVERE,
2811 "Failed to move HDU" ) ;
2812 return 0L ;
2813 }
2814 }
2815 else {
2816 for ( int ihdu = 1 ; ihdu < numHDU ; ihdu++ ) {
2817 if ( iRow <= cNRowCum[ihdu] ) {
2818 jRow = iRow - cNRowCum[ihdu-1] ;
2819 sprintf( cMsg, "Moving %d(%d)", ihdu, hduId[ihdu] ) ;
2820 log(LogOrigin( className, methodName, WHERE ), LogIO::DEBUGGING, cMsg) ;
2821 if ( fits_movabs_hdu( cSDptr, hduId[ihdu], NULL, &cStatus ) ) {
2822 log( LogOrigin( className, methodName, WHERE ), LogIO::SEVERE,
2823 "Failed to move HDU" ) ;
2824 return 0L ;
2825 }
2826 break ;
2827 }
2828 }
2829 }
2830 sprintf( cMsg, "jRow = %ld", jRow );
2831 log(LogOrigin( className, methodName, WHERE ), LogIO::DEBUGGING, cMsg) ;
2832
2833 return jRow ;
2834}
2835
2836//----------------------------------------------------- GBTFITSreader::nIFPerHDU
2837
2838// get number of IFs per HDU
2839
2840Vector<Int> GBTFITSreader::nIFPerHDU()
2841{
2842 const string methodName = "nIFPerHDU()" ;
2843 //cNumKeys.resize( numHDU ) ;
2844 cNumKeys.resize( 0 ) ;
2845 cRefValKey.resize( 0 ) ;
2846 cRefVal.resize( 0 ) ;
2847 cIncVal.resize( 0 ) ;
2848 cIFNO.resize( 0 ) ;
2849 Vector<Int> nIF_per_HDU( numHDU, 0 ) ;
2850 for ( int ihdu = 0 ; ihdu < numHDU ; ihdu++ ) {
2851 Vector<String> sampler( 0 ) ;
2852 Vector<Double> rfreq( 0 ) ;
2853 Vector<String> samprf( 0 ) ;
2854 if ( fits_movabs_hdu( cSDptr, hduId[ihdu], NULL, &cStatus ) ) {
2855 log( LogOrigin( className, methodName, WHERE ), LogIO::SEVERE,
2856 "Failed to move HDU" ) ;
2857 nIF_per_HDU = -1 ;
2858 return nIF_per_HDU ;
2859 }
2860 int startpos = (int) cRefValKey.size() ;
2861 int numkeys = 0 ;
2862 for ( int irow = 0 ; irow < cNRowT[ihdu] ; irow++ ) {
2863 char samp[24] ;
2864 double rf ;
2865 readData( SAMPLER, irow+1, samp ) ;
2866 readData( RESTFRQ, irow+1, &rf ) ;
2867 char sarf[48] ;
2868 sprintf( sarf, "%d/%s/%lf", ihdu, samp, rf ) ;
2869 String sr = String( sarf ) ;
2870 if ( allNE( sr, samprf ) ) {
2871 int oldsize = samprf.size() ;
2872 samprf.resize( oldsize+1, True ) ;
2873 samprf[oldsize] = sr ;
2874 sprintf( cMsg, "added samprf = %s", sarf ) ;
2875 log( LogOrigin( className, methodName, WHERE ), LogIO::DEBUGGING, cMsg ) ;
2876 Double refval = getRefValLSR( irow+1 ) ;
2877 double increm ;
2878 readData(FqDelt, irow+1, &increm);
2879 cRefValKey.push_back( sr ) ;
2880 cRefVal.push_back( (double)refval ) ;
2881 cIncVal.push_back( increm ) ;
2882 numkeys++ ;
2883 }
2884 }
2885 cNumKeys.push_back( numkeys ) ;
2886 nIF_per_HDU[ihdu] = mergeIFPerHDU( startpos, numkeys ) ;
2887 sprintf( cMsg, "number of IF for HDU %d = %d",
2888 hduId[ihdu], nIF_per_HDU[ihdu] ) ;
2889 log( LogOrigin( className, methodName, WHERE ), LogIO::NORMAL, cMsg ) ;
2890 }
2891
2892 return nIF_per_HDU ;
2893}
2894
2895Double GBTFITSreader::getRefValLSR( int irow )
2896{
2897 const string methodName = "getRefValLSR()" ;
2898
2899 char keyw[9], ctype[9];
2900 sprintf(keyw, "CTYPE%ld", cFreqAxis+1);
2901 readParm(keyw, TSTRING, ctype);
2902 double fqRefVal ;
2903 readData(FqRefVal, irow, &fqRefVal);
2904 if ( strstr( ctype, "-LSR" ) == NULL && strstr( ctype, "LSRK" ) == NULL ) {
2905 // need conversion
2906 Int year, month ;
2907 Double dd, hour, min, sec ;
2908 char dateObs[32] ;
2909 readData(DATE_OBS, irow, dateObs);
2910 sscanf( dateObs, "%4d-%2d-%2lfT%lf:%lf:%lf", &year, &month, &dd, &hour, &min, &sec ) ;
2911 dd = dd + ( hour * 3600.0 + min * 60.0 + sec ) / 86400.0 ;
2912 dd = MVTime(year, month, dd).day() ;
2913 double duration ;
2914 readData(DURATION, irow, &duration) ;
2915 dd += 0.5 * duration / 86400.0 ;
2916 MEpoch me( Quantity( dd, "d" ), MEpoch::UTC ) ;
2917 double antPos[3] ;
2918 readParm( "SITELONG", TDOUBLE, antPos ) ;
2919 readParm( "SITELAT", TDOUBLE, antPos+1 ) ;
2920 readParm( "SITEELEV", TDOUBLE, antPos+2 ) ;
2921 Vector<Double> pos( 2 ) ;
2922 pos[0] = Double( antPos[0] ) ;
2923 pos[1] = Double( antPos[1] ) ;
2924 MPosition mp( MVPosition( Quantity( Double(antPos[2]), "m" ),
2925 Quantum< Vector<Double> >( pos, "deg" ) ),
2926 MPosition::WGS84 ) ;
2927 double ra, dec ;
2928 readData(RA, irow, &ra);
2929 readData(DEC, irow, &dec);
2930 MDirection md( Quantity( Double(ra), "deg" ), Quantity( Double(dec), "deg" ), MDirection::J2000 ) ;
2931 MeasFrame mf( me, mp, md ) ;
2932 MFrequency::Types fromType = MFrequency::DEFAULT ;
2933 if ( strstr( ctype, "-OBS" ) != NULL ) {
2934 fromType = MFrequency::TOPO ;
2935 }
2936 else if ( strstr( ctype, "-HEL" ) != NULL ) {
2937 fromType = MFrequency::BARY ;
2938 }
2939 else if ( strstr( ctype, "-GEO" ) != NULL ) {
2940 fromType = MFrequency::GEO ;
2941 }
2942 else if ( strstr( ctype, "-GAL" ) != NULL ) {
2943 fromType = MFrequency::GALACTO ;
2944 }
2945 else if ( strstr( ctype, "REST" ) != NULL ) {
2946 fromType = MFrequency::REST ;
2947 }
2948 MFrequency::Convert tolsr( fromType, MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
2949 Double world = tolsr( Double( fqRefVal ) ).get( "Hz" ).getValue() ;
2950 char sampler[10];
2951 readData(SAMPLER,irow,sampler);
2952 int hduId = 0 ;
2953 fits_get_hdu_num( cSDptr, &hduId ) ;
2954 sprintf( cMsg, "HDU %d Row %d (%s): before conversion %.7f, after conversion %.7f", hduId, irow, sampler, fqRefVal, world ) ;
2955 log( LogOrigin( className, methodName, WHERE ), LogIO::DEBUGGING, cMsg ) ;
2956 fqRefVal = (double)world ;
2957 }
2958
2959 return Double(fqRefVal) ;
2960}
2961
2962Double GBTFITSreader::getRefValLSR( long iRow, int idx, char *sampler, double restfreq )
2963{
2964 const string methodName = "getRefValLSR()" ;
2965
2966 // get ihdu
2967// int hduid = 0 ;
2968// while ( iRow > cNRowCum[hduid] && hduid < numHDU) {
2969// hduid++ ;
2970// }
2971
2972 // create key
2973// char sarf[48] ;
2974// sprintf( sarf, "%d/%s/%lf", hduid, sampler, restfreq ) ;
2975// String sr = String( sarf ) ;
2976
2977 sprintf( cMsg, "Row %ld (%s): reference frequency %lf", iRow, sampler, cRefVal[idx] ) ;
2978 log( LogOrigin( className, methodName, WHERE ), LogIO::DEBUGGING, cMsg ) ;
2979
2980 if ( idx == (int)cRefValKey.size() ) {
2981 log( LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Failed to get reference frequency value" ) ;
2982 return 0.0 ;
2983 }
2984
2985 return cRefVal[idx] ;
2986}
2987
2988Int GBTFITSreader::mergeIFPerHDU( int start, int n )
2989{
2990 Int if0 = 0 ;
2991 if ( cIFNO.size() > 0 ) {
2992 if0 = max( cIFNO ) + 1 ;
2993 }
2994 cIFNO.resize( start+n, True ) ;
2995 Vector< Vector<Double> > lIncVal( 0 ) ;
2996 for ( int i = start ; i < start+n ; i++ ) {
2997 Vector<Double> v( 2 ) ;
2998 v[0] = cRefVal[i] ;
2999 v[1] = cIncVal[i] ;
3000 int idx = -1 ;
3001 for ( int j = 0 ; j < (int)lIncVal.size() ; j++ ) {
3002 if ( allEQ( v, lIncVal[j] ) ) {
3003 idx = j ;
3004 break ;
3005 }
3006 }
3007 if ( idx == -1 ) {
3008 int oldsize = lIncVal.size() ;
3009 lIncVal.resize( oldsize+1, True ) ;
3010 lIncVal[oldsize] = v ;
3011 cIFNO[i] = if0 + oldsize ;
3012 }
3013 else {
3014 cIFNO[i] = if0 + idx ;
3015 }
3016 }
3017 return max(cIFNO) - if0 + 1 ;
3018}
Note: See TracBrowser for help on using the repository browser.