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

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

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