source: trunk/external/atnf/PKSIO/MBFITSreader.cc@ 1555

//#---------------------------------------------------------------------------
//# MBFITSreader.cc: ATNF single-dish RPFITS reader.
//#---------------------------------------------------------------------------
//# Copyright (C) 2000-2008
//# Mark Calabretta, ATNF
//#
//# This library is free software; you can redistribute it and/or modify it
//# under the terms of the GNU Library General Public License as published by
//# the Free Software Foundation; either version 2 of the License, or (at your
//# option) any later version.
//#
//# This library is distributed in the hope that it will be useful, but WITHOUT
//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
//# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
//# License for more details.
//#
//# You should have received a copy of the GNU Library General Public License
//# along with this library; if not, write to the Free Software Foundation,
//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
//#
//# Correspondence concerning this software should be addressed as follows:
//#        Internet email: mcalabre@atnf.csiro.au.
//#        Postal address: Dr. Mark Calabretta,
//#                        Australia Telescope National Facility,
//#                        P.O. Box 76,
//#                        Epping, NSW, 2121,
//#                        AUSTRALIA
//#
//# $Id: MBFITSreader.cc,v 19.54 2008-11-17 06:51:55 cal103 Exp $
//#---------------------------------------------------------------------------
//# The MBFITSreader class reads single dish RPFITS files (such as Parkes
//# Multibeam MBFITS files).
//#
//# Original: 2000/07/28 Mark Calabretta
//#---------------------------------------------------------------------------

#include <atnf/pks/pks_maths.h>
#include <atnf/PKSIO/MBFITSreader.h>
#include <atnf/PKSIO/MBrecord.h>

#include <casa/math.h>
#include <casa/iostream.h>
#include <casa/stdio.h>
#include <casa/stdlib.h>
#include <casa/string.h>
#include <unistd.h>

#include <RPFITS.h>

using namespace std;

// Numerical constants.
const double PI = 3.141592653589793238462643;
const double TWOPI = 2.0 * PI;
const double HALFPI = PI / 2.0;
const double R2D = 180.0 / PI;

//------------------------------------------------- MBFITSreader::MBFITSreader

// Default constructor.

MBFITSreader::MBFITSreader(
        const int retry,
        const int interpolate)
{
  cRetry = retry;
  if (cRetry > 10) {
    cRetry = 10;
  }

  cInterp = interpolate;
  if (cInterp < 0 || cInterp > 2) {
    cInterp = 1;
  }

  // Initialize pointers.
  cBeams     = 0x0;
  cIFs       = 0x0;
  cNChan     = 0x0;
  cNPol      = 0x0;
  cHaveXPol  = 0x0;
  cStartChan = 0x0;
  cEndChan   = 0x0;
  cRefChan   = 0x0;

  cVis = 0x0;
  cWgt = 0x0;

  cBeamSel   = 0x0;
  cIFSel     = 0x0;
  cChanOff   = 0x0;
  cXpolOff   = 0x0;
  cBuffer    = 0x0;
  cPosUTC    = 0x0;

  cMBopen = 0;

  // Tell RPFITSIN not to report errors directly.
  iostat_.errlun = -1;

  // By default, messages are written to stderr.
  initMsg();
}

//------------------------------------------------ MBFITSreader::~MBFITSreader

// Destructor.

MBFITSreader::~MBFITSreader()
{
  close();
}

//--------------------------------------------------------- MBFITSreader::open

// Open the RPFITS file for reading.

int MBFITSreader::open(
        char *rpname,
        int  &nBeam,
        int* &beams,
        int  &nIF,
        int* &IFs,
        int* &nChan,
        int* &nPol,
        int* &haveXPol,
        int  &haveBase,
        int  &haveSpectra,
        int  &extraSysCal)
{
  // Clear the message stack.
  clearMsg();

  if (cMBopen) {
    close();
  }

  strcpy(names_.file, rpname);

  // Open the RPFITS file.
  int jstat = -3;
  if (rpfitsin(jstat)) {
    sprintf(cMsg, "ERROR: failed to open MBFITS file\n       %s", rpname);
    logMsg(cMsg);
    return 1;
  }

  cMBopen = 1;

  // Tell RPFITSIN that we want the OBSTYPE card.
  int j;
  param_.ncard = 1;
  for (j = 0; j < 80; j++) {
    names_.card[j] = ' ';
  }
  strncpy(names_.card, "OBSTYPE", 7);

  // Read the first header.
  jstat = -1;
  if (rpfitsin(jstat)) {
    sprintf(cMsg, "ERROR: failed to read MBFITS header in file\n"
                  "       %s", rpname);
    logMsg(cMsg);
    close();
    return 1;
  }

  // Mopra data has some peculiarities.
  cMopra = strncmp(names_.instrument, "ATMOPRA", 7) == 0;

  // Non-ATNF data may not store the position in (u,v,w).
  if (strncmp(names_.sta, "tid", 3) == 0) {
    sprintf(cMsg, "WARNING: Found Tidbinbilla data");
    cSUpos = 1;
  } else if (strncmp(names_.sta, "HOB", 3) == 0) {
    sprintf(cMsg, "WARNING: Found Hobart data");
    cSUpos = 1;
  } else if (strncmp(names_.sta, "CED", 3) == 0) {
    sprintf(cMsg, "WARNING: Found Ceduna data");
    cSUpos = 1;
  } else {
    cSUpos = 0;
  }

  if (cSUpos) {
    strcat(cMsg, ", using telescope position\n         from SU table.");
    logMsg(cMsg);
    cInterp = 0;
  }

  // Mean scan rate (for timestamp repairs).
  cNRate = 0;
  cAvRate[0] = 0.0;
  cAvRate[1] = 0.0;
  cCode5 = 0;


  // Find the maximum beam number.
  cNBeam = 0;
  for (int iBeam = 0; iBeam < anten_.nant; iBeam++) {
    if (anten_.ant_num[iBeam] > cNBeam) {
      cNBeam = anten_.ant_num[iBeam];
    }
  }

  if (cNBeam <= 0) {
    logMsg("ERROR, couldn't determine number of beams.");
    close();
    return 1;
  }

  // Construct the beam mask.
  cBeams = new int[cNBeam];
  for (int iBeam = 0; iBeam < cNBeam; iBeam++) {
    cBeams[iBeam] = 0;
  }

  // ...beams present in the data.
  for (int iBeam = 0; iBeam < anten_.nant; iBeam++) {
    // Guard against dubious beam numbers, e.g. zeroes in
    // 1999-09-29_1632_024848p14_071b.hpf and the four scans following.
    // Note that the actual beam number is decoded from the 'baseline' random
    // parameter for each spectrum and is only used for beam selection.
    int beamNo = anten_.ant_num[iBeam];
    if (beamNo != iBeam+1) {
      char sta[8];
      strncpy(sta, names_.sta+(8*iBeam), 8);
      char *cp = sta + 7;
      while (*cp == ' ') *(cp--) = '\0';

      sprintf(cMsg,
        "WARNING: RPFITSIN returned beam number %2d for AN table\n"
        "         entry %2d with name '%.8s'", beamNo, iBeam+1, sta);

      char text[8];
      sprintf(text, "MB%2.2d", iBeam+1);
      cp = cMsg + strlen(cMsg);
      if (strncmp(sta, text, 8) == 0) {
        beamNo = iBeam + 1;
        sprintf(cp, "; using beam number %2d.", beamNo);
      } else {
        sprintf(cp, ".");
      }

      logMsg(cMsg);
    }

    if (0 < beamNo && beamNo <= cNBeam) {
      cBeams[beamNo-1] = 1;
    }
  }

  // Passing back the address of the array allows PKSFITSreader::select() to
  // modify its elements directly.
  nBeam = cNBeam;
  beams = cBeams;


  // Number of IFs.
  cNIF = if_.n_if;
  cIFs = new int[cNIF];
  for (int iIF = 0; iIF < cNIF; iIF++) {
    cIFs[iIF] = 1;
  }

  // Passing back the address of the array allows PKSFITSreader::select() to
  // modify its elements directly.
  nIF = cNIF;
  IFs = cIFs;


  // Number of channels and polarizations.
  cNChan    = new int[cNIF];
  cNPol     = new int[cNIF];
  cHaveXPol = new int[cNIF];
  cGetXPol  = 0;

  int maxProd = 0;
  for (int iIF = 0; iIF < cNIF; iIF++) {
    cNChan[iIF] = if_.if_nfreq[iIF];
    cNPol[iIF]  = if_.if_nstok[iIF];
    cNChan[iIF] -= cNChan[iIF]%2;

    // Do we have cross-polarization data?
    if ((cHaveXPol[iIF] = cNPol[iIF] > 2)) {
      // Cross-polarization data is handled separately.
      cNPol[iIF] = 2;

      // Default is to get it if we have it.
      cGetXPol = 1;
    }

    // Maximum number of spectral products in any IF.
    int nProd = if_.if_nfreq[iIF] * if_.if_nstok[iIF];
    if (maxProd < nProd) maxProd = nProd;
  }

  // Allocate memory for RPFITSIN subroutine arguments.
  if (cVis) delete [] cVis;
  if (cWgt) delete [] cWgt;
  cVis = new float[2*maxProd];
  cWgt = new float[maxProd];

  nChan    = cNChan;
  nPol     = cNPol;
  haveXPol = cHaveXPol;


  // Default channel range selection.
  cStartChan = new int[cNIF];
  cEndChan   = new int[cNIF];
  cRefChan   = new int[cNIF];

  for (int iIF = 0; iIF < cNIF; iIF++) {
    cStartChan[iIF] = 1;
    cEndChan[iIF] = cNChan[iIF];
    cRefChan[iIF] = cNChan[iIF]/2 + 1;
  }

  cGetSpectra = 1;


  // No baseline parameters in MBFITS.
  haveBase = 0;

  // Always have spectra in MBFITS.
  haveSpectra = cHaveSpectra = 1;


  // Integration cycle time (s).
  cIntTime = param_.intime;

  // Can't deduce binning mode till later.
  cNBin = 0;


  // Read the first syscal record.
  if (rpget(1, cEOS)) {
    logMsg("ERROR, failed to read first syscal record.");
    close();
    return 1;
  }

  // Additional information for Parkes Multibeam data?
  extraSysCal = (sc_.sc_ant > anten_.nant);


  cFirst = 1;
  cEOF = 0;
  cFlushing = 0;

  return 0;
}

//---------------------------------------------------- MBFITSreader::getHeader

// Get parameters describing the data.

int MBFITSreader::getHeader(
        char   observer[32],
        char   project[32],
        char   telescope[32],
        double antPos[3],
        char   obsType[32],
        char   bunit[32],
        float  &equinox,
        char   radecsys[32],
        char   dopplerFrame[32],
        char   datobs[32],
        double &utc,
        double &refFreq,
        double &bandwidth)
{
  if (!cMBopen) {
    logMsg("ERROR, an MBFITS file has not been opened.");
    return 1;
  }

  sprintf(observer,  "%-16.16s", names_.rp_observer);
  sprintf(project,   "%-16.16s", names_.object);
  sprintf(telescope, "%-16.16s", names_.instrument);

  // Observatory coordinates (ITRF), in m.
  antPos[0] = doubles_.x[0];
  antPos[1] = doubles_.y[0];
  antPos[2] = doubles_.z[0];

  // This is the only sure way to identify the telescope, maybe.
  if (strncmp(names_.sta, "MB0", 3) == 0) {
    // Parkes Multibeam.
    sprintf(telescope, "%-16.16s", "ATPKSMB");
    antPos[0] = -4554232.087;
    antPos[1] =  2816759.046;
    antPos[2] = -3454035.950;

  } else if (strncmp(names_.sta, "HOH", 3) == 0) {
    // Parkes HOH receiver.
    sprintf(telescope, "%-16.16s", "ATPKSHOH");
    antPos[0] = -4554232.087;
    antPos[1] =  2816759.046;
    antPos[2] = -3454035.950;

  } else if (strncmp(names_.sta, "CA0", 3) == 0) {
    // An ATCA antenna, use the array centre position.
    sprintf(telescope, "%-16.16s", "ATCA");
    antPos[0] = -4750915.837;
    antPos[1] =  2792906.182;
    antPos[2] = -3200483.747;

    // ATCA-104.  Updated position at epoch 2007/06/24 from Chris Phillips.
    // antPos[0] = -4751640.182; // ± 0.008
    // antPos[1] =  2791700.322; // ± 0.006
    // antPos[2] = -3200490.668; // ± 0.007
    //
  } else if (strncmp(names_.sta, "MOP", 3) == 0) {
    // Mopra.  Updated position at epoch 2007/06/24 from Chris Phillips.
    sprintf(telescope, "%-16.16s", "ATMOPRA");
    antPos[0] = -4682769.444; // ± 0.009
    antPos[1] =  2802618.963; // ± 0.006
    antPos[2] = -3291758.864; // ± 0.008

  } else if (strncmp(names_.sta, "HOB", 3) == 0) {
    // Hobart.
    sprintf(telescope, "%-16.16s", "HOBART");
    antPos[0] = -3950236.735;
    antPos[1] =  2522347.567;
    antPos[2] = -4311562.569;

  } else if (strncmp(names_.sta, "CED", 3) == 0) {
    // Ceduna.  Updated position at epoch 2007/06/24 from Chris Phillips.
    sprintf(telescope, "%-16.16s", "CEDUNA");
    antPos[0] = -3753443.168; // ± 0.017
    antPos[1] =  3912709.794; // ± 0.017
    antPos[2] = -3348067.060; // ± 0.016

  } else if (strncmp(names_.sta, "tid", 3) == 0) {
    // DSS.
    sprintf(telescope, "%-16.16s", "DSS-43");
    antPos[0] = -4460894.727;
    antPos[1] =  2682361.530;
    antPos[2] = -3674748.424;
  }

  // Observation type.
  int j;
  for (j = 0; j < 31; j++) {
    obsType[j] = names_.card[11+j];
    if (obsType[j] == '\'') break;
  }
  obsType[j] = '\0';

  // Brightness unit.
  sprintf(bunit, "%-16.16s", names_.bunit);
  if (strcmp(bunit, "JY") == 0) {
    bunit[1] = 'y';
  } else if (strcmp(bunit, "JY/BEAM") == 0) {
    strcpy(bunit, "Jy/beam");
  }

  // Coordinate frames.
  equinox = 2000.0f;
  strcpy(radecsys, "FK5");
  strcpy(dopplerFrame, "TOPOCENT");

  // Time at start of observation.
  sprintf(datobs, "%-10.10s", names_.datobs);
  utc = cUTC;

  // Spectral parameters.
  refFreq   = doubles_.if_freq[0];
  bandwidth = doubles_.if_bw[0];

  return 0;
}

//-------------------------------------------------- MBFITSreader::getFreqInfo

// Get frequency parameters for each IF.

int MBFITSreader::getFreqInfo(
        int     &nIF,
        double* &startFreq,
        double* &endFreq)
{
  // This is RPFITS - can't do it!
  return 1;
}

//---------------------------------------------------- MBFITSreader::findRange

// Find the range of the data selected in time and position.

int MBFITSreader::findRange(
        int    &nRow,
        int    &nSel,
        char   dateSpan[2][32],
        double utcSpan[2],
        double* &positions)
{
  // This is RPFITS - can't do it!
  return 1;
}

//--------------------------------------------------------- MBFITSreader::read

// Read the next data record (if you're feeling lucky).

int MBFITSreader::read(
        MBrecord &MBrec)
{
  int beamNo = -1;
  int haveData, pCode = 0, status;
  double raRate = 0.0, decRate = 0.0, paRate = 0.0;
  MBrecord *iMBuff = 0x0;

  if (!cMBopen) {
    logMsg("ERROR, an MBFITS file has not been opened.");
    return 1;
  }

  // Positions recorded in the input records usually do not coincide with the
  // midpoint of the integration and hence the input must be buffered so that
  // true positions may be interpolated.
  //
  // On the first call nBeamSel buffers of length nBin, are allocated and
  // filled, where nBin is the number of time bins.
  //
  // The input records for binned, single beam data with multiple simultaneous
  // IFs are ordered by IF within each integration rather than by bin number
  // and hence are not in time order.  No multibeam data exists with
  // nBin > 1 but the likelihood that the input records would be in beam/IF
  // order and the requirement that output records be in time order would
  // force an elaborate double-buffering system and we do not support it.
  //
  // Once all buffers are filled, the next record for each beam pertains to
  // the next integration and should contain new position information allowing
  // the proper position for each spectrum in the buffer to be interpolated.
  // The buffers are then flushed in time order.  For single beam data there
  // is only one buffer and reads from the MBFITS file are suspended while the
  // flush is in progress.  For multibeam data each buffer is of unit length
  // so the flush completes immediately and the new record takes its place.

  haveData = 0;
  while (!haveData) {
    int iBeamSel = -1, iIFSel = -1;

    if (!cFlushing) {
      if (cEOF) {
        return -1;
      }

      // Read the next record.
      pCode = 0;
      if ((status = rpget(0, cEOS)) == -1) {
        // EOF.
        cEOF = 1;
        cFlushing = 1;
        cFlushBin = 0;
        cFlushIF  = 0;

#ifdef PKSIO_DEBUG
        fprintf(stderr, "\nEnd-of-file detected, flushing last cycle.\n");
#endif

      } else if (status) {
        // IO error.
        return 1;

      } else {
        if (cFirst) {
          // First data; cBeamSel[] stores the buffer index for each beam.
          cNBeamSel = 0;
          cBeamSel = new int[cNBeam];

          for (int iBeam = 0; iBeam < cNBeam; iBeam++) {
            if (cBeams[iBeam]) {
              // Buffer offset for this beam.
              cBeamSel[iBeam] = cNBeamSel++;
            } else {
              // Signal that the beam is not selected.
              cBeamSel[iBeam] = -1;
            }
          }

          // Set up bookkeeping arrays for IFs.
          cIFSel   = new int[cNIF];
          cChanOff = new int[cNIF];
          cXpolOff = new int[cNIF];

          int simulIF = 0;
          int maxChan = 0;
          int maxXpol = 0;

          for (int iIF = 0; iIF < cNIF; iIF++) {
            if (cIFs[iIF]) {
              // Buffer index for each IF within each simultaneous set.
              cIFSel[iIF] = 0;

              // Array offsets for each IF within each simultaneous set.
              cChanOff[iIF] = 0;
              cXpolOff[iIF] = 0;

              // Look for earlier IFs in the same simultaneous set.
              for (int jIF = 0; jIF < iIF; jIF++) {
                if (!cIFs[jIF]) continue;

                if (if_.if_simul[jIF] == if_.if_simul[iIF]) {
                  // Got one, increment indices.
                  cIFSel[iIF]++;

                  cChanOff[iIF] += cNChan[jIF] * cNPol[jIF];
                  if (cHaveXPol[jIF]) {
                    cXpolOff[iIF] += 2 * cNChan[jIF];
                  }
                }
              }

              // Maximum number of selected IFs in any simultaneous set.
              simulIF = max(simulIF, cIFSel[iIF]+1);

              // Maximum memory required for any simultaneous set.
              maxChan = max(maxChan, cChanOff[iIF] + cNChan[iIF]*cNPol[iIF]);
              if (cHaveXPol[iIF]) {
                maxXpol = max(maxXpol, cXpolOff[iIF] + 2*cNChan[iIF]);
              }

            } else {
              // Signal that the IF is not selected.
              cIFSel[iIF] = -1;
            }
          }

          // Check for binning mode observations.
          if (param_.intbase > 0.0f) {
            cNBin = int((cIntTime / param_.intbase) + 0.5);

            // intbase sometimes contains rubbish.
            if (cNBin == 0) {
              cNBin = 1;
            }
          } else {
            cNBin = 1;
          }

          if (cNBin > 1 && cNBeamSel > 1) {
            logMsg("ERROR, cannot handle binning mode for multiple beams.");
            close();
            return 1;
          }

          // Allocate buffer data storage; the MBrecord constructor zeroes
          // class members such as cycleNo that are tested in the first pass
          // below.
          int nBuff = cNBeamSel * cNBin;
          cBuffer = new MBrecord[nBuff];

          // Allocate memory for spectral arrays.
          for (int ibuff = 0; ibuff < nBuff; ibuff++) {
            cBuffer[ibuff].setNIFs(simulIF);
            cBuffer[ibuff].allocate(0, maxChan, maxXpol);

            // Signal that this IF in this buffer has been flushed.
            for (int iIF = 0; iIF < simulIF; iIF++) {
              cBuffer[ibuff].IFno[iIF] = 0;
            }
          }

          cPosUTC = new double[cNBeamSel];

          cFirst = 0;
          cScanNo  = 1;
          cCycleNo = 0;
          cPrevUTC = -1.0;
        }

        // Check for end-of-scan.
        if (cEOS) {
          cScanNo++;
          cCycleNo = 0;
          cPrevUTC = -1.0;
        }

        // Check for change-of-day.
        double cod = 0.0;
        if ((cUTC + 86400.0) < (cPrevUTC + 600.0)) {
          // cUTC should continue to increase past 86400 during a single scan.
          // However, if the RPFITS file contains multiple scans that straddle
          // midnight then cUTC can jump backwards from the end of one scan to
          // the start of the next.
#ifdef PKSIO_DEBUG
          fprintf(stderr, "Change-of-day on cUTC: %.1f -> %.1f",
            cPrevUTC, cUTC);
#endif
          // Can't change the recorded value of cUTC directly (without also
          // changing dateobs) so change-of-day must be recorded separately as
          // an offset to be applied when comparing integration timestamps.
          cod = 86400.0;

        } else if (cUTC < cPrevUTC - 1.0) {
          sprintf(cMsg,
            "WARNING: Cycle %d:%03d-%03d, UTC went backwards from\n"
            "         %.1f to %.1f!  Incrementing day number,\n"
            "         positions may be unreliable.", cScanNo, cCycleNo,
            cCycleNo+1, cPrevUTC, cUTC);
          logMsg(cMsg);
          cUTC += 86400.0;
        }

        if (cNBin > 1) {
          // Binning mode: correct the time.
          cUTC += param_.intbase * (cBin - (cNBin + 1)/2.0);
        }

        // New integration cycle?
        if ((cUTC+cod) > cPrevUTC) {
          cCycleNo++;
          cPrevUTC = cUTC + 0.0001;
        }

        // Apply beam selection.
        beamNo = int(cBaseline / 256.0);
        if (beamNo == 1) {
          // Store the position of beam 1 for grid convergence corrections.
          cRA0  = cU;
          cDec0 = cV;
        }
        iBeamSel = cBeamSel[beamNo-1];
        if (iBeamSel < 0) continue;

        // Sanity check (mainly for MOPS).
        if (cIFno > cNIF) continue;

        // Apply IF selection.
        iIFSel = cIFSel[cIFno - 1];
        if (iIFSel < 0) continue;

        sprintf(cDateObs, "%-10.10s", names_.datobs);
        cDateObs[10] = '\0';

        // Compute buffer number.
        iMBuff = cBuffer + iBeamSel;
        if (cNBin > 1) iMBuff += cNBeamSel*(cBin-1);

        if (cCycleNo < iMBuff->cycleNo) {
          // Note that if the first beam and IF are not both selected cEOS
          // will be cleared by rpget() when the next beam/IF is read.
          cEOS = 1;
        }

        // Begin flush cycle?
        if (cEOS || (iMBuff->nIF && (cUTC+cod) > (iMBuff->utc+0.0001))) {
          cFlushing = 1;
          cFlushBin = 0;
          cFlushIF  = 0;
        }

#ifdef PKSIO_DEBUG
        char rel = '=';
        double dt = utcDiff(cUTC, cW);
        if (dt < 0.0) {
          rel = '<';
        } else if (dt > 0.0) {
          rel = '>';
        }

        fprintf(stderr, "\n In:%4d%4d%3d%3d  %.3f %c %.3f (%+.3fs) - "
          "%sflushing\n", cScanNo, cCycleNo, beamNo, cIFno, cUTC, rel, cW, dt,
          cFlushing ? "" : "not ");
        if (cEOS) {
          fprintf(stderr, "Start of new scan, flushing previous scan.\n");
        }
#endif
      }
    }


    if (cFlushing) {
      // Find the oldest integration to flush, noting that the last
      // integration cycle may be incomplete.
      beamNo = 0;
      int cycleNo = 0;
      for (; cFlushBin < cNBin; cFlushBin++) {
        for (iBeamSel = 0; iBeamSel < cNBeamSel; iBeamSel++) {
          iMBuff = cBuffer + iBeamSel + cNBeamSel*cFlushBin;

          // iMBuff->nIF is set to zero (below) to signal that all IFs in
          // an integration have been flushed.
          if (iMBuff->nIF) {
            if (cycleNo == 0 || iMBuff->cycleNo < cycleNo) {
              beamNo  = iMBuff->beamNo;
              cycleNo = iMBuff->cycleNo;
            }
          }
        }

        if (beamNo) {
          // Found an integration to flush.
          break;
        }
      }

      if (beamNo) {
        iBeamSel = cBeamSel[beamNo-1];
        iMBuff = cBuffer + iBeamSel + cNBeamSel*cFlushBin;

        // Find the IF to flush.
        for (; cFlushIF < iMBuff->nIF; cFlushIF++) {
          if (iMBuff->IFno[cFlushIF]) break;
        }

      } else {
        // Flush complete.
        cFlushing = 0;
        if (cEOF) {
          return -1;
        }

        // The last record read must have been the first of a new cycle.
        beamNo = int(cBaseline / 256.0);
        iBeamSel = cBeamSel[beamNo-1];

        // Compute buffer number.
        iMBuff = cBuffer + iBeamSel;
        if (cNBin > 1) iMBuff += cNBeamSel*(cBin-1);
      }
    }


    if (cFlushing && cFlushBin == 0 && cFlushIF == 0 && cInterp) {
      // Start of flush cycle, interpolate the beam position.
      //
      // The position is measured by the control system at a time returned by
      // RPFITSIN as the 'w' visibility coordinate.  The ra and dec, returned
      // as the 'u' and 'v' visibility coordinates, must be interpolated to
      // the integration time which RPFITSIN returns as 'cUTC', this usually
      // being a second or two later.  The interpolation method used here is
      // based on the scan rate.
      //
      // "This" RA, Dec, and UTC refers to the position currently stored in
      // the buffer marked for output (iMBuff).  This position is interpolated
      // to the midpoint of that integration using either
      //   a) the rate currently sitting in iMBuff, which was computed from
      //      the previous integration, otherwise
      //   b) from the position recorded in the "next" integration which is
      //      currently sitting in the RPFITS commons,
      // so that the position timestamps straddle the midpoint of the
      // integration and is thereby interpolated rather than extrapolated.
      //
      // At the end of a scan, or if the next position has not been updated
      // or its timestamp does not advance sufficiently, the most recent
      // determination of the scan rate will be used for extrapolation which
      // is quantified by the "rate age" measured in seconds beyond the
      // interval defined by the position timestamps.

      // At this point, iMBuff contains cU, cV, cW, parAngle and focusRot
      // stored from the previous call to rpget() for this beam (i.e. "this"),
      // and also raRate, decRate and paRate computed from that integration
      // and the previous one.
      double thisRA  = iMBuff->ra;
      double thisDec = iMBuff->dec;
      double thisUTC = cPosUTC[iBeamSel];
      double thisPA  = iMBuff->parAngle + iMBuff->focusRot;

#ifdef PKSIO_DEBUG
      fprintf(stderr, "This (%d) ra, dec, UTC: %9.4f %9.4f %10.3f %9.4f\n",
        iMBuff->cycleNo, thisRA*R2D, thisDec*R2D, thisUTC, thisPA*R2D);
#endif

      if (cEOF || cEOS) {
        // Use rates from the last cycle.
        raRate  = iMBuff->raRate;
        decRate = iMBuff->decRate;
        paRate  = iMBuff->paRate;

      } else {
        if (cW == thisUTC) {
          // The control system at Mopra typically does not update the
          // positions between successive integration cycles at the end of a
          // scan (nor are they flagged).  In this case we use the previously
          // computed rates, even if from the previous scan since these are
          // likely to be a better guess than anything else.
          raRate  = iMBuff->raRate;
          decRate = iMBuff->decRate;
          paRate  = iMBuff->paRate;

          if (cU == thisRA && cV == thisDec) {
            // Position and timestamp unchanged.
            pCode = 1;

          } else if (fabs(cU-thisRA) < 0.0001 && fabs(cV-thisDec) < 0.0001) {
            // Allow small rounding errors (seen infrequently).
            pCode = 1;

          } else {
            // (cU,cV) are probably rubbish (not yet seen in practice).
            pCode = 2;
            cU = thisRA;
            cV = thisDec;
          }

#ifdef PKSIO_DEBUG
          fprintf(stderr, "Next (%d) ra, dec, UTC: %9.4f %9.4f %10.3f "
            "(0.000s)\n", cCycleNo, cU*R2D, cV*R2D, cW);
#endif

        } else {
          double nextRA  = cU;
          double nextDec = cV;

          // Check and, if necessary, repair the position timestamp,
          // remembering that pCode refers to the NEXT cycle.
          pCode = fixw(cDateObs, cCycleNo, beamNo, cAvRate, thisRA, thisDec,
                       thisUTC, nextRA, nextDec, cW);
          if (pCode > 0) pCode += 3;
          double nextUTC = cW;

#ifdef PKSIO_DEBUG
          fprintf(stderr, "Next (%d) ra, dec, UTC: %9.4f %9.4f %10.3f "
            "(%+.3fs)\n", cCycleNo, nextRA*R2D, nextDec*R2D, nextUTC,
            utcDiff(nextUTC, thisUTC));
#endif

          // Compute the scan rate for this beam.
          double dUTC = utcDiff(nextUTC, thisUTC);
          if ((0.0 < dUTC) && (dUTC < 600.0)) {
            scanRate(cRA0, cDec0, thisRA, thisDec, nextRA, nextDec, dUTC,
                     raRate, decRate);

            // Update the mean scan rate.
            cAvRate[0] = (cAvRate[0]*cNRate +  raRate) / (cNRate + 1);
            cAvRate[1] = (cAvRate[1]*cNRate + decRate) / (cNRate + 1);
            cNRate++;

            // Rate of change of position angle.
            if (sc_.sc_ant <= anten_.nant) {
              paRate = 0.0;
            } else {
              int iOff = sc_.sc_q * (sc_.sc_ant - 1) - 1;
              double nextPA = sc_.sc_cal[iOff + 4] + sc_.sc_cal[iOff + 7];
              double paDiff = nextPA - thisPA;
              if (paDiff > PI) {
                paDiff -= TWOPI;
              } else if (paDiff < -PI) {
                paDiff += TWOPI;
              }
              paRate = paDiff / dUTC;
            }

            if (cInterp == 2) {
              // Use the same interpolation scheme as the original pksmbfits
              // client.  This incorrectly assumed that (nextUTC - thisUTC) is
              // equal to the integration time and interpolated by computing a
              // weighted sum of the positions before and after the required
              // time.

              double utc = iMBuff->utc;
              double tw1 = 1.0 - utcDiff(utc, thisUTC) / iMBuff->exposure;
              double tw2 = 1.0 - utcDiff(nextUTC, utc) / iMBuff->exposure;
              double gamma = (tw2 / (tw1 + tw2)) * dUTC / (utc - thisUTC);

              // Guard against RA cycling through 24h in either direction.
              if (fabs(nextRA - thisRA) > PI) {
                if (nextRA < thisRA) {
                  nextRA += TWOPI;
                } else {
                  nextRA -= TWOPI;
                }
              }

              raRate  = gamma * (nextRA  - thisRA)  / dUTC;
              decRate = gamma * (nextDec - thisDec) / dUTC;
            }

          } else {
            if (cCycleNo == 2 && fabs(utcDiff(cUTC,cW)) < 600.0) {
              // thisUTC (i.e. cW for the first cycle) is rubbish, and
              // probably the position as well (extremely rare in practice,
              // e.g. 97-12-19_1029_235708-18_586e.hpf which actually has the
              // t/1000 scaling bug in the first cycle).
              iMBuff->pCode = 3;
              thisRA  = cU;
              thisDec = cV;
              thisUTC = cW;
              raRate  = 0.0;
              decRate = 0.0;
              paRate  = 0.0;

            } else {
              // cW is rubbish and probably (cU,cV), and possibly the
              // parallactic angle and everything else as well (rarely seen
              // in practice, e.g. 97-12-09_0743_235707-58_327c.hpf and
              // 97-09-01_0034_123717-42_242b.hpf, the latter with bad
              // parallactic angle).
              pCode = 3;
              cU = thisRA;
              cV = thisDec;
              cW = thisUTC;
              raRate  = iMBuff->raRate;
              decRate = iMBuff->decRate;
              paRate  = iMBuff->paRate;
            }
          }
        }
      }


      // Choose the closest rate determination.
      if (cCycleNo == 1) {
        // Scan containing a single integration.
        iMBuff->raRate  = 0.0;
        iMBuff->decRate = 0.0;
        iMBuff->paRate  = 0.0;

      } else {
        double dUTC = iMBuff->utc - cPosUTC[iBeamSel];

        if (dUTC >= 0.0) {
          // In HIPASS/ZOA, the position timestamp, which should always occur
          // on the whole second, normally precedes an integration midpoint
          // falling on the half-second.  Consequently, positive ages are
          // always half-integral.
          dUTC = utcDiff(iMBuff->utc, cW);
          if (dUTC > 0.0) {
            iMBuff->rateAge = dUTC;
          } else {
            iMBuff->rateAge = 0.0f;
          }

          iMBuff->raRate  =  raRate;
          iMBuff->decRate = decRate;
          iMBuff->paRate  =  paRate;

        } else {
          // In HIPASS/ZOA, negative ages occur when the integration midpoint,
          // occurring on the whole second, precedes the position timestamp.
          // Thus negative ages are always an integral number of seconds.
          // They have only been seen to occur sporadically in the period
          // 1999/05/31 to 1999/11/01, e.g. 1999-07-26_1821_005410-74_007c.hpf
          //
          // In recent (2008/10/07) Mopra data, small negative ages (~10ms,
          // occasionally up to ~300ms) seem to be the norm, with both the
          // position timestamp and integration midpoint falling close to but
          // not on the integral second.
          if (cCycleNo == 2) {
            // We have to start with something!
            iMBuff->rateAge = dUTC;

          } else {
            // Although we did not record the relevant position timestamp
            // explicitly, it can easily be deduced.
            double w = iMBuff->utc - utcDiff(cUTC, iMBuff->utc) -
                       iMBuff->rateAge;
            dUTC = utcDiff(iMBuff->utc, w);

            if (dUTC > 0.0) {
              iMBuff->rateAge = 0.0f;
            } else {
              iMBuff->rateAge = dUTC;
            }
          }

          iMBuff->raRate  =  raRate;
          iMBuff->decRate = decRate;
          iMBuff->paRate  =  paRate;
        }
      }

#ifdef PKSIO_DEBUG
      double avRate = sqrt(cAvRate[0]*cAvRate[0] + cAvRate[1]*cAvRate[1]);
      fprintf(stderr, "RA, Dec, Av & PA rates: %8.4f %8.4f %8.4f %8.4f "
        "pCode %d\n", raRate*R2D, decRate*R2D, avRate*R2D, paRate*R2D, pCode);
#endif


      // Compute the position of this beam for all bins.
      for (int idx = 0; idx < cNBin; idx++) {
        int jbuff = iBeamSel + cNBeamSel*idx;

        cBuffer[jbuff].raRate  = iMBuff->raRate;
        cBuffer[jbuff].decRate = iMBuff->decRate;
        cBuffer[jbuff].paRate  = iMBuff->paRate;

        double dUTC = utcDiff(cBuffer[jbuff].utc, thisUTC);
        if (dUTC > 100.0) {
          // Must have cycled through midnight.
          dUTC -= 86400.0;
        }

        applyRate(cRA0, cDec0, thisRA, thisDec,
          cBuffer[jbuff].raRate, cBuffer[jbuff].decRate, dUTC,
          cBuffer[jbuff].ra, cBuffer[jbuff].dec);

#ifdef PKSIO_DEBUG
        fprintf(stderr, "Intp (%d) ra, dec, UTC: %9.4f %9.4f %10.3f (pCode, "
          "age: %d %.1fs)\n", iMBuff->cycleNo, cBuffer[jbuff].ra*R2D,
          cBuffer[jbuff].dec*R2D, cBuffer[jbuff].utc, iMBuff->pCode,
          iMBuff->rateAge);
#endif
      }
    }


    if (cFlushing) {
      // Copy buffer location out one IF at a time.
      MBrec.extract(*iMBuff, cFlushIF);
      haveData = 1;

#ifdef PKSIO_DEBUG
      fprintf(stderr, "Out:%4d%4d%3d%3d\n", MBrec.scanNo, MBrec.cycleNo,
        MBrec.beamNo, MBrec.IFno[0]);
#endif

      // Signal that this IF in this buffer location has been flushed.
      iMBuff->IFno[cFlushIF] = 0;

      if (cFlushIF == iMBuff->nIF - 1) {
        // Signal that all IFs in this buffer location have been flushed.
        iMBuff->nIF = 0;

        // Stop cEOS being set when the next integration is read.
        iMBuff->cycleNo = 0;

      } else {
        // Carry on flushing the other IFs.
        continue;
      }

      // Has the whole buffer been flushed?
      if (cFlushBin == cNBin - 1) {
        if (cEOS || cEOF) {
          // Carry on flushing other buffers.
          cFlushIF = 0;
          continue;
        }

        cFlushing = 0;

        beamNo = int(cBaseline / 256.0);
        iBeamSel = cBeamSel[beamNo-1];

        // Compute buffer number.
        iMBuff = cBuffer + iBeamSel;
        if (cNBin > 1) iMBuff += cNBeamSel*(cBin-1);
      }
    }

    if (!cFlushing) {
      // Buffer this MBrec.
      if ((cScanNo > iMBuff->scanNo) && iMBuff->IFno[0]) {
        // Sanity check on the number of IFs in the new scan.
        if (if_.n_if != cNIF) {
          sprintf(cMsg, "WARNING: Scan %d has %d IFs instead of %d, "
            "continuing.", cScanNo, if_.n_if, cNIF);
          logMsg(cMsg);
        }
      }

      // Sanity check on incomplete integrations within a scan.
      if (iMBuff->nIF && (iMBuff->cycleNo != cCycleNo)) {
        // Force the incomplete integration to be flushed before proceeding.
        cFlushing = 1;
        continue;
      }

#ifdef PKSIO_DEBUG
      fprintf(stderr, "Buf:%4d%4d%3d%3d\n", cScanNo, cCycleNo, beamNo, cIFno);
#endif

      // Store IF-independent parameters only for the first IF of a new cycle,
      // particularly because this is the only one for which the scan rates
      // are computed above.
      int firstIF = (iMBuff->nIF == 0);
      if (firstIF) {
        iMBuff->scanNo  = cScanNo;
        iMBuff->cycleNo = cCycleNo;

        // Times.
        strcpy(iMBuff->datobs, cDateObs);
        iMBuff->utc = cUTC;
        iMBuff->exposure = param_.intbase;

        // Source identification.
        sprintf(iMBuff->srcName, "%-16.16s",
                names_.su_name + (cSrcNo-1)*16);
        iMBuff->srcName[16] = '\0';
        iMBuff->srcRA  = doubles_.su_ra[cSrcNo-1];
        iMBuff->srcDec = doubles_.su_dec[cSrcNo-1];

        // Rest frequency of the line of interest.
        iMBuff->restFreq = doubles_.rfreq;
        if (strncmp(names_.instrument, "ATPKSMB", 7) == 0) {
          // Fix the HI rest frequency recorded for Parkes multibeam data.
          double reffreq  = doubles_.freq;
          double restfreq = doubles_.rfreq;
          if ((restfreq == 0.0 || fabs(restfreq - reffreq) == 0.0) &&
               fabs(reffreq - 1420.405752e6) < 100.0) {
            iMBuff->restFreq = 1420.405752e6;
          }
        }

        // Observation type.
        int j;
        for (j = 0; j < 15; j++) {
          iMBuff->obsType[j] = names_.card[11+j];
          if (iMBuff->obsType[j] == '\'') break;
        }
        iMBuff->obsType[j] = '\0';

        // Beam-dependent parameters.
        iMBuff->beamNo = beamNo;

        // Beam position at the specified time.
        if (cSUpos) {
          // Non-ATNF data that does not store the position in (u,v,w).
          iMBuff->ra  = doubles_.su_ra[cSrcNo-1];
          iMBuff->dec = doubles_.su_dec[cSrcNo-1];
        } else {
          iMBuff->ra  = cU;
          iMBuff->dec = cV;
        }
        cPosUTC[iBeamSel] = cW;
        iMBuff->pCode = pCode;

        // Store rates for next time.
        iMBuff->raRate  =  raRate;
        iMBuff->decRate = decRate;
        iMBuff->paRate  =  paRate;
      }

      // IF-dependent parameters.
      int iIF = cIFno - 1;
      int startChan = cStartChan[iIF];
      int endChan   = cEndChan[iIF];
      int refChan   = cRefChan[iIF];

      int nChan = abs(endChan - startChan) + 1;

      iIFSel = cIFSel[iIF];
      if (iMBuff->IFno[iIFSel] == 0) {
        iMBuff->nIF++;
        iMBuff->IFno[iIFSel] = cIFno;
      } else {
        // Integration cycle written to the output file twice (the only known
        // example is 1999-05-22_1914_000-031805_03v.hpf).
        sprintf(cMsg, "WARNING: Integration cycle %d:%d, beam %2d, \n"
                      "         IF %d was duplicated.", cScanNo, cCycleNo-1,
                      beamNo, cIFno);
        logMsg(cMsg);
      }
      iMBuff->nChan[iIFSel] = nChan;
      iMBuff->nPol[iIFSel]  = cNPol[iIF];

      iMBuff->fqRefPix[iIFSel] = doubles_.if_ref[iIF];
      iMBuff->fqRefVal[iIFSel] = doubles_.if_freq[iIF];
      iMBuff->fqDelt[iIFSel]   =
        if_.if_invert[iIF] * fabs(doubles_.if_bw[iIF] /
          (if_.if_nfreq[iIF] - 1));

      // Adjust for channel selection.
      if (iMBuff->fqRefPix[iIFSel] != refChan) {
        iMBuff->fqRefVal[iIFSel] +=
          (refChan - iMBuff->fqRefPix[iIFSel]) *
            iMBuff->fqDelt[iIFSel];
        iMBuff->fqRefPix[iIFSel] = refChan;
      }

      if (endChan < startChan) {
        iMBuff->fqDelt[iIFSel] = -iMBuff->fqDelt[iIFSel];
      }


      // System temperature.
      int iBeam = beamNo - 1;
      int scq = sc_.sc_q;
      float TsysPol1 = sc_.sc_cal[scq*iBeam + 3];
      float TsysPol2 = sc_.sc_cal[scq*iBeam + 4];
      iMBuff->tsys[iIFSel][0] = TsysPol1*TsysPol1;
      iMBuff->tsys[iIFSel][1] = TsysPol2*TsysPol2;

      // Calibration factor; may be changed later if the data is recalibrated.
      if (scq > 14) {
        // Will only be present for Parkes Multibeam or LBA data.
        iMBuff->calfctr[iIFSel][0] = sc_.sc_cal[scq*iBeam + 14];
        iMBuff->calfctr[iIFSel][1] = sc_.sc_cal[scq*iBeam + 15];
      } else {
        iMBuff->calfctr[iIFSel][0] = 0.0f;
        iMBuff->calfctr[iIFSel][1] = 0.0f;
      }

      // Cross-polarization calibration factor (unknown to MBFITS).
      for (int j = 0; j < 2; j++) {
        iMBuff->xcalfctr[iIFSel][j] = 0.0f;
      }

      // Baseline parameters (unknown to MBFITS).
      iMBuff->haveBase = 0;

      // Data (always present in MBFITS).
      iMBuff->haveSpectra = 1;

      // Flag:  bit 0 set if off source.
      //        bit 1 set if loss of sync in A polarization.
      //        bit 2 set if loss of sync in B polarization.
      unsigned char rpflag =
        (unsigned char)(sc_.sc_cal[scq*iBeam + 12] + 0.5f);

      // The baseline flag may be set independently.
      if (rpflag == 0) rpflag = cFlag;

      // Copy and scale data.
      int inc = 2 * if_.if_nstok[iIF];
      if (endChan < startChan) inc = -inc;

      float TsysF;
      iMBuff->spectra[iIFSel] = iMBuff->spectra[0] + cChanOff[iIF];
      iMBuff->flagged[iIFSel] = iMBuff->flagged[0] + cChanOff[iIF];

      float *spectra = iMBuff->spectra[iIFSel];
      unsigned char *flagged = iMBuff->flagged[iIFSel];
      for (int ipol = 0; ipol < cNPol[iIF]; ipol++) {
        if (sc_.sc_cal[scq*iBeam + 3 + ipol] > 0.0f) {
          // The correlator has already applied the calibration.
          TsysF = 1.0f;
        } else {
          // The correlator has normalized cVis[k] to a Tsys of 500K.
          TsysF = iMBuff->tsys[iIFSel][ipol] / 500.0f;
        }

        int k = 2 * (if_.if_nstok[iIF]*(startChan - 1) + ipol);
        for (int ichan = 0; ichan < nChan; ichan++) {
          *(spectra++) = TsysF * cVis[k];
          *(flagged++) = rpflag;
          k += inc;
        }
      }

      if (cHaveXPol[iIF]) {
        int k = 2 * (3*(startChan - 1) + 2);
        iMBuff->xpol[iIFSel] = iMBuff->xpol[0] + cXpolOff[iIF];
        float *xpol = iMBuff->xpol[iIFSel];
        for (int ichan = 0; ichan < nChan; ichan++) {
          *(xpol++) = cVis[k];
          *(xpol++) = cVis[k+1];
          k += inc;
        }
      }


      // Calibration factor applied to the data by the correlator.
      if (scq > 14) {
        // Will only be present for Parkes Multibeam or LBA data.
        iMBuff->tcal[iIFSel][0] = sc_.sc_cal[scq*iBeam + 14];
        iMBuff->tcal[iIFSel][1] = sc_.sc_cal[scq*iBeam + 15];
      } else {
        iMBuff->tcal[iIFSel][0] = 0.0f;
        iMBuff->tcal[iIFSel][1] = 0.0f;
      }

      if (firstIF) {
        if (sc_.sc_ant <= anten_.nant) {
          // No extra syscal information present.
          iMBuff->extraSysCal = 0;
          iMBuff->azimuth   = 0.0f;
          iMBuff->elevation = 0.0f;
          iMBuff->parAngle  = 0.0f;
          iMBuff->focusAxi  = 0.0f;
          iMBuff->focusTan  = 0.0f;
          iMBuff->focusRot  = 0.0f;
          iMBuff->temp      = 0.0f;
          iMBuff->pressure  = 0.0f;
          iMBuff->humidity  = 0.0f;
          iMBuff->windSpeed = 0.0f;
          iMBuff->windAz    = 0.0f;
          strcpy(iMBuff->tcalTime, "                ");
          iMBuff->refBeam = 0;

        } else {
          // Additional information for Parkes Multibeam data.
          int iOff = scq*(sc_.sc_ant - 1) - 1;
          iMBuff->extraSysCal = 1;

          iMBuff->azimuth   = sc_.sc_cal[iOff + 2];
          iMBuff->elevation = sc_.sc_cal[iOff + 3];
          iMBuff->parAngle  = sc_.sc_cal[iOff + 4];

          iMBuff->focusAxi  = sc_.sc_cal[iOff + 5] * 1e-3;
          iMBuff->focusTan  = sc_.sc_cal[iOff + 6] * 1e-3;
          iMBuff->focusRot  = sc_.sc_cal[iOff + 7];

          iMBuff->temp      = sc_.sc_cal[iOff + 8];
          iMBuff->pressure  = sc_.sc_cal[iOff + 9];
          iMBuff->humidity  = sc_.sc_cal[iOff + 10];
          iMBuff->windSpeed = sc_.sc_cal[iOff + 11];
          iMBuff->windAz    = sc_.sc_cal[iOff + 12];

          char *tcalTime = iMBuff->tcalTime;
          sprintf(tcalTime, "%-16.16s", (char *)(&sc_.sc_cal[iOff+13]));
          tcalTime[16] = '\0';

#ifndef AIPS_LITTLE_ENDIAN
          // Do byte swapping on the ASCII date string.
          for (int j = 0; j < 16; j += 4) {
            char ctmp;
            ctmp = tcalTime[j];
            tcalTime[j]   = tcalTime[j+3];
            tcalTime[j+3] = ctmp;
            ctmp = tcalTime[j+1];
            tcalTime[j+1] = tcalTime[j+2];
            tcalTime[j+2] = ctmp;
          }
#endif

          // Reference beam number.
          float refbeam = sc_.sc_cal[iOff + 17];
          if (refbeam > 0.0f || refbeam < 100.0f) {
            iMBuff->refBeam = int(refbeam);
          } else {
            iMBuff->refBeam = 0;
          }
        }
      }
    }
  }

  return 0;
}

//-------------------------------------------------------- MBFITSreader::rpget

// Read the next data record from the RPFITS file.

int MBFITSreader::rpget(int syscalonly, int &EOS)
{
  EOS = 0;

  int retries = 0;

  // Allow 10 read errors.
  int numErr = 0;

  int jstat = 0;
  while (numErr < 10) {
    int lastjstat = jstat;

    switch(rpfitsin(jstat)) {
    case -1:
      // Read failed; retry.
      numErr++;
      logMsg("WARNING: RPFITS read failed - retrying.");
      jstat = 0;
      break;

    case 0:
      // Successful read.
      if (lastjstat == 0) {
        if (cBaseline == -1) {
          // Syscal data.
          if (syscalonly) {
            return 0;
          }

        } else {
          if (!syscalonly) {
            return 0;
          }
        }
      }

      // Last operation was to read header or FG table; now read data.
      break;

    case 1:
      // Encountered header while trying to read data; read it.
      EOS = 1;
      jstat = -1;
      break;

    case 2:
      // End of scan; read past it.
      jstat = 0;
      break;

    case 3:
      // End-of-file; retry applies to real-time mode.
      if (retries++ >= cRetry) {
        return -1;
      }

      sleep(10);
      jstat = 0;
      break;

    case 4:
      // Encountered FG table while trying to read data; read it.
      jstat = -1;
      break;

    case 5:
      // Illegal data at end of block after close/reopen operation; retry.
      jstat = 0;
      break;

    default:
      // Shouldn't reach here.
      sprintf(cMsg, "WARNING: Unrecognized RPFITSIN return code: %d "
                    "(retrying).", jstat);
      logMsg(cMsg);
      jstat = 0;
      break;
    }
  }

  logMsg("ERROR: RPFITS read failed too many times.");
  return 2;
}

//----------------------------------------------------- MBFITSreader::rpfitsin

// Wrapper around RPFITSIN that reports errors.  Returned RPFITSIN subroutine
// arguments are captured as MBFITSreader member variables.

int MBFITSreader::rpfitsin(int &jstat)

{
  rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
            &cBin, &cIFno, &cSrcNo);

  // Handle messages from RPFITSIN.
  if (names_.errmsg[0] != ' ') {
    int i;
    for (i = 80; i > 0; i--) {
      if (names_.errmsg[i-1] != ' ') break;
    }

    sprintf(cMsg, "WARNING: Cycle %d:%03d, RPFITSIN reported -\n"
                  "         %.*s", cScanNo, cCycleNo, i, names_.errmsg);
    logMsg(cMsg);
  }

  return jstat;
}

//------------------------------------------------------- MBFITSreader::fixPos

// Check and, if necessary, repair a position timestamp.
//
// Problems with the position timestamp manifest themselves via the scan rate:
//
//   1) Zero scan rate pairs, 1997/02/28 to 1998/01/07
//
//      These occur because the position timestamp for the first integration
//      of the pair is erroneous; the value recorded is t/1000, where t is the
//      true value.
//        Earliest known: 97-02-28_1725_132653-42_258a.hpf
//          Latest known: 98-01-02_1923_095644-50_165c.hpf
//        (time range chosen to encompass observing runs).
//
//   2) Slow-fast scan rate pairs (0.013 - 0.020 deg/s),
//        1997/03/28 to 1998/01/07.
//
//      The UTC position timestamp is 1.0s later than it should be (never
//      earlier), almost certainly arising from an error in the telescope
//      control system.
//        Earliest known: 97-03-28_0150_010420-74_008d.hpf
//          Latest known: 98-01-04_1502_065150-02_177c.hpf
//        (time range chosen to encompass observing runs).
//
//   3) Slow-fast scan rate pairs (0.015 - 0.018 deg/s),
//        1999/05/20 to 2001/07/12 (HIPASS and ZOA),
//        2001/09/02 to 2001/12/04 (HIPASS and ZOA),
//        2002/03/28 to 2002/05/13 (ZOA only),
//        2003/04/26 to 2003/06/09 (ZOA only).
//        Earliest known: 1999-05-20_1818_175720-50_297e.hpf
//          Latest known: 2001-12-04_1814_065531p14_173e.hpf (HIPASS)
//                        2003-06-09_1924_352-085940_-6c.hpf (ZOA)
//
//      Caused by the Linux signalling NaN problem.  IEEE "signalling" NaNs
//      are silently transformed to "quiet" NaNs during assignment by setting
//      bit 22.  This affected RPFITS because of its use of VAX-format
//      floating-point numbers which, with their permuted bytes, may sometimes
//      appear as signalling NaNs.
//
//      The problem arose when the linux correlator came online and was
//      fixed with a workaround to the RPFITS library (repeated episodes
//      are probably due to use of an older version of the library).  It
//      should not have affected the data significantly because of the
//      low relative error, which ranges from 0.0000038 to 0.0000076, but
//      it is important for the computation of scan rates which requires
//      taking the difference of two large UTC timestamps, one or other
//      of which will have 0.5s added to it.
//
// The return value identifies which, if any, of these problems was repaired.

int MBFITSreader::fixw(
  const char *datobs,
  int    cycleNo,
  int    beamNo,
  double avRate[2],
  double thisRA,
  double thisDec,
  double thisUTC,
  double nextRA,
  double nextDec,
  float &nextUTC)
{
  if (strcmp(datobs, "2003-06-09") > 0) {
    return 0;

  } else if (strcmp(datobs, "1998-01-07") <= 0) {
    if (nextUTC < thisUTC && (nextUTC + 86400.0) > (thisUTC + 600.0)) {
      // Possible scaling problem.
      double diff = nextUTC*1000.0 - thisUTC;
      if (0.0 < diff && diff < 600.0) {
        nextUTC *= 1000.0;
        return 1;
      } else {
        // Irreparable.
        return -1;
      }
    }

    if (cycleNo > 2) {
      if (beamNo == 1) {
        // This test is only reliable for beam 1.
        double dUTC = nextUTC - thisUTC;
        if (dUTC < 0.0) dUTC += 86400.0;

        // Guard against RA cycling through 24h in either direction.
        if (fabs(nextRA - thisRA) > PI) {
          if (nextRA < thisRA) {
            nextRA += TWOPI;
          } else {
            nextRA -= TWOPI;
          }
        }

        double  dRA = (nextRA  - thisRA) * cos(nextDec);
        double dDec =  nextDec - thisDec;
        double  arc = sqrt(dRA*dRA + dDec*dDec);

        double averate = sqrt(avRate[0]*avRate[0] + avRate[1]*avRate[1]);
        double diff1 = fabs(averate - arc/(dUTC-1.0));
        double diff2 = fabs(averate - arc/dUTC);
        if ((diff1 < diff2) && (diff1 < 0.05*averate)) {
          nextUTC -= 1.0;
          cCode5 = cycleNo;
          return 2;
        } else {
          cCode5 = 0;
        }

      } else {
        if (cycleNo == cCode5) {
          nextUTC -= 1.0;
          return 2;
        }
      }
    }

  } else if ((strcmp(datobs, "1999-05-20") >= 0 &&
              strcmp(datobs, "2001-07-12") <= 0) ||
             (strcmp(datobs, "2001-09-02") >= 0 &&
              strcmp(datobs, "2001-12-04") <= 0) ||
             (strcmp(datobs, "2002-03-28") >= 0 &&
              strcmp(datobs, "2002-05-13") <= 0) ||
             (strcmp(datobs, "2003-04-26") >= 0 &&
              strcmp(datobs, "2003-06-09") <= 0)) {
    // Signalling NaN problem, e.g. 1999-07-26_1839_011106-74_009c.hpf.
    // Position timestamps should always be an integral number of seconds.
    double resid = nextUTC - int(nextUTC);
    if (resid == 0.5) {
      nextUTC -= 0.5;
      return 3;
    }
  }

  return 0;
}

//-------------------------------------------------------- MBFITSreader::close

// Close the input file.

void MBFITSreader::close(void)
{
  if (cMBopen) {
    int jstat = 1;
    rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
              &cBin, &cIFno, &cSrcNo);

    if (cBeams)     delete [] cBeams;
    if (cIFs)       delete [] cIFs;
    if (cNChan)     delete [] cNChan;
    if (cNPol)      delete [] cNPol;
    if (cHaveXPol)  delete [] cHaveXPol;
    if (cStartChan) delete [] cStartChan;
    if (cEndChan)   delete [] cEndChan;
    if (cRefChan)   delete [] cRefChan;

    if (cVis) delete [] cVis;
    if (cWgt) delete [] cWgt;

    if (cBeamSel)   delete [] cBeamSel;
    if (cIFSel)     delete [] cIFSel;
    if (cChanOff)   delete [] cChanOff;
    if (cXpolOff)   delete [] cXpolOff;
    if (cBuffer)    delete [] cBuffer;
    if (cPosUTC)    delete [] cPosUTC;

    cMBopen = 0;
  }
}

//-------------------------------------------------------------------- utcDiff

// Subtract two UTCs (s) allowing for any plausible number of cycles through
// 86400s, returning a result in the range [-43200, +43200]s.

double MBFITSreader::utcDiff(double utc1, double utc2)
{
  double diff = utc1 - utc2;

  if (diff > 43200.0) {
    diff -= 86400.0;
    while (diff > 43200.0) diff -= 86400.0;
  } else if (diff < -43200.0) {
    diff += 86400.0;
    while (diff < -43200.0) diff += 86400.0;
  }

  return diff;
}

//------------------------------------------------------- scanRate & applyRate

// Compute and apply the scan rate corrected for grid convergence.  (ra0,dec0)
// are the coordinates of the central beam, assumed to be the tracking centre.
// The rate computed in RA will be a rate of change of angular distance in the
// direction of increasing RA at the position of the central beam.  Similarly
// for declination.  Angles in radian, time in s.

void MBFITSreader::scanRate(
  double ra0,
  double dec0,
  double ra1,
  double dec1,
  double ra2,
  double dec2,
  double dt,
  double &raRate,
  double &decRate)
{
  // Transform to a system where the central beam lies on the equator at 12h.
  eulerx(ra1, dec1, ra0+HALFPI, -dec0, -HALFPI, ra1, dec1);
  eulerx(ra2, dec2, ra0+HALFPI, -dec0, -HALFPI, ra2, dec2);

  raRate  = (ra2  - ra1)  / dt;
  decRate = (dec2 - dec1) / dt;
}


void MBFITSreader::applyRate(
  double ra0,
  double dec0,
  double ra1,
  double dec1,
  double raRate,
  double decRate,
  double dt,
  double &ra2,
  double &dec2)
{
  // Transform to a system where the central beam lies on the equator at 12h.
  eulerx(ra1, dec1, ra0+HALFPI, -dec0, -HALFPI, ra1, dec1);

  ra2  = ra1  + (raRate  * dt);
  dec2 = dec1 + (decRate * dt);

  // Transform back.
  eulerx(ra2, dec2, -HALFPI, dec0, ra0+HALFPI, ra2, dec2);
}

//--------------------------------------------------------------------- eulerx

void MBFITSreader::eulerx(
  double lng0,
  double lat0,
  double phi0,
  double theta,
  double phi,
  double &lng1,
  double &lat1)

// Applies the Euler angle based transformation of spherical coordinates.
//
//     phi0  Longitude of the ascending node in the old system, radians.  The
//           ascending node is the point of intersection of the equators of
//           the two systems such that the equator of the new system crosses
//           from south to north as viewed in the old system.
//
//    theta  Angle between the poles of the two systems, radians.  THETA is
//           positive for a positive rotation about the ascending node.
//
//      phi  Longitude of the ascending node in the new system, radians.

{
  // Compute intermediaries.
  double lng0p  = lng0 - phi0;
  double slng0p = sin(lng0p);
  double clng0p = cos(lng0p);
  double slat0  = sin(lat0);
  double clat0  = cos(lat0);
  double ctheta = cos(theta);
  double stheta = sin(theta);

  double x = clat0*clng0p;
  double y = clat0*slng0p*ctheta + slat0*stheta;

  // Longitude in the new system.
  if (x != 0.0 || y != 0.0) {
    lng1 = phi + atan2(y, x);
  } else {
    // Longitude at the poles in the new system is consistent with that
    // specified in the old system.
    lng1 = phi + lng0p;
  }
  lng1 = fmod(lng1, TWOPI);
  if (lng1 < 0.0) lng1 += TWOPI;

  lat1 = asin(slat0*ctheta - clat0*stheta*slng0p);
}