Changeset 1635 for trunk/external


Ignore:
Timestamp:
09/25/09 11:47:11 (15 years ago)
Author:
Malte Marquarding
Message:

Update from livedata CVS

Location:
trunk/external/atnf/PKSIO
Files:
14 edited

Legend:

Unmodified
Added
Removed
  • trunk/external/atnf/PKSIO/FITSreader.h

    r1452 r1635  
    2727//#                        AUSTRALIA
    2828//#
    29 //# $Id: FITSreader.h,v 19.9 2008-11-17 06:28:04 cal103 Exp $
     29//# $Id: FITSreader.h,v 19.10 2008-11-27 04:28:24 cal103 Exp $
    3030//#---------------------------------------------------------------------------
    3131//# The FITSreader class is an abstract base class for the Parkes Multibeam
     
    9696    // Coordinate systems are
    9797    //   0: equatorial (RA,Dec),
    98     //   1: vertical (Az,El),
    99     //   2: feed-plane.
     98    //   1: horizontal (Az,El),
     99    //   2: feed-plane,
     100    //   3: zenithal position angle of feed and elevation, (ZPA,El).
    100101    int select(
    101102        const int startChan[],
  • trunk/external/atnf/PKSIO/MBFITSreader.cc

    r1452 r1635  
    2727//#                        AUSTRALIA
    2828//#
    29 //# $Id: MBFITSreader.cc,v 19.54 2008-11-17 06:51:55 cal103 Exp $
     29//# $Id: MBFITSreader.cc,v 19.55 2009-01-20 06:45:33 cal103 Exp $
    3030//#---------------------------------------------------------------------------
    3131//# The MBFITSreader class reads single dish RPFITS files (such as Parkes
     
    141141  int jstat = -3;
    142142  if (rpfitsin(jstat)) {
    143     sprintf(cMsg, "ERROR: failed to open MBFITS file\n       %s", rpname);
     143    sprintf(cMsg, "ERROR: Failed to open MBFITS file\n       %s", rpname);
    144144    logMsg(cMsg);
    145145    return 1;
     
    159159  jstat = -1;
    160160  if (rpfitsin(jstat)) {
    161     sprintf(cMsg, "ERROR: failed to read MBFITS header in file\n"
     161    sprintf(cMsg, "ERROR: Failed to read MBFITS header in file\n"
    162162                  "       %s", rpname);
    163163    logMsg(cMsg);
     
    205205
    206206  if (cNBeam <= 0) {
    207     logMsg("ERROR, couldn't determine number of beams.");
     207    logMsg("ERROR: Couldn't determine number of beams.");
    208208    close();
    209209    return 1;
     
    337337  // Read the first syscal record.
    338338  if (rpget(1, cEOS)) {
    339     logMsg("ERROR, failed to read first syscal record.");
     339    logMsg("ERROR: Failed to read first syscal record.");
    340340    close();
    341341    return 1;
     
    373373{
    374374  if (!cMBopen) {
    375     logMsg("ERROR, an MBFITS file has not been opened.");
     375    logMsg("ERROR: An MBFITS file has not been opened.");
    376376    return 1;
    377377  }
     
    515515
    516516  if (!cMBopen) {
    517     logMsg("ERROR, an MBFITS file has not been opened.");
     517    logMsg("ERROR: An MBFITS file has not been opened.");
    518518    return 1;
    519519  }
     
    588588          cXpolOff = new int[cNIF];
    589589
    590           int simulIF = 0;
    591590          int maxChan = 0;
    592591          int maxXpol = 0;
    593592
     593          cSimulIF = 0;
    594594          for (int iIF = 0; iIF < cNIF; iIF++) {
    595595            if (cIFs[iIF]) {
     
    617617
    618618              // Maximum number of selected IFs in any simultaneous set.
    619               simulIF = max(simulIF, cIFSel[iIF]+1);
     619              cSimulIF = max(cSimulIF, cIFSel[iIF]+1);
    620620
    621621              // Maximum memory required for any simultaneous set.
     
    644644
    645645          if (cNBin > 1 && cNBeamSel > 1) {
    646             logMsg("ERROR, cannot handle binning mode for multiple beams.");
     646            logMsg("ERROR: Cannot handle binning mode for multiple beams.\n"
     647                   "       Select a single beam for input.");
    647648            close();
    648649            return 1;
     
    657658          // Allocate memory for spectral arrays.
    658659          for (int ibuff = 0; ibuff < nBuff; ibuff++) {
    659             cBuffer[ibuff].setNIFs(simulIF);
     660            cBuffer[ibuff].setNIFs(cSimulIF);
    660661            cBuffer[ibuff].allocate(0, maxChan, maxXpol);
    661662
    662663            // Signal that this IF in this buffer has been flushed.
    663             for (int iIF = 0; iIF < simulIF; iIF++) {
     664            for (int iIF = 0; iIF < cSimulIF; iIF++) {
    664665              cBuffer[ibuff].IFno[iIF] = 0;
    665666            }
     
    679680          cCycleNo = 0;
    680681          cPrevUTC = -1.0;
     682        }
     683
     684        // Apply beam and IF selection before the change-of-day test to allow
     685        // a single selected beam and IF to be handled in binning-mode.
     686        beamNo = int(cBaseline / 256.0);
     687        if (beamNo == 1) {
     688          // Store the position of beam 1 for grid convergence corrections.
     689          cRA0  = cU;
     690          cDec0 = cV;
     691        }
     692        iBeamSel = cBeamSel[beamNo-1];
     693        if (iBeamSel < 0) continue;
     694
     695        // Sanity check (mainly for MOPS).
     696        if (cIFno > cNIF) continue;
     697
     698        // Apply IF selection.
     699        iIFSel = cIFSel[cIFno - 1];
     700        if (iIFSel < 0) continue;
     701
     702
     703        if (cNBin > 1) {
     704          // Binning mode: correct the time.
     705          cUTC += param_.intbase * (cBin - (cNBin + 1)/2.0);
    681706        }
    682707
     
    689714          // the start of the next.
    690715#ifdef PKSIO_DEBUG
    691           fprintf(stderr, "Change-of-day on cUTC: %.1f -> %.1f",
     716          fprintf(stderr, "Change-of-day on cUTC: %.1f -> %.1f\n",
    692717            cPrevUTC, cUTC);
    693718#endif
     
    707732        }
    708733
    709         if (cNBin > 1) {
    710           // Binning mode: correct the time.
    711           cUTC += param_.intbase * (cBin - (cNBin + 1)/2.0);
    712         }
    713 
    714734        // New integration cycle?
    715735        if ((cUTC+cod) > cPrevUTC) {
     
    717737          cPrevUTC = cUTC + 0.0001;
    718738        }
    719 
    720         // Apply beam selection.
    721         beamNo = int(cBaseline / 256.0);
    722         if (beamNo == 1) {
    723           // Store the position of beam 1 for grid convergence corrections.
    724           cRA0  = cU;
    725           cDec0 = cV;
    726         }
    727         iBeamSel = cBeamSel[beamNo-1];
    728         if (iBeamSel < 0) continue;
    729 
    730         // Sanity check (mainly for MOPS).
    731         if (cIFno > cNIF) continue;
    732 
    733         // Apply IF selection.
    734         iIFSel = cIFSel[cIFno - 1];
    735         if (iIFSel < 0) continue;
    736739
    737740        sprintf(cDateObs, "%-10.10s", names_.datobs);
     
    784787          iMBuff = cBuffer + iBeamSel + cNBeamSel*cFlushBin;
    785788
    786           // iMBuff->nIF is set to zero (below) to signal that all IFs in
    787           // an integration have been flushed.
     789          // iMBuff->nIF is decremented (below) and if zero signals that all
     790          // IFs in an integration have been flushed.
    788791          if (iMBuff->nIF) {
    789792            if (cycleNo == 0 || iMBuff->cycleNo < cycleNo) {
     
    805808
    806809        // Find the IF to flush.
    807         for (; cFlushIF < iMBuff->nIF; cFlushIF++) {
     810        for (; cFlushIF < cSimulIF; cFlushIF++) {
    808811          if (iMBuff->IFno[cFlushIF]) break;
    809812        }
     
    827830
    828831
    829     if (cFlushing && cFlushBin == 0 && cFlushIF == 0 && cInterp) {
     832    if (cInterp && cFlushing == 1) {
    830833      // Start of flush cycle, interpolate the beam position.
    831834      //
     
    10981101#endif
    10991102      }
     1103
     1104      cFlushing = 2;
    11001105    }
    11011106
     
    11141119      iMBuff->IFno[cFlushIF] = 0;
    11151120
    1116       if (cFlushIF == iMBuff->nIF - 1) {
    1117         // Signal that all IFs in this buffer location have been flushed.
    1118         iMBuff->nIF = 0;
    1119 
    1120         // Stop cEOS being set when the next integration is read.
     1121      iMBuff->nIF--;
     1122      if (iMBuff->nIF == 0) {
     1123        // All IFs in this buffer location have been flushed.  Stop cEOS
     1124        // being set when the next integration is read.
    11211125        iMBuff->cycleNo = 0;
    11221126
  • trunk/external/atnf/PKSIO/MBFITSreader.h

    r1452 r1635  
    2727//#                        AUSTRALIA
    2828//#
    29 //# $Id: MBFITSreader.h,v 19.21 2008-11-17 06:33:10 cal103 Exp $
     29//# $Id: MBFITSreader.h,v 19.22 2009-01-20 06:36:03 cal103 Exp $
    3030//#---------------------------------------------------------------------------
    3131//# The MBFITSreader class reads single dish RPFITS files (such as Parkes
     
    116116    char   cDateObs[12];
    117117    int    *cBeamSel, *cChanOff, cFirst, *cIFSel, cInterp, cIntTime, cMBopen,
    118            cMopra, cNBeamSel, cNBin, cRetry, cSUpos, *cXpolOff;
     118           cMopra, cNBeamSel, cNBin, cRetry, cSimulIF, cSUpos, *cXpolOff;
    119119
    120120    // The data has to be bufferred to allow positions to be interpolated.
  • trunk/external/atnf/PKSIO/MBrecord.cc

    r1452 r1635  
    22//# MBrecord.cc: Class to store an MBFITS single-dish data record.
    33//#---------------------------------------------------------------------------
    4 //# Copyright (C) 2000-2008
     4//# Copyright (C) 2000-2009
    55//# Mark Calabretta, ATNF
    66//#
     
    2727//#                        AUSTRALIA
    2828//#
    29 //# $Id: MBrecord.cc,v 19.12 2008-11-17 06:52:35 cal103 Exp $
     29//# $Id: MBrecord.cc,v 19.13 2009-03-24 06:15:33 cal103 Exp $
    3030//#---------------------------------------------------------------------------
    3131//# The MBrecord class stores an MBFITS single-dish data record.
     
    9393    xcalfctr = new float[nif][2];
    9494    baseLin  = new float[nif][2][2];
    95     baseSub  = new float[nif][2][9];
     95    baseSub  = new float[nif][2][24];
    9696    spectra  = new float*[nif];
    9797    flagged  = new unsigned char*[nif];
     
    260260      baseLin[iIF][ipol][1] = other.baseLin[iIF][ipol][1];
    261261
    262       for (int j = 0; j < 9; j++) {
     262      for (int j = 0; j < 24; j++) {
    263263        baseSub[iIF][ipol][j] = other.baseSub[iIF][ipol][j];
    264264      }
     
    380380    baseLin[0][ipol][1] = other.baseLin[iIF][ipol][1];
    381381
    382     for (int j = 0; j < 9; j++) {
     382    for (int j = 0; j < 24; j++) {
    383383      baseSub[0][ipol][j] = other.baseSub[iIF][ipol][j];
    384384    }
  • trunk/external/atnf/PKSIO/MBrecord.h

    r1452 r1635  
    22//# MBrecord.h: Class to store an MBFITS single-dish data record.
    33//#---------------------------------------------------------------------------
    4 //# Copyright (C) 2000-2008
     4//# Copyright (C) 2000-2009
    55//# Mark Calabretta, ATNF
    66//#
     
    2727//#                        AUSTRALIA
    2828//#
    29 //# $Id: MBrecord.h,v 19.14 2008-11-17 06:36:12 cal103 Exp $
     29//# $Id: MBrecord.h,v 19.15 2009-03-24 06:15:33 cal103 Exp $
    3030//#---------------------------------------------------------------------------
    3131//# The MBrecord class stores an MBFITS single-dish data record.
     
    134134    int    haveBase;            // Are baseline parameters present?
    135135    float  (*baseLin)[2][2];    // Linear baseline fit for each polarization.
    136     float  (*baseSub)[2][9];    // Polynomial baseline subtracted.
     136    float  (*baseSub)[2][24];   // Polynomial baseline subtracted.
    137137    int    haveSpectra;         // Is spectral data present?
    138138    float* *spectra;            // Spectra for each polarization, Jy.
  • trunk/external/atnf/PKSIO/PKSFITSreader.cc

    r1452 r1635  
    11//# PKSFITSreader.cc: Class to read Parkes multibeam data from a FITS file.
    22//#---------------------------------------------------------------------------
    3 //# Copyright (C) 2000-2008
     3//# Copyright (C) 2000-2009
    44//# Associated Universities, Inc. Washington DC, USA.
    55//#
     
    2525//#                        Charlottesville, VA 22903-2475 USA
    2626//#
    27 //# $Id: PKSFITSreader.cc,v 19.21 2008-11-17 06:54:25 cal103 Exp $
     27//# $Id: PKSFITSreader.cc,v 19.22 2009-03-24 06:15:33 cal103 Exp $
    2828//#---------------------------------------------------------------------------
    2929//# Original: 2000/08/02, Mark Calabretta, ATNF
     
    446446  if (cMBrec.haveBase) {
    447447    pksrec.baseLin.resize(2,nPol);
    448     pksrec.baseSub.resize(9,nPol);
     448    pksrec.baseSub.resize(24,nPol);
    449449
    450450    for (uInt ipol = 0; ipol < nPol; ipol++) {
     
    452452      pksrec.baseLin(1,ipol) = cMBrec.baseLin[0][ipol][1];
    453453
    454       for (uInt j = 0; j < 9; j++) {
     454      for (uInt j = 0; j < 24; j++) {
    455455        pksrec.baseSub(j,ipol) = cMBrec.baseSub[0][ipol][j];
    456456      }
  • trunk/external/atnf/PKSIO/PKSMS2reader.cc

    r1452 r1635  
    22//# PKSMS2reader.cc: Class to read Parkes Multibeam data from a v2 MS.
    33//#---------------------------------------------------------------------------
    4 //# Copyright (C) 2000-2008
     4//# Copyright (C) 2000-2009
    55//# Associated Universities, Inc. Washington DC, USA.
    66//#
     
    2626//#                        Charlottesville, VA 22903-2475 USA
    2727//#
    28 //# $Id: PKSMS2reader.cc,v 19.21 2008-11-17 06:55:18 cal103 Exp $
     28//# $Id: PKSMS2reader.cc,v 19.22 2009-03-24 06:15:33 cal103 Exp $
    2929//#---------------------------------------------------------------------------
    3030//# Original: 2000/08/03, Mark Calabretta, ATNF
     
    656656    cBaseLinCol.get(cIdx, pksrec.baseLin);
    657657
    658     pksrec.baseSub.resize(9,cNPol(iIF));
     658    pksrec.baseSub.resize(24,cNPol(iIF));
    659659    cBaseSubCol.get(cIdx, pksrec.baseSub);
    660660
  • trunk/external/atnf/PKSIO/PKSMS2writer.cc

    r1452 r1635  
    22//# PKSMS2writer.cc: Class to write Parkes multibeam data to a measurementset.
    33//#---------------------------------------------------------------------------
    4 //# Copyright (C) 2000-2008
     4//# Copyright (C) 2000-2009
    55//# Associated Universities, Inc. Washington DC, USA.
    66//#
     
    2626//#                        Charlottesville, VA 22903-2475 USA
    2727//#
    28 //# $Id: PKSMS2writer.cc,v 19.14 2008-11-17 06:56:13 cal103 Exp $
     28//# $Id: PKSMS2writer.cc,v 19.15 2009-03-24 06:15:33 cal103 Exp $
    2929//#---------------------------------------------------------------------------
    3030
     
    123123                IPosition(2,2,maxNPol), ColumnDesc::Direct));
    124124    pksDesc.addColumn(ArrayColumnDesc<Float>("BASESUB", "Baseline subtracted",
    125                 IPosition(2,9,maxNPol), ColumnDesc::Direct));
     125                IPosition(2,24,maxNPol), ColumnDesc::Direct));
    126126  }
    127127
     
    469469    cBaseLinCol->put(irow, pksrec.baseLin);
    470470
    471     if (pksrec.baseSub.nrow() == 9) {
     471    if (pksrec.baseSub.nrow() == 24) {
    472472      cBaseSubCol->put(irow, pksrec.baseSub);
    473473
    474474    } else {
    475       Matrix<Float> tmp(9, 2, 0.0f);
     475      Matrix<Float> tmp(24, 2, 0.0f);
    476476      for (Int ipol = 0; ipol < nPol; ipol++) {
    477477        for (uInt j = 0; j < pksrec.baseSub.nrow(); j++) {
  • trunk/external/atnf/PKSIO/PKSSDwriter.cc

    r1452 r1635  
    22//# PKSSDwriter.cc: Class to write Parkes multibeam data to an SDFITS file.
    33//#---------------------------------------------------------------------------
    4 //# Copyright (C) 2000-2008
     4//# Copyright (C) 2000-2009
    55//# Associated Universities, Inc. Washington DC, USA.
    66//#
     
    2626//#                        Charlottesville, VA 22903-2475 USA
    2727//#
    28 //# $Id: PKSSDwriter.cc,v 19.15 2008-11-17 06:56:50 cal103 Exp $
     28//# $Id: PKSSDwriter.cc,v 19.16 2009-03-24 06:15:33 cal103 Exp $
    2929//#---------------------------------------------------------------------------
    3030
     
    239239        mbrec.baseSub[0][ipol][j] = pksrec.baseSub(j,ipol);
    240240      }
    241       for (uInt j = pksrec.baseSub.nrow(); j < 9; j++) {
     241      for (uInt j = pksrec.baseSub.nrow(); j < 24; j++) {
    242242        mbrec.baseSub[0][ipol][j] = 0.0f;
    243243      }
  • trunk/external/atnf/PKSIO/PKSmsg.cc

    r1466 r1635  
    3636
    3737#include <string>
    38 #include <cstring>
    3938
    4039//------------------------------------------------------------- PKSmsg::PKSmsg
  • trunk/external/atnf/PKSIO/PKSreader.h

    r1452 r1635  
    2626//#                        Charlottesville, VA 22903-2475 USA
    2727//#
    28 //# $Id: PKSreader.h,v 19.22 2008-11-17 06:44:34 cal103 Exp $
     28//# $Id: PKSreader.h,v 19.23 2008-11-27 04:28:24 cal103 Exp $
    2929//#---------------------------------------------------------------------------
    3030//# Original: 2000/08/02, Mark Calabretta, ATNF
     
    134134    // Coordinate system selection (only supported for SDFITS input):
    135135    //   0: equatorial (RA,Dec),
    136     //   1: vertical (Az,El),
    137     //   2: feed-plane.
     136    //   1: horizontal (Az,El),
     137    //   2: feed-plane,
     138    //   3: zenithal position angle of feed and elevation, (ZPA,El).
    138139    virtual uInt select(
    139140        const Vector<Bool> beamSel,
  • trunk/external/atnf/PKSIO/SDFITSreader.cc

    r1509 r1635  
    11//#---------------------------------------------------------------------------
    2 //# SDFITSreader.cc: ATNF CFITSIO interface class for SDFITS input.
     2//# SDFITSreader.cc: ATNF interface class for SDFITS input using CFITSIO.
    33//#---------------------------------------------------------------------------
    4 //# Copyright (C) 2000-2008
     4//# Copyright (C) 2000-2009
    55//# Associated Universities, Inc. Washington DC, USA.
    66//#
     
    2626//#                        Charlottesville, VA 22903-2475 USA
    2727//#
    28 //# $Id: SDFITSreader.cc,v 19.33 2008-11-17 06:58:34 cal103 Exp $
     28//# $Id: SDFITSreader.cc,v 19.43 2009-05-06 03:28:17 cal103 Exp $
    2929//#---------------------------------------------------------------------------
    3030//# The SDFITSreader class reads single dish FITS files such as those written
     
    4141#include <casa/math.h>
    4242#include <casa/stdio.h>
    43 #include <cstring>
    4443
    4544#include <algorithm>
     
    6463const double D2R = PI / 180.0;
    6564
     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
    6675//------------------------------------------------- SDFITSreader::SDFITSreader
    6776
     
    6978{
    7079  // Default constructor.
    71   cSDptr = 0;
     80  cSDptr = 0x0;
    7281
    7382  // Allocate space for data descriptors.
     
    157166    // Arecibo ALFA data of some kind.
    158167    cALFA = 1;
    159     for (int iBeam = 0; iBeam < 8; iBeam++) {
    160       for (int iPol = 0; iPol < 2; iPol++) {
    161         cALFAcalOn[iBeam][iPol]  = 0.0f;
    162         cALFAcalOff[iBeam][iPol] = 0.0f;
    163 
    164         // Nominal factor to calibrate spectra in Jy.
    165         cALFAcal[iBeam][iPol] = 3.0f;
    166       }
     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;
    167180    }
    168181  }
     
    174187         strncmp(telescope, "NRAO_GBT", 8) == 0;
    175188
    176   cRow = 0;
    177 
    178189
    179190  // Check that the DATA array column is present.
     
    181192  haveSpectra = cHaveSpectra = cData[DATA].colnum > 0;
    182193
     194  cNAxisTime = 0;
    183195  if (cHaveSpectra) {
    184196    // Find the number of data axes (must be the same for each IF).
    185     cNAxis = 5;
    186     if (readDim(DATA, 1, &cNAxis, cNAxes)) {
     197    cNAxes = 5;
     198    if (readDim(DATA, 1, &cNAxes, cNAxis)) {
    187199      logMsg();
    188200      close();
     
    193205      // ALFA BDFITS: variable length arrays don't actually vary and there is
    194206      // no TDIM (or MAXISn) card; use the LAGS_IN value.
    195       cNAxis = 5;
    196       readParm("LAGS_IN", TLONG, cNAxes);
    197       cNAxes[1] = 1;
    198       cNAxes[2] = 1;
    199       cNAxes[3] = 1;
    200       cNAxes[4] = 1;
    201       cData[DATA].nelem = cNAxes[0];
    202     }
    203 
    204     if (cNAxis < 4) {
     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) {
    205217      // Need at least four axes (for now).
    206218      logMsg("ERROR: DATA array contains fewer than four axes.");
    207219      close();
    208220      return 1;
    209     } else if (cNAxis > 5) {
     221    } else if (cNAxes > 5) {
    210222      // We support up to five axes.
    211223      logMsg("ERROR: DATA array contains more than five axes.");
     
    229241    readParm("DATAXED", TSTRING, dataxed);
    230242
    231     for (int iaxis = 0; iaxis < 5; iaxis++) cNAxes[iaxis] = 0;
    232     sscanf(dataxed, "(%ld,%ld,%ld,%ld,%ld)", cNAxes, cNAxes+1, cNAxes+2,
    233       cNAxes+3, cNAxes+4);
     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);
    234246    for (int iaxis = 4; iaxis > -1; iaxis--) {
    235       if (cNAxes[iaxis] == 0) cNAxis = iaxis;
     247      if (cNAxis[iaxis] == 0) cNAxes = iaxis;
    236248    }
    237249  }
     
    244256  // Find required DATA array axes.
    245257  char ctype[5][72];
    246   for (int iaxis = 0; iaxis < cNAxis; iaxis++) {
     258  for (int iaxis = 0; iaxis < cNAxes; iaxis++) {
    247259    strcpy(ctype[iaxis], "");
    248260    readParm(CTYPE[iaxis], TSTRING, ctype[iaxis]);      // Core.
     
    255267  }
    256268
    257   char *fqCRPIX  = 0;
    258269  char *fqCRVAL  = 0;
    259270  char *fqCDELT  = 0;
     271  char *fqCRPIX  = 0;
    260272  char *raCRVAL  = 0;
    261273  char *decCRVAL = 0;
    262274  char *timeCRVAL = 0;
     275  char *timeCDELT = 0;
     276  char *timeCRPIX = 0;
    263277  char *beamCRVAL = 0;
    264278
    265   for (int iaxis = 0; iaxis < cNAxis; iaxis++) {
     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++) {
    266287    if (strncmp(ctype[iaxis], "FREQ", 4) == 0) {
    267       cReqax[0] = iaxis;
    268       fqCRPIX  = CRPIX[iaxis];
    269       fqCRVAL  = CRVAL[iaxis];
    270       fqCDELT  = CDELT[iaxis];
     288      cFreqAxis = iaxis;
     289      fqCRVAL   = CRVAL[iaxis];
     290      fqCDELT   = CDELT[iaxis];
     291      fqCRPIX   = CRPIX[iaxis];
    271292
    272293    } else if (strncmp(ctype[iaxis], "STOKES", 6) == 0) {
    273       cReqax[1] = iaxis;
     294      cStokesAxis = iaxis;
    274295
    275296    } else if (strncmp(ctype[iaxis], "RA", 2) == 0) {
    276       cReqax[2] = iaxis;
    277       raCRVAL  = CRVAL[iaxis];
     297      cRaAxis  = iaxis;
     298      raCRVAL   = CRVAL[iaxis];
    278299
    279300    } else if (strncmp(ctype[iaxis], "DEC", 3) == 0) {
    280       cReqax[3] = iaxis;
    281       decCRVAL = CRVAL[iaxis];
     301      cDecAxis = iaxis;
     302      decCRVAL  = CRVAL[iaxis];
    282303
    283304    } else if (strcmp(ctype[iaxis], "TIME") == 0) {
    284       // TIME (UTC seconds since midnight) can be a keyword or axis type.
     305      // TIME (UTC seconds since midnight); axis type, if present, takes
     306      // precedence over keyword.
     307      cTimeAxis = iaxis;
    285308      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      }
    286318
    287319    } else if (strcmp(ctype[iaxis], "BEAM") == 0) {
    288320      // BEAM can be a keyword or axis type.
     321      cBeamAxis = iaxis;
    289322      beamCRVAL = CRVAL[iaxis];
    290323    }
     
    293326  if (cALFA_BD) {
    294327    // Fixed in ALFA CIMAFITS.
    295     cReqax[2] = 2;
     328    cRaAxis = 2;
    296329    raCRVAL = "CRVAL2A";
    297330
    298     cReqax[3] = 3;
     331    cDecAxis = 3;
    299332    decCRVAL = "CRVAL3A";
    300333  }
    301334
    302   // Check that all are present.
    303   for (int iaxis = 0; iaxis < 4; iaxis++) {
    304     if (cReqax[iaxis] < 0) {
    305       logMsg("ERROR: Could not find required DATA array axes.");
    306       close();
    307       return 1;
    308     }
     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;
    309340  }
    310341
     
    313344  findData(CYCLE,    "CYCLE",    TINT);         // Additional.
    314345  findData(DATE_OBS, "DATE-OBS", TSTRING);      // Core.
    315   findData(TIME,     "TIME",     TDOUBLE);      // 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
    316363  findData(EXPOSURE, "EXPOSURE", TFLOAT);       // Core.
    317364  findData(OBJECT,   "OBJECT",   TSTRING);      // Core.
     
    323370  findData(BEAM,     "BEAM",     TSHORT);       // Additional.
    324371  findData(IF,       "IF",       TSHORT);       // Additional.
    325   findData(FqRefPix,  fqCRPIX,   TFLOAT);       // Frequency reference pixel.
    326372  findData(FqRefVal,  fqCRVAL,   TDOUBLE);      // Frequency reference value.
    327373  findData(FqDelt,    fqCDELT,   TDOUBLE);      // Frequency increment.
     374  findData(FqRefPix,  fqCRPIX,   TFLOAT);       // Frequency reference pixel.
    328375  findData(RA,        raCRVAL,   TDOUBLE);      // Right ascension.
    329376  findData(DEC,      decCRVAL,   TDOUBLE);      // Declination.
     
    367414      findData(SCAN,  "SCAN_ID", TINT);
    368415      if (cALFA_CIMA > 1) {
     416        // Note that RECNUM increases by cNAxisTime per row.
    369417        findData(CYCLE, "RECNUM", TINT);
    370418      } else {
     
    410458    }
    411459    cIF_1rel = 0;
    412   }
    413 
    414   if (cData[TIME].colnum < 0) {
    415     if (timeCRVAL) {
    416       // There is a TIME axis.
    417       findData(TIME, timeCRVAL, TDOUBLE);
    418     }
    419460  }
    420461
     
    607648          if (cData[DATA].nelem < 0) {
    608649            // Variable dimension array.
    609             if (readDim(DATA, irow+1, &cNAxis, cNAxes)) {
     650            if (readDim(DATA, irow+1, &cNAxes, cNAxis)) {
    610651              logMsg();
    611652              close();
     
    619660            readParm("DATAXED", TSTRING, dataxed);
    620661
    621             sscanf(dataxed, "(%ld,%ld,%ld,%ld,%ld)", cNAxes, cNAxes+1,
    622               cNAxes+2, cNAxes+3, cNAxes+4);
     662            sscanf(dataxed, "(%ld,%ld,%ld,%ld,%ld)", cNAxis, cNAxis+1,
     663              cNAxis+2, cNAxis+3, cNAxis+4);
    623664          }
    624665        }
    625666
    626667        // Number of channels and polarizations.
    627         cNChan[iIF]    = cNAxes[cReqax[0]];
    628         cNPol[iIF]     = cNAxes[cReqax[1]];
     668        cNChan[iIF]    = cNAxis[cFreqAxis];
     669        cNPol[iIF]     = cNAxis[cStokesAxis];
    629670        cHaveXPol[iIF] = 0;
    630671
     
    665706
    666707    // Number of channels and polarizations.
    667     cNChan[0] = cNAxes[cReqax[0]];
    668     cNPol[0]  = cNAxes[cReqax[1]];
     708    cNChan[0] = cNAxis[cFreqAxis];
     709    cNPol[0]  = cNAxis[cStokesAxis];
    669710    cHaveXPol[0] = 0;
    670711  }
     
    726767
    727768  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;
    728781
    729782  return 0;
     
    825878      // No, try digging it out of the CTYPE card (AIPS convention).
    826879      char keyw[9], ctype[9];
    827       sprintf(keyw, "CTYPE%ld", cReqax[0]+1);
     880      sprintf(keyw, "CTYPE%ld", cFreqAxis+1);
    828881      readParm(keyw, TSTRING, ctype);
    829882
     
    860913
    861914  // Get parameters from first row of table.
    862   readData(DATE_OBS, 1, datobs);
    863   readData(TIME,     1, &utc);
     915  readTime(1, 1, datobs, utc);
    864916  readData(FqRefVal, 1, &refFreq);
    865917  readParm("BANDWID", TDOUBLE, &bandwidth);             // Core.
    866 
    867   if (cALFA_BD) utc *= 3600.0;
    868918
    869919  if (cStatus) {
    870920    logMsg();
    871921    return 1;
    872   }
    873 
    874   // Check DATE-OBS format.
    875   if (datobs[2] == '/') {
    876     // Translate an old-format DATE-OBS.
    877     datobs[9] = datobs[1];
    878     datobs[8] = datobs[0];
    879     datobs[2] = datobs[6];
    880     datobs[5] = datobs[3];
    881     datobs[3] = datobs[7];
    882     datobs[6] = datobs[4];
    883     datobs[7] = '-';
    884     datobs[4] = '-';
    885     datobs[1] = '9';
    886     datobs[0] = '1';
    887     datobs[10] = '\0';
    888 
    889   } else if (datobs[10] == 'T' && cData[TIME].colnum < 0) {
    890     // Dig UTC out of a new-format DATE-OBS.
    891     int   hh, mm;
    892     float ss;
    893     sscanf(datobs+11, "%d:%d:%f", &hh, &mm, &ss);
    894     utc = (hh*60 + mm)*60 + ss;
    895     datobs[10] = '\0';
    896922  }
    897923
     
    9931019
    9941020  // Find the number of rows selected.
    995   short *sel = new short[nRow];
    996   for (int irow = 0; irow < nRow; irow++) {
     1021  short *sel = new short[cNRow];
     1022  for (int irow = 0; irow < cNRow; irow++) {
    9971023    sel[irow] = 1;
    9981024  }
     
    10101036    }
    10111037
    1012     for (int irow = 0; irow < nRow; irow++) {
     1038    for (int irow = 0; irow < cNRow; irow++) {
    10131039      if (!cBeams[beamCol[irow]-cBeam_1rel]) {
    10141040        sel[irow] = 0;
     
    10301056    }
    10311057
    1032     for (int irow = 0; irow < nRow; irow++) {
     1058    for (int irow = 0; irow < cNRow; irow++) {
    10331059      if (!cIFs[IFCol[irow]-cIF_1rel]) {
    10341060        sel[irow] = 0;
     
    10401066
    10411067  nSel = 0;
    1042   for (int irow = 0; irow < nRow; irow++) {
     1068  for (int irow = 0; irow < cNRow; irow++) {
    10431069    nSel += sel[irow];
    10441070  }
     
    10461072
    10471073  // Find the time range assuming the data is in chronological order.
    1048   readData(DATE_OBS, 1,    dateSpan[0]);
    1049   readData(DATE_OBS, nRow, dateSpan[1]);
    1050   readData(TIME, 1,    utcSpan);
    1051   readData(TIME, nRow, utcSpan+1);
    1052 
    1053   if (cALFA_BD) {
    1054     utcSpan[0] *= 3600.0;
    1055     utcSpan[1] *= 3600.0;
    1056   }
    1057 
    1058   // Check DATE-OBS format.
    1059   for (int i = 0; i < 2; i++) {
    1060     if (dateSpan[0][2] == '/') {
    1061       // Translate an old-format DATE-OBS.
    1062       dateSpan[i][9] = dateSpan[i][1];
    1063       dateSpan[i][8] = dateSpan[i][0];
    1064       dateSpan[i][2] = dateSpan[i][6];
    1065       dateSpan[i][5] = dateSpan[i][3];
    1066       dateSpan[i][3] = dateSpan[i][7];
    1067       dateSpan[i][6] = dateSpan[i][4];
    1068       dateSpan[i][7] = '-';
    1069       dateSpan[i][4] = '-';
    1070       dateSpan[i][1] = '9';
    1071       dateSpan[i][0] = '1';
    1072       dateSpan[i][10] = '\0';
    1073     }
    1074 
    1075     if (dateSpan[i][10] == 'T' && cData[TIME].colnum < 0) {
    1076       // Dig UTC out of a new-format DATE-OBS.
    1077       int   hh, mm;
    1078       float ss;
    1079       sscanf(dateSpan[i]+11, "%d:%d:%f", &hh, &mm, &ss);
    1080       utcSpan[i] = (hh*60 + mm)*60 + ss;
    1081     }
    1082   }
     1074  readTime(1, 1, dateSpan[0], utcSpan[0]);
     1075  readTime(cNRow, cNAxisTime, dateSpan[1], utcSpan[1]);
    10831076
    10841077
     
    10881081
    10891082  if (cCoordSys == 1) {
    1090     // Vertical (Az,El).
    1091     if (cData[AZIMUTH].colnum  < 1 ||
    1092         cData[ELEVATIO].colnum < 1) {
     1083    // Horizontal (Az,El).
     1084    if (cData[AZIMUTH].colnum  < 0 ||
     1085        cData[ELEVATIO].colnum < 0) {
    10931086      logMsg("WARNING: Azimuth/elevation information absent.");
    10941087      cStatus = -1;
     
    10971090      float *az = new float[cNRow];
    10981091      float *el = new float[cNRow];
    1099       fits_read_col(cSDptr, TFLOAT, cData[AZIMUTH].colnum,  1, 1, nRow, 0, az,
    1100                     &anynul, &cStatus);
    1101       fits_read_col(cSDptr, TFLOAT, cData[ELEVATIO].colnum, 1, 1, nRow, 0, el,
    1102                     &anynul, &cStatus);
     1092      readCol(AZIMUTH,  az);
     1093      readCol(ELEVATIO, el);
    11031094
    11041095      if (!cStatus) {
    1105         for (int irow = 0; irow < nRow; irow++) {
     1096        for (int irow = 0; irow < cNRow; irow++) {
    11061097          if (sel[irow]) {
    11071098            positions[isel++] = az[irow] * D2R;
     
    11151106    }
    11161107
     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
    11171155  } else {
    11181156    double *ra  = new double[cNRow];
    11191157    double *dec = new double[cNRow];
    1120     fits_read_col(cSDptr, TDOUBLE, cData[RA].colnum,  1, 1, nRow, 0, ra,
    1121                   &anynul, &cStatus);
    1122     fits_read_col(cSDptr, TDOUBLE, cData[DEC].colnum, 1, 1, nRow, 0, dec,
    1123                   &anynul, &cStatus);
     1158    readCol(RA,  ra);
     1159    readCol(DEC, dec);
     1160
    11241161    if (cStatus) {
    11251162      delete [] ra;
     
    11291166
    11301167    if (cALFA_BD) {
    1131       for (int irow = 0; irow < nRow; irow++) {
     1168      for (int irow = 0; irow < cNRow; irow++) {
    11321169        // Convert hours to degrees.
    11331170        ra[irow] *= 15.0;
     
    11371174    if (cCoordSys == 0) {
    11381175      // Equatorial (RA,Dec).
    1139       for (int irow = 0; irow < nRow; irow++) {
     1176      for (int irow = 0; irow < cNRow; irow++) {
    11401177        if (sel[irow]) {
    11411178          positions[isel++] =  ra[irow] * D2R;
     
    11481185      if (cData[OBJ_RA].colnum   < 0 ||
    11491186          cData[OBJ_DEC].colnum  < 0 ||
    1150           cData[PARANGLE].colnum < 1 ||
    1151           cData[FOCUSROT].colnum < 1) {
     1187          cData[PARANGLE].colnum < 0 ||
     1188          cData[FOCUSROT].colnum < 0) {
    11521189        logMsg("WARNING: Insufficient information to compute feed-plane\n"
    11531190               "         coordinates.");
     
    11601197        float  *rot = new float[cNRow];
    11611198
    1162         if (cData[OBJ_RA].colnum == 0) {
    1163           // Header keyword.
    1164           readData(OBJ_RA, 0, srcRA);
    1165           for (int irow = 1; irow < nRow; irow++) {
    1166             srcRA[irow] = *srcRA;
    1167           }
    1168         } else {
    1169           // Table column.
    1170           fits_read_col(cSDptr, TDOUBLE, cData[OBJ_RA].colnum,   1, 1, nRow,
    1171                         0, srcRA,  &anynul, &cStatus);
    1172         }
    1173 
    1174         if (cData[OBJ_DEC].colnum == 0) {
    1175           // Header keyword.
    1176           readData(OBJ_DEC, 0, srcDec);
    1177           for (int irow = 1; irow < nRow; irow++) {
    1178             srcDec[irow] = *srcDec;
    1179           }
    1180         } else {
    1181           // Table column.
    1182           fits_read_col(cSDptr, TDOUBLE, cData[OBJ_DEC].colnum,  1, 1, nRow,
    1183                         0, srcDec, &anynul, &cStatus);
    1184         }
    1185 
    1186         fits_read_col(cSDptr, TFLOAT,  cData[PARANGLE].colnum, 1, 1, nRow, 0,
    1187                       par,    &anynul, &cStatus);
    1188         fits_read_col(cSDptr, TFLOAT,  cData[FOCUSROT].colnum, 1, 1, nRow, 0,
    1189                       rot,    &anynul, &cStatus);
     1199        readCol(OBJ_RA,   srcRA);
     1200        readCol(OBJ_DEC,  srcDec);
     1201        readCol(PARANGLE, par);
     1202        readCol(FOCUSROT, rot);
    11901203
    11911204        if (!cStatus) {
    1192           for (int irow = 0; irow < nRow; irow++) {
     1205          for (int irow = 0; irow < cNRow; irow++) {
    11931206            if (sel[irow]) {
    11941207              // Convert to feed-plane coordinates.
     
    12471260  // Find the next selected beam and IF.
    12481261  short iBeam = 0, iIF = 0;
    1249   while (++cRow <= cNRow) {
     1262  while (1) {
     1263    if (++cTimeIdx > cNAxisTime) {
     1264      if (++cRow > cNRow) break;
     1265      cTimeIdx = 1;
     1266    }
     1267
    12501268    if (cData[BEAM].colnum > 0) {
    12511269      readData(BEAM, cRow, &iBeam);
     
    12691287          char chars[32];
    12701288          readData(OBSMODE, cRow, chars);
    1271           if (strcmp(chars, "CAL") == 0) {
     1289          if (strcmp(chars, "DROP") == 0) {
     1290            // Completely flagged integration.
     1291            continue;
     1292
     1293          } else if (strcmp(chars, "CAL") == 0) {
     1294            sReset = 1;
    12721295            if (cALFA_CIMA > 1) {
    12731296              for (short iPol = 0; iPol < cNPol[iIF]; iPol++) {
     
    12801303              continue;
    12811304            }
     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            }
    12821321          }
    12831322        }
     
    13121351  // Times.
    13131352  char datobs[32];
    1314   readData(DATE_OBS, cRow,  datobs);
    1315   readData(TIME,     cRow, &mbrec.utc);
    1316   if (cALFA_BD) mbrec.utc *= 3600.0;
    1317 
    1318   if (datobs[2] == '/') {
    1319     // Translate an old-format DATE-OBS.
    1320     datobs[9] = datobs[1];
    1321     datobs[8] = datobs[0];
    1322     datobs[2] = datobs[6];
    1323     datobs[5] = datobs[3];
    1324     datobs[3] = datobs[7];
    1325     datobs[6] = datobs[4];
    1326     datobs[7] = '-';
    1327     datobs[4] = '-';
    1328     datobs[1] = '9';
    1329     datobs[0] = '1';
    1330 
    1331   } else if (datobs[10] == 'T' && cData[TIME].colnum < 0) {
    1332     // Dig UTC out of a new-format DATE-OBS.
    1333     int   hh, mm;
    1334     float ss;
    1335     sscanf(datobs+11, "%d:%d:%f", &hh, &mm, &ss);
    1336     mbrec.utc = (hh*60 + mm)*60 + ss;
    1337   }
    1338 
    1339   datobs[10] = '\0';
     1353  readTime(cRow, cTimeIdx, datobs, mbrec.utc);
    13401354  strcpy(mbrec.datobs, datobs);
    13411355
    13421356  if (cData[CYCLE].colnum > 0) {
    13431357    readData(CYCLE, cRow, &mbrec.cycleNo);
     1358    mbrec.cycleNo += cTimeIdx - 1;
    13441359    if (cALFA_BD) mbrec.cycleNo++;
    13451360  } else {
     
    14711486
    14721487  // Read data, sectioning and transposing it in the process.
    1473   long *blc = new long[cNAxis+1];
    1474   long *trc = new long[cNAxis+1];
    1475   long *inc = new long[cNAxis+1];
    1476   for (int iaxis = 0; iaxis <= cNAxis; iaxis++) {
     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++) {
    14771492    blc[iaxis] = 1;
    14781493    trc[iaxis] = 1;
     
    14801495  }
    14811496
    1482   blc[cReqax[0]] = std::min(startChan, endChan);
    1483   trc[cReqax[0]] = std::max(startChan, endChan);
    1484   blc[cNAxis] = cRow;
    1485   trc[cNAxis] = cRow;
     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;
    14861505
    14871506  mbrec.haveSpectra = cGetSpectra;
     
    14891508    int  anynul;
    14901509
    1491     for (int ipol = 0; ipol < nPol; ipol++) {
    1492       blc[cReqax[1]] = ipol+1;
    1493       trc[cReqax[1]] = ipol+1;
     1510    for (int iPol = 0; iPol < nPol; iPol++) {
     1511      blc[cStokesAxis] = iPol+1;
     1512      trc[cStokesAxis] = iPol+1;
    14941513
    14951514      if (cALFA && cALFA_CIMA < 2) {
    14961515        // ALFA data: polarizations are stored in successive rows.
    1497         blc[cReqax[1]] = 1;
    1498         trc[cReqax[1]] = 1;
    1499 
    1500         if (ipol) {
     1516        blc[cStokesAxis] = 1;
     1517        trc[cStokesAxis] = 1;
     1518
     1519        if (iPol) {
    15011520          if (++cRow > cNRow) {
    15021521            return -1;
    15031522          }
    15041523
    1505           blc[cNAxis] = cRow;
    1506           trc[cNAxis] = cRow;
     1524          blc[cNAxes] = cRow;
     1525          trc[cNAxes] = cRow;
    15071526        }
    15081527
    15091528      } else if (cData[DATA].nelem < 0) {
    15101529        // Variable dimension array; get axis lengths.
    1511         int  naxis = 5, status;
    1512 
    1513         if ((status = readDim(DATA, cRow, &naxis, cNAxes))) {
     1530        int naxes = 5, status;
     1531
     1532        if ((status = readDim(DATA, cRow, &naxes, cNAxis))) {
    15141533          logMsg();
    15151534
    1516         } else if ((status = (naxis != cNAxis))) {
     1535        } else if ((status = (naxes != cNAxes))) {
    15171536          logMsg("ERROR: DATA array dimensions changed.");
    15181537        }
     
    15261545      }
    15271546
    1528       if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAxis, cNAxes,
    1529           blc, trc, inc, 0, mbrec.spectra[0] + ipol*nChan, &anynul,
     1547      if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAxes, cNAxis,
     1548          blc, trc, inc, 0, mbrec.spectra[0] + iPol*nChan, &anynul,
    15301549          &cStatus)) {
    15311550        logMsg();
     
    15381557      if (endChan < startChan) {
    15391558        // Reverse the spectrum.
    1540         float *iptr = mbrec.spectra[0] + ipol*nChan;
     1559        float *iptr = mbrec.spectra[0] + iPol*nChan;
    15411560        float *jptr = iptr + nChan - 1;
    15421561        float *mid  = iptr + nChan/2;
     
    15501569      if (cALFA) {
    15511570        // ALFA data, rescale the spectrum.
    1552         float *chan  = mbrec.spectra[0] + ipol*nChan;
     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;
    15531587        float *chanN = chan + nChan;
    15541588        while (chan < chanN) {
    15551589          // Approximate conversion to Jy.
    1556           *(chan++) *= cALFAcal[iBeam][iIF];
     1590          *(chan++) *= factor;
    15571591        }
    15581592      }
    15591593
    1560       if (mbrec.tsys[0][ipol] == 0.0) {
     1594      if (mbrec.tsys[0][iPol] == 0.0) {
    15611595        // Compute Tsys as the average across the spectrum.
    1562         float *chan  = mbrec.spectra[0] + ipol*nChan;
     1596        float *chan  = mbrec.spectra[0] + iPol*nChan;
    15631597        float *chanN = chan + nChan;
    1564         float *tsys = mbrec.tsys[0] + ipol;
     1598        float *tsys = mbrec.tsys[0] + iPol;
    15651599        while (chan < chanN) {
    15661600          *tsys += *(chan++);
     
    15721606      // Read data flags.
    15731607      if (cData[FLAGGED].colnum > 0) {
    1574         if (fits_read_subset_byt(cSDptr, cData[FLAGGED].colnum, cNAxis,
    1575             cNAxes, blc, trc, inc, 0, mbrec.flagged[0] + ipol*nChan, &anynul,
     1608        if (fits_read_subset_byt(cSDptr, cData[FLAGGED].colnum, cNAxes,
     1609            cNAxis, blc, trc, inc, 0, mbrec.flagged[0] + iPol*nChan, &anynul,
    15761610            &cStatus)) {
    15771611          logMsg();
     
    15841618        if (endChan < startChan) {
    15851619          // Reverse the flag vector.
    1586           unsigned char *iptr = mbrec.flagged[0] + ipol*nChan;
     1620          unsigned char *iptr = mbrec.flagged[0] + iPol*nChan;
    15871621          unsigned char *jptr = iptr + nChan - 1;
    15881622          for (int ichan = 0; ichan < nChan/2; ichan++) {
     
    15951629      } else {
    15961630        // All channels are unflagged by default.
    1597         unsigned char *iptr = mbrec.flagged[0] + ipol*nChan;
     1631        unsigned char *iptr = mbrec.flagged[0] + iPol*nChan;
    15981632        for (int ichan = 0; ichan < nChan; ichan++) {
    15991633          *(iptr++) = 0;
     
    17081742    int status = 0;
    17091743    fits_close_file(cSDptr, &status);
    1710     cSDptr = 0;
     1744    cSDptr = 0x0;
    17111745
    17121746    if (cBeams)     delete [] cBeams;
     
    17911825}
    17921826
     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
    17931855//------------------------------------------------------ SDFITSreader::readDim
    17941856
     
    17981860        int  iData,
    17991861        long iRow,
    1800         int *naxis,
    1801         long naxes[])
     1862        int *naxes,
     1863        long naxis[])
    18021864{
    18031865  int colnum = cData[iData].colnum;
     
    18061868  }
    18071869
    1808   int maxdim = *naxis;
     1870  int maxdim = *naxes;
    18091871  if (cData[iData].tdimcol < 0) {
    18101872    // No TDIMnnn column for this array.
    18111873    if (cData[iData].nelem < 0) {
    18121874      // Variable length array; read the array descriptor.
    1813       *naxis = 1;
     1875      *naxes = 1;
    18141876      long dummy;
    1815       if (fits_read_descript(cSDptr, colnum, iRow, naxes, &dummy, &cStatus)) {
     1877      if (fits_read_descript(cSDptr, colnum, iRow, naxis, &dummy, &cStatus)) {
    18161878        return 1;
    18171879      }
     
    18191881    } else {
    18201882      // Read the repeat count from TFORMnnn.
    1821       if (fits_read_tdim(cSDptr, colnum, maxdim, naxis, naxes, &cStatus)) {
     1883      if (fits_read_tdim(cSDptr, colnum, maxdim, naxes, naxis, &cStatus)) {
    18221884        return 1;
    18231885      }
     
    18371899
    18381900    tp++;
    1839     *naxis = 0;
     1901    *naxes = 0;
    18401902    for (size_t j = 1; j < strlen(tdimval); j++) {
    18411903      if (tdimval[j] == ',' || tdimval[j] == ')') {
    1842         sscanf(tp, "%ld", naxes + (*naxis)++);
     1904        sscanf(tp, "%ld", naxis + (*naxes)++);
    18431905        if (tdimval[j] == ')') break;
    18441906        tp = tdimval + j + 1;
     
    18751937  findCol(name, &colnum);
    18761938
    1877   if (colnum > 0) {
     1939  if (colnum > 0 && iRow > 0) {
    18781940    // Read the first value from the specified row of the table.
    18791941    int  coltype;
     
    19382000        void *value)
    19392001{
    1940   char *name  = cData[iData].name;
    19412002  int  type   = cData[iData].type;
    19422003  int  colnum = cData[iData].colnum;
    1943   long nelem  = cData[iData].nelem;
    1944 
    1945   if (colnum > 0) {
     2004
     2005  if (colnum > 0 && iRow > 0) {
    19462006    // Read the required number of values from the specified row of the table.
     2007    long nelem = cData[iData].nelem;
    19472008    int anynul;
    19482009    if (type == TSTRING) {
     
    19732034  } else if (colnum == 0) {
    19742035    // Read keyword value.
     2036    char *name  = cData[iData].name;
    19752037    fits_read_key(cSDptr, type, name, value, 0, &cStatus);
    19762038
     
    19932055}
    19942056
    1995 //------------------------------------------------------ SDFITSreader::findCol
    1996 
    1997 // Locate a parameter in the SDFITS file.
    1998 
    1999 void SDFITSreader::findCol(
    2000         char *name,
    2001         int *colnum)
     2057//------------------------------------------------------ SDFITSreader::readCol
     2058
     2059// Read a scalar column from the SDFITS file.
     2060
     2061int SDFITSreader::readCol(
     2062        int  iData,
     2063        void *value)
    20022064{
    2003   *colnum = 0;
    2004   int status = 0;
    2005   fits_get_colnum(cSDptr, CASESEN, name, colnum, &status);
    2006 
    2007   if (status) {
    2008     // Not a real column - maybe it's virtual.
    2009     char card[81];
    2010 
    2011     status = 0;
    2012     fits_read_card(cSDptr, name, card, &status);
    2013     if (status) {
    2014       // Not virtual either.
    2015       *colnum = -1;
    2016     }
    2017 
    2018     // Clear error messages.
    2019     fits_clear_errmsg();
    2020   }
     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;
    20212140}
    20222141
     
    20472166
    20482167  // Read cal data.
    2049   long *blc = new long[cNAxis+1];
    2050   long *trc = new long[cNAxis+1];
    2051   long *inc = new long[cNAxis+1];
    2052   for (int iaxis = 0; iaxis <= cNAxis; iaxis++) {
     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++) {
    20532172    blc[iaxis] = 1;
    20542173    trc[iaxis] = 1;
     
    20602179  int endChan   = cEndChan[iIF];
    20612180
    2062   blc[cNAxis] = cRow;
    2063   trc[cNAxis] = cRow;
    2064   blc[cReqax[0]] = std::min(startChan, endChan);
    2065   trc[cReqax[0]] = std::max(startChan, endChan);
     2181  blc[cFreqAxis] = std::min(startChan, endChan);
     2182  trc[cFreqAxis] = std::max(startChan, endChan);
    20662183  if (cALFA_CIMA > 1) {
    20672184    // CIMAFITS 2.x has a legitimate STOKES axis...
    2068     blc[cReqax[1]] = iPol+1;
    2069     trc[cReqax[1]] = iPol+1;
     2185    blc[cStokesAxis] = iPol+1;
     2186    trc[cStokesAxis] = iPol+1;
    20702187  } else {
    20712188    // ...older ALFA data does not.
    2072     blc[cReqax[1]] = 1;
    2073     trc[cReqax[1]] = 1;
    2074   }
     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;
    20752198
    20762199  float spectrum[endChan];
    20772200  int anynul;
    2078   if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAxis, cNAxes,
     2201  if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAxes, cNAxis,
    20792202      blc, trc, inc, 0, spectrum, &anynul, &cStatus)) {
    20802203    logMsg();
     
    20852208  }
    20862209
     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
    20872220  // Average the spectrum.
    20882221  float mean = 1e9f;
     
    20982231      if (*chan < discrim) {
    20992232        nChan++;
    2100         sum += *chan;
     2233        sum += *chan * factor;
    21012234      }
    21022235    }
     
    21062239
    21072240  if (calOn) {
    2108     cALFAcalOn[iBeam][iPol]  += mean;
     2241    sALFAcalOn[iBeam][iPol]  *= sALFAcalNon[iBeam][iPol];
     2242    sALFAcalOn[iBeam][iPol]  += mean;
     2243    sALFAcalOn[iBeam][iPol]  /= ++sALFAcalNon[iBeam][iPol];
    21092244  } else {
    2110     cALFAcalOff[iBeam][iPol] += mean;
    2111   }
    2112 
    2113   if (cALFAcalOn[iBeam][iPol] != 0.0f &&
    2114       cALFAcalOff[iBeam][iPol] != 0.0f) {
     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]) {
    21152251    // Tcal should come from the TCAL table, it varies weakly with beam,
    21162252    // polarization, and frequency.  However, TCAL is not written properly.
    21172253    float Tcal = 12.0f;
    2118     cALFAcal[iBeam][iPol] = Tcal / (cALFAcalOn[iBeam][iPol] -
    2119                                     cALFAcalOff[iBeam][iPol]);
     2254    sALFAcal[iBeam][iPol] = Tcal / (sALFAcalOn[iBeam][iPol] -
     2255                                    sALFAcalOff[iBeam][iPol]);
    21202256
    21212257    // Scale from K to Jy; the gain also varies weakly with beam,
    21222258    // polarization, frequency, and zenith angle.
    21232259    float fluxCal = 10.0f;
    2124     cALFAcal[iBeam][iPol] /= fluxCal;
    2125 
    2126     cALFAcalOn[iBeam][iPol]  = 0.0f;
    2127     cALFAcalOff[iBeam][iPol] = 0.0f;
     2260    sALFAcal[iBeam][iPol] /= fluxCal;
    21282261  }
    21292262
    21302263  return 0;
    21312264}
     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
  • trunk/external/atnf/PKSIO/SDFITSreader.h

    r1452 r1635  
    22//# SDFITSreader.h: ATNF CFITSIO interface class for SDFITS input.
    33//#---------------------------------------------------------------------------
    4 //# Copyright (C) 2000-2008
     4//# Copyright (C) 2000-2009
    55//# Associated Universities, Inc. Washington DC, USA.
    66//#
     
    2626//#                        Charlottesville, VA 22903-2475 USA
    2727//#
    28 //# $Id: SDFITSreader.h,v 19.16 2008-11-17 06:47:05 cal103 Exp $
     28//# $Id: SDFITSreader.h,v 19.21 2009-03-18 07:11:51 cal103 Exp $
    2929//#---------------------------------------------------------------------------
    3030//# The SDFITSreader class reads single dish FITS files such as those written
     
    108108
    109109  private:
    110     int      cCycleNo, cExtraSysCal, cNAxis, cStatus;
    111     long     cNAxes[5], cNRow, cReqax[4], cRow;
     110    int      cCycleNo, cExtraSysCal, cNAxes, cStatus;
     111    long     cBeamAxis, cDecAxis, cFreqAxis, cNAxis[5], cNAxisTime, cNRow,
     112             cRaAxis, cRow, cStokesAxis, cTimeAxis, cTimeIdx;
    112113    double   cLastUTC;
    113114    fitsfile *cSDptr;
     
    118119
    119120    enum {SCAN, CYCLE, DATE_OBS, TIME, EXPOSURE, OBJECT, OBJ_RA, OBJ_DEC,
    120           RESTFRQ, OBSMODE, BEAM, IF, FqRefPix, FqRefVal, FqDelt, RA, DEC,
    121           SCANRATE, TSYS, CALFCTR, XCALFCTR, BASELIN, BASESUB, DATA, FLAGGED,
    122           DATAXED, XPOLDATA, REFBEAM, TCAL, TCALTIME, AZIMUTH, ELEVATIO,
    123           PARANGLE, FOCUSAXI, FOCUSTAN, FOCUSROT, TAMBIENT, PRESSURE,
    124           HUMIDITY, WINDSPEE, WINDDIRE, NDATA};
     121          RESTFRQ, OBSMODE, BEAM, IF, FqRefVal, FqDelt, FqRefPix, RA, DEC,
     122          TimeRefVal, TimeDelt, TimeRefPix, SCANRATE, TSYS, CALFCTR, XCALFCTR,
     123          BASELIN, BASESUB, DATA, FLAGGED, DATAXED, XPOLDATA, REFBEAM, TCAL,
     124          TCALTIME, AZIMUTH, ELEVATIO, PARANGLE, FOCUSAXI, FOCUSTAN, FOCUSROT,
     125          TAMBIENT, PRESSURE, HUMIDITY, WINDSPEE, WINDDIRE, NDATA};
    125126
    126127    // Message handling.
     
    128129
    129130    void findData(int iData, char *name, int type);
     131    void  findCol(char *name, int *colnum);
    130132    int   readDim(int iData, long iRow, int *naxis, long naxes[]);
    131133    int  readParm(char *name, int type, void *value);
    132134    int  readData(char *name, int type, long iRow, void *value);
    133135    int  readData(int iData, long iRow, void *value);
    134     void  findCol(char *name, int *colnum);
     136    int  readCol(int iData, void *value);
     137    int  readTime(long iRow, int iPix, char *datobs, double &utc);
    135138
    136     // These are for ALFA data: "BDFITS" or "CIMAFITS".
     139    // These are for ALFA data: "BDFITS" or "CIMAFITS".  Statics are required
     140    // for CIMAFITS v2.0 because CAL ON/OFF data is split into separate files.
     141    static int  sInit, sReset;
     142    static int  (*sALFAcalNon)[2], (*sALFAcalNoff)[2];
     143    static float (*sALFAcal)[2], (*sALFAcalOn)[2], (*sALFAcalOff)[2];
     144
    137145    int   cALFA, cALFA_BD, cALFA_CIMA, cALFAscan, cScanNo;
    138     float cALFAcal[8][2], cALFAcalOn[8][2], cALFAcalOff[8][2];
     146    float cALFAacc;
    139147    int   alfaCal(short iBeam, short iIF, short iPol);
     148    float alfaGain(float zd);
    140149
    141150    // These are for GBT data.
  • trunk/external/atnf/PKSIO/SDFITSwriter.cc

    r1466 r1635  
    22//# SDFITSwriter.cc: ATNF CFITSIO interface class for SDFITS output.
    33//#---------------------------------------------------------------------------
    4 //# Copyright (C) 2000-2008
     4//# Copyright (C) 2000-2009
    55//# Mark Calabretta, ATNF
    66//#
     
    2727//#                        AUSTRALIA
    2828//#
    29 //# $Id: SDFITSwriter.cc,v 19.15 2008-11-17 06:58:58 cal103 Exp $
     29//# $Id: SDFITSwriter.cc,v 19.17 2009-09-18 07:29:26 cal103 Exp $
    3030//#---------------------------------------------------------------------------
    3131//# Original: 2000/07/24, Mark Calabretta, ATNF
     
    3737
    3838#include <casa/iostream.h>
    39 #include <cstring>
    4039
    4140#include <algorithm>
     
    174173  char version[7];
    175174  char date[11];
    176   sscanf("$Revision: 19.15 $", "%*s%s", version);
    177   sscanf("$Date: 2008-11-17 06:58:58 $", "%*s%s", date);
     175  sscanf("$Revision: 19.17 $", "%*s%s", version);
     176  sscanf("$Date: 2009-09-18 07:29:26 $", "%*s%s", date);
    178177  sprintf(text, "SDFITSwriter (v%s, %s)", version, date);
    179178  fits_write_key_str(cSDptr, "ORIGIN", text, "output class", &cStatus);
     
    231230
    232231  // CYCLE (additional, real).
    233   fits_insert_col(cSDptr, ++ncol, "CYCLE", "1I", &cStatus);
     232  fits_insert_col(cSDptr, ++ncol, "CYCLE", "1J", &cStatus);
    234233
    235234  // DATE-OBS (core, real).
     
    391390
    392391    // BASESUB (additional, real).
    393     sprintf(tform, "%dE", 9*maxNPol);
     392    sprintf(tform, "%dE", 24*maxNPol);
    394393    fits_insert_col(cSDptr, ++ncol, "BASESUB", tform, &cStatus);
    395     tdim[0] = 9;
     394    tdim[0] = 24;
    396395    fits_write_tdim(cSDptr, ncol, 2, tdim, &cStatus);
    397396  }
     
    681680
    682681    // BASESUB.
    683     fits_write_col_flt(cSDptr, ++icol, cRow, 1, 9*nPol, mbrec.baseSub[0][0],
     682    fits_write_col_flt(cSDptr, ++icol, cRow, 1, 24*nPol, mbrec.baseSub[0][0],
    684683                       &cStatus);
    685684  }
Note: See TracChangeset for help on using the changeset viewer.