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

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

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