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

Last change on this file since 1635 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.