Context Navigation

source: trunk/external-alma/atnf/PKSIO/SDFITSreader.cc@ 2960

Visit:

Last change on this file since 2960 was 2940, checked in by Malte Marquarding, 11 years ago
Finally fix aall warning: deprecated conversion from string constant to 'char*'
File size: 67.2 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: