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

Last change on this file since 1454 was 1452, checked in by Malte Marquarding, 16 years ago

update from livedata CVS
File size: 54.0 KB
RevLine 
[1325]1//#---------------------------------------------------------------------------
2//# SDFITSreader.cc: ATNF CFITSIO interface class for SDFITS input.
3//#---------------------------------------------------------------------------
[1427]4//# Copyright (C) 2000-2008
[1325]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//#
[1452]28//# $Id: SDFITSreader.cc,v 19.33 2008-11-17 06:58:34 cal103 Exp $
[1325]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
[1452]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>
[1325]40
41#include <casa/math.h>
42#include <casa/stdio.h>
43
[1452]44#include <algorithm>
45#include <strings.h>
[1325]46
47class 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.
[1399]56 char units[32]; // Units from TUNITn keyword.
[1325]57};
58
59// Numerical constants.
60const double PI = 3.141592653589793238462643;
61
62// Factor to convert radians to degrees.
63const double D2R = PI / 180.0;
64
65//------------------------------------------------- SDFITSreader::SDFITSreader
66
67SDFITSreader::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;
[1452]85
86 // By default, messages are written to stderr.
87 initMsg();
[1325]88}
89
90//------------------------------------------------ SDFITSreader::~SDFITSreader
91
92SDFITSreader::~SDFITSreader()
93{
94 close();
95
96 delete [] cData;
97}
98
99//--------------------------------------------------------- SDFITSreader::open
100
101// Open an SDFITS file for reading.
102
103int 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{
[1452]116 // Clear the message stack.
117 clearMsg();
118
[1325]119 if (cSDptr) {
120 close();
121 }
122
123 // Open the SDFITS file.
124 cStatus = 0;
125 if (fits_open_file(&cSDptr, sdName, READONLY, &cStatus)) {
[1452]126 sprintf(cMsg, "ERROR: Failed to open SDFITS file\n %s", sdName);
127 logMsg(cMsg);
[1325]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
[1452]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
[1325]149 } else {
[1452]150 logMsg("ERROR: Failed to locate SDFITS binary table.");
[1325]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)) {
[1452]186 logMsg();
[1325]187 close();
188 return 1;
189 }
190
191 if (cALFA_BD) {
[1399]192 // ALFA BDFITS: variable length arrays don't actually vary and there is
[1325]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).
[1452]205 logMsg("ERROR: DATA array contains fewer than four axes.");
[1325]206 close();
207 return 1;
208 } else if (cNAxis > 5) {
209 // We support up to five axes.
[1452]210 logMsg("ERROR: DATA array contains more than five axes.");
[1325]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) {
[1452]221 logMsg("ERROR: DATA array column absent from binary table.");
[1325]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) {
[1452]251 logMsg();
[1325]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) {
[1452]304 logMsg("ERROR: Could not find required DATA array axes.");
[1325]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) {
[1452]354 logMsg();
[1325]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;
[1452]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) {
[1325]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;
[1452]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) {
[1325]391 if (beamCRVAL) {
392 // There is a BEAM axis.
393 findData(BEAM, beamCRVAL, TDOUBLE);
394 } else {
[1452]395 // ms2sdfits output, 0-relative "feed" number.
396 findData(BEAM, "MAIN_FEED1", TSHORT);
[1325]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.
[1452]405 if (cALFA_CIMA > 1) {
406 findData(IF, "IFN", TSHORT);
407 } else {
408 findData(IF, "IFVAL", TSHORT);
409 }
[1325]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) {
[1452]492 logMsg("ERROR: Table contains no entries.");
[1325]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;
[1452]506 logMsg();
[1325]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;
[1452]521 logMsg("ERROR: SDFITS file contains invalid beam number.");
[1325]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;
[1452]562 logMsg();
[1325]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;
[1452]577 logMsg("ERROR: SDFITS file contains invalid IF number.");
[1325]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)) {
[1452]609 logMsg();
[1325]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)) {
[1452]637 logMsg();
[1325]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
[1452]671 if (cALFA && cALFA_CIMA < 2) {
672 // Older ALFA data labels each polarization as a separate IF.
[1325]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
735int SDFITSreader::getHeader(
736 char observer[32],
737 char project[32],
738 char telescope[32],
739 double antPos[3],
740 char obsMode[32],
[1399]741 char bunit[32],
[1325]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
[1399]790 // Brightness unit.
[1452]791 if (cData[DATAXED].colnum >= 0) {
792 strcpy(bunit, "Jy");
793 } else {
794 strcpy(bunit, cData[DATA].units);
795 }
796
[1399]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
[1325]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) {
[1452]856 logMsg();
[1325]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) {
[1452]869 logMsg();
[1325]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 = (hh*60 + 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
904int 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;
[1452]923 logMsg();
[1325]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
979int 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
[1452]999 int anynul;
[1325]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;
[1452]1007 logMsg();
[1325]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;
[1452]1027 logMsg();
[1325]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] = (hh*60 + mm)*60 + ss;
1080 }
1081 }
1082
1083
1084 // Retrieve positions for selected data.
1085 int isel = 0;
1086 positions = new double[2*nSel];
1087
[1452]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;
[1325]1094
[1452]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);
[1325]1102
[1452]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 }
[1325]1111
[1452]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;
[1325]1133 }
1134 }
1135
[1452]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 }
[1325]1144
[1452]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.0*PI) spin -= 2.0*PI;
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;
[1325]1212 }
1213 }
[1452]1214
1215 delete [] ra;
1216 delete [] dec;
[1325]1217 }
1218
[1452]1219cleanup:
[1325]1220 delete [] sel;
1221
[1452]1222 if (cStatus) {
1223 nSel = 0;
1224 delete [] positions;
1225 logMsg();
1226 cStatus = 0;
1227 return 1;
1228 }
1229
1230 return 0;
[1325]1231}
1232
1233
1234//--------------------------------------------------------- SDFITSreader::read
1235
1236// Read the next data record.
1237
1238int SDFITSreader::read(
[1452]1239 MBrecord &mbrec)
[1325]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) {
[1452]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 }
[1325]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 = (hh*60 + 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.
[1399]1355 readData(OBJECT, cRow, mbrec.srcName);
[1325]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 }
[1452]1392 mbrec.paRate = 0.0f;
[1325]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) {
[1452]1430 logMsg();
[1325]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) {
[1452]1467 logMsg();
[1325]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
[1452]1494 if (cALFA && cALFA_CIMA < 2) {
[1325]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))) {
[1452]1513 logMsg();
[1325]1514
1515 } else if ((status = (naxis != cNAxis))) {
[1452]1516 logMsg("ERROR: DATA array dimensions changed.");
[1325]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)) {
[1452]1530 logMsg();
[1325]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] + ipol*nChan;
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] + ipol*nChan;
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] + ipol*nChan;
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)) {
[1452]1576 logMsg();
[1325]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] + ipol*nChan;
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] + ipol*nChan;
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)) {
[1452]1627 logMsg();
[1325]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) {
[1452]1658 logMsg();
[1325]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);
[1452]1666
[1325]1667 readData(AZIMUTH, cRow, &mbrec.azimuth);
1668 readData(ELEVATIO, cRow, &mbrec.elevation);
1669 readData(PARANGLE, cRow, &mbrec.parAngle);
[1452]1670
[1325]1671 readData(FOCUSAXI, cRow, &mbrec.focusAxi);
1672 readData(FOCUSTAN, cRow, &mbrec.focusTan);
1673 readData(FOCUSROT, cRow, &mbrec.focusRot);
[1452]1674
[1325]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) {
[1452]1693 logMsg();
[1325]1694 return 1;
1695 }
1696
1697 return 0;
1698}
1699
1700//-------------------------------------------------------- SDFITSreader::close
1701
1702// Close the SDFITS file.
1703
1704void 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
[1452]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
1724void 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
[1325]1738//----------------------------------------------------- SDFITSreader::findData
1739
1740// Locate a data item in the SDFITS file.
1741
1742void 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);
[1399]1759 fits_get_bcolparms(cSDptr, colnum, 0x0, cData[iData].units, 0x0, 0x0, 0x0,
1760 0x0, 0x0, 0x0, &cStatus);
[1325]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
1796int 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
1855int 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
1867int 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
1934int 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
1998void 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}
[1452]2021
2022//------------------------------------------------------ SDFITSreader::alfaCal
2023
2024// Process ALFA calibration data.
2025
2026int 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.