Context Navigation

GBTFITSreader.cc @ 3067

Visit:

Last change on this file since 3067 was 3067, checked in by Takeshi Nakazato, 8 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: Yes/No?

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

Make PKS classes warning free.

File size: 85.6 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
61	class 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.
74	const double PI = 3.141592653589793238462643;
75
76	// Factor to convert radians to degrees.
77	const double D2R = PI / 180.0;
78
79	// Class name
80	const string className = "GBTFITSreader" ;
81
82	//---------------------------------------------------- GBTFITSreader::(statics)
83
84	int GBTFITSreader::sInit = 1;
85	int GBTFITSreader::sReset = 0;
86	int (GBTFITSreader::sALFAcalNon)[2] = (int ()[2])(new float[16]);
87	int (GBTFITSreader::sALFAcalNoff)[2] = (int ()[2])(new float[16]);
88	float (GBTFITSreader::sALFAcalOn)[2] = (float ()[2])(new float[16]);
89	float (GBTFITSreader::sALFAcalOff)[2] = (float ()[2])(new float[16]);
90	float (GBTFITSreader::sALFAcal)[2] = (float ()[2])(new float[16]);
91
92	//------------------------------------------------- GBTFITSreader::GBTFITSreader
93
94	GBTFITSreader::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
121	GBTFITSreader::~GBTFITSreader()
122	{
123	close();
124
125	delete [] cData;
126	}
127
128	//--------------------------------------------------------- GBTFITSreader::open
129
130	// Open an SDFITS file for reading.
131
132	int 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 = (char )"SINGLE DISH" ;
173	if (fits_movnam_hdu(cSDptr, BINARY_TBL, extName, 0, &cStatus)) {
174	// No SDFITS table, look for BDFITS or CIMAFITS.
175	cStatus = 0;
176	if (fits_movnam_hdu(cSDptr, BINARY_TBL, (char *)"BDFITS", 0, &cStatus) == 0) {
177	cALFA_BD = 1;
178	extName = (char *)"BDFITS" ;
179
180	} else {
181	cStatus = 0;
182	if (fits_movnam_hdu(cSDptr, BINARY_TBL, (char *)"CIMAFITS", 0, &cStatus) == 0) {
183	cALFA_CIMA = 1;
184	extName = (char *)"CIMAFITS" ;
185
186	// Check for later versions of CIMAFITS.
187	float version;
188	readParm((char )(char )"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((char )(char )"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, (char )(char )"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((char *)"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, (char )(char )"FLAGGED", TBYTE);
328
329	} else {
330	// DATA column not present, check for a DATAXED keyword.
331	findData(DATAXED, (char )(char )"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((char )(char )"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] = {(char )"CTYPE1", (char )"CTYPE2", (char )"CTYPE3", (char )"CTYPE4", (char )"CTYPE5"};
351	char CRPIX[5] = {(char )"CRPIX1", (char )"CRPIX2", (char )"CRPIX3", (char )"CRPIX4", (char )"CRPIX5"};
352	char CRVAL[5] = {(char )"CRVAL1", (char )"CRVAL2", (char )"CRVAL3", (char )"CRVAL4", (char )"CRVAL5"};
353	char CDELT[5] = {(char )"CDELT1", (char )"CDELT2", (char )"CDELT3", (char )"CDELT4", (char )"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 = (char *)"CRVAL2A";
432
433	cDecAxis = 3;
434	decCRVAL = (char *)"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, (char )(char )"SCAN", TINT); // Shared.
447	findData(CYCLE, (char )(char )"CYCLE", TINT); // Additional.
448	findData(DATE_OBS, (char *)"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, (char )(char )"TIME", TDOUBLE); // Core.
464	}
465
466	findData(EXPOSURE, (char )(char )"EXPOSURE", TFLOAT); // Core.
467	findData(OBJECT, (char )(char )"OBJECT", TSTRING); // Core.
468	findData(OBJ_RA, (char *)"OBJ-RA", TDOUBLE); // Additional.
469	findData(OBJ_DEC, (char *)"OBJ-DEC", TDOUBLE); // Additional.
470	findData(RESTFRQ, (char )(char )"RESTFRQ", TDOUBLE); // Additional.
471	findData(OBSMODE, (char )(char )"OBSMODE", TSTRING); // Shared.
472
473	findData(BEAM, (char )(char )"BEAM", TSHORT); // Additional.
474	findData(IF, (char )(char )"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, (char )(char )"SCANRATE", TFLOAT); // Additional.
481
482	findData(TSYS, (char )(char )"TSYS", TFLOAT); // Core.
483	findData(CALFCTR, (char )(char )"CALFCTR", TFLOAT); // Additional.
484	findData(XCALFCTR, (char )(char )"XCALFCTR", TFLOAT); // Additional.
485	findData(BASELIN, (char )(char )"BASELIN", TFLOAT); // Additional.
486	findData(BASESUB, (char )(char )"BASESUB", TFLOAT); // Additional.
487	findData(XPOLDATA, (char )(char )"XPOLDATA", TFLOAT); // Additional.
488
489	findData(REFBEAM, (char )(char )"REFBEAM", TSHORT); // Additional.
490	findData(TCAL, (char )(char )"TCAL", TFLOAT); // Shared.
491	findData(TCALTIME, (char )(char )"TCALTIME", TSTRING); // Additional.
492	findData(AZIMUTH, (char )(char )"AZIMUTH", TFLOAT); // Shared.
493	findData(ELEVATIO, (char )(char )"ELEVATIO", TFLOAT); // Shared.
494	findData(PARANGLE, (char )(char )"PARANGLE", TFLOAT); // Additional.
495	findData(FOCUSAXI, (char )(char )"FOCUSAXI", TFLOAT); // Additional.
496	findData(FOCUSTAN, (char )(char )"FOCUSTAN", TFLOAT); // Additional.
497	findData(FOCUSROT, (char )(char )"FOCUSROT", TFLOAT); // Additional.
498	findData(TAMBIENT, (char )(char )"TAMBIENT", TFLOAT); // Shared.
499	findData(PRESSURE, (char )(char )"PRESSURE", TFLOAT); // Shared.
500	findData(HUMIDITY, (char )(char )"HUMIDITY", TFLOAT); // Shared.
501	findData(WINDSPEE, (char )(char )"WINDSPEE", TFLOAT); // Shared.
502	findData(WINDDIRE, (char )(char )"WINDDIRE", TFLOAT); // Shared.
503
504	findData(STOKES, polCRVAL, TINT);
505	findData(SIG, (char )(char )"SIG", TSTRING);
506	findData(CAL, (char )(char )"CAL", TSTRING);
507
508	findData(RVSYS, (char )(char )"RVSYS", TDOUBLE);
509	findData(VFRAME, (char )(char )"VFRAME", TDOUBLE);
510	findData(VELDEF, (char )(char )"VELDEF", TSTRING);
511
512	findData(TIMESTAMP, (char )(char )"TIMESTAMP", TSTRING);
513	findData(DURATION, (char )(char )"DURATION", TDOUBLE);
514	findData(SAMPLER, (char )(char )"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, (char *)"SCAN_ID", TINT);
530	if (cALFA_CIMA > 1) {
531	// Note that RECNUM increases by cNAxisTime per row.
532	findData(CYCLE, (char )(char )"RECNUM", TINT);
533	} else {
534	findData(CYCLE, (char )(char )"SUBSCAN", TINT);
535	}
536	} else if (cALFA_BD) {
537	findData(SCAN, (char *)"SCAN_NUMBER", TINT);
538	findData(CYCLE, (char *)"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, (char *)"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, (char *)"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, (char )(char )"IFN", TSHORT);
576	} else {
577	findData(IF, (char )(char )"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((char *)"SYSCAL_TSYS", &colnum);
585	if (colnum > 0) {
586	// This contains the vector Tsys.
587	findData(TSYS, (char *)"SYSCAL_TSYS", TFLOAT);
588	}
589
590	// XPOLDATA?
591
592	if (cData[SCANRATE].colnum < 0) {
593	findData(SCANRATE, (char *)"FIELD_POINTING_DIR_RATE", TFLOAT);
594	}
595
596	if (cData[RESTFRQ].colnum < 0) {
597	findData(RESTFRQ, (char )(char )"RESTFREQ", TDOUBLE);
598	if (cData[RESTFRQ].colnum < 0) {
599	findData(RESTFRQ, (char *)"SPECTRAL_WINDOW_REST_FREQUENCY", TDOUBLE);
600	}
601	}
602
603	if (cData[OBJ_RA].colnum < 0) {
604	findData(OBJ_RA, (char *)"SOURCE_DIRECTION", TDOUBLE);
605	}
606	if (cData[OBJ_DEC].colnum < 0) {
607	findData(OBJ_DEC, (char *)"SOURCE_DIRECTION", TDOUBLE);
608	}
609
610	// REFBEAM?
611
612	if (cData[TCAL].colnum < 0) {
613	findData(TCAL, (char *)"SYSCAL_TCAL", TFLOAT);
614	} else if (cALFA_BD) {
615	// ALFA BDFITS has a different TCAL with 64 elements - kill it!
616	findData(TCAL, (char *)"NO NO NO", TFLOAT);
617	}
618
619	if (cALFA_BD) {
620	// ALFA BDFITS.
621	findData(AZIMUTH, (char )(char )"CRVAL2B", TFLOAT);
622	findData(ELEVATIO, (char )(char )"CRVAL3B", TFLOAT);
623	}
624
625	if (cALFA) {
626	// ALFA data.
627	findData(PARANGLE, (char *)"PARA_ANG", TFLOAT);
628	}
629
630	if (cData[TAMBIENT].colnum < 0) {
631	findData(TAMBIENT, (char *)"WEATHER_TEMPERATURE", TFLOAT);
632	}
633
634	if (cData[PRESSURE].colnum < 0) {
635	findData(PRESSURE, (char *)"WEATHER_PRESSURE", TFLOAT);
636	}
637
638	if (cData[HUMIDITY].colnum < 0) {
639	findData(HUMIDITY, (char *)"WEATHER_REL_HUMIDITY", TFLOAT);
640	}
641
642	if (cData[WINDSPEE].colnum < 0) {
643	findData(WINDSPEE, (char *)"WEATHER_WIND_SPEED", TFLOAT);
644	}
645
646	if (cData[WINDDIRE].colnum < 0) {
647	findData(WINDDIRE, (char *)"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((char )(char )"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((char *)"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
1030	int 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((char )(char )"OBSERVER", TSTRING, observer); // Shared.
1054	readParm((char )(char )"PROJID", TSTRING, project); // Shared.
1055	readParm((char )(char )"TELESCOP", TSTRING, telescope); // Core.
1056
1057	antPos[0] = 0.0;
1058	antPos[1] = 0.0;
1059	antPos[2] = 0.0;
1060	if (readParm((char *)"ANTENNA_POSITION", TDOUBLE, antPos)) {
1061	readParm((char *)"OBSGEO-X", TDOUBLE, antPos); // Additional.
1062	readParm((char *)"OBSGEO-Y", TDOUBLE, antPos + 1); // Additional.
1063	readParm((char *)"OBSGEO-Z", TDOUBLE, antPos + 2); // Additional.
1064	}
1065
1066	if (antPos[0] == 0.0) {
1067	if (cGBT) {
1068	readParm( (char )(char )"SITELONG", TDOUBLE, antPos ) ;
1069	readParm( (char )(char )"SITELAT", TDOUBLE, antPos+1 ) ;
1070	readParm( (char )(char )"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((char )(char )"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((char )(char )"EQUINOX", TSTRING, strtmp);
1124	sscanf(strtmp, "%f", &equinox);
1125	}
1126
1127	if (readParm((char )(char )"RADESYS", TSTRING, radecsys)) { // Additional.
1128	if (readParm((char )(char )"RADECSYS", TSTRING, radecsys)) { // Additional.
1129	strcpy(radecsys, "");
1130	}
1131	}
1132
1133	if (readParm((char )(char )"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((char )(char )"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((char )(char )"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
1227	int 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((char *)"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
1304	int 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.0PI) spin -= 2.0PI;
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
1534	cleanup:
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
1553	int 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 && numHDU > 1 ) {
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((char *)"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] + iPolnChan;
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((char *)"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] + iPolnChan;
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] + iPolnChan;
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] + iPolnChan;
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] + iPolnChan;
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
2182	void 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
2205	void 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
2228	void 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
2282	void 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
2310	int 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
2369	int 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
2381	int 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
2448	int 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
2512	int 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
2565	int 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 = (hh60 + 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
2623	int 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((char *)"OBS_NAME", TSTRING, cRow, chars);
2634	} else {
2635	readData((char *)"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((char *)"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
2749	float 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
2786	long 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
2840	Vector<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
2895	Double 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( (char )(char )"SITELONG", TDOUBLE, antPos ) ;
2919	readParm( (char )(char )"SITELAT", TDOUBLE, antPos+1 ) ;
2920	readParm( (char )(char )"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
2962	Double 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
2988	Int 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.

Context Navigation

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

Download in other formats: