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

Last change on this file since 3116 was 1737, checked in by Malte Marquarding, 15 years ago

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