Context Navigation

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

Visit:

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

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

Download in other formats: