Context Navigation

source: trunk/external/atnf/PKSIO/SDFITSreader.cc @ 1509

Visit:

Last change on this file since 1509 was 1509, checked in by Malte Marquarding, 15 years ago
make gcc-4.3 compliant; Mark C. still needs to fix char* cast deprecation warnings
File size: 54.0 KB

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

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

Download in other formats: