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

Last change on this file since 1627 was 1509, checked in by Malte Marquarding, 16 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
48class 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.
61const double PI = 3.141592653589793238462643;
62
63// Factor to convert radians to degrees.
64const double D2R = PI / 180.0;
65
66//------------------------------------------------- SDFITSreader::SDFITSreader
67
68SDFITSreader::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
93SDFITSreader::~SDFITSreader()
94{
95 close();
96
97 delete [] cData;
98}
99
100//--------------------------------------------------------- SDFITSreader::open
101
102// Open an SDFITS file for reading.
103
104int 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
736int 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 = (hh*60 + 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
905int 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
980int 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] = (hh*60 + 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.0*PI) spin -= 2.0*PI;
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
1220cleanup:
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
1239int 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 = (hh*60 + 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] + ipol*nChan;
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] + ipol*nChan;
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] + ipol*nChan;
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] + ipol*nChan;
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] + ipol*nChan;
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
1705void 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
1725void 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
1743void 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
1797int 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
1856int 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
1868int 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
1935int 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
1999void 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
2027int 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.