source: trunk/external-alma/atnf/PKSIO/SDFITSreader.cc @ 3067

Last change on this file since 3067 was 3067, checked in by Takeshi Nakazato, 8 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: Yes/No?

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...


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