Changeset 1757 for branches/alma/external/atnf/PKSIO/MBFITSreader.cc

Timestamp:

06/09/10 19:03:06 (14 years ago)

Author:

Kana Sugimoto

Message:

New Development: Yes

JIRA Issue: Yes (CAS-2211)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: ASAP 3.0.0 interface changes

Test Programs:

Put in Release Notes: Yes

Module(s): all the CASA sd tools and tasks are affected.

Description: Merged ATNF-ASAP 3.0.0 developments to CASA (alma) branch.

Note you also need to update casa/code/atnf.

Location:

branches/alma

Files:

• . (modified) (2 props)
• external/atnf/PKSIO/MBFITSreader.cc (modified) (54 diffs)

branches/alma/external/atnf/PKSIO/MBFITSreader.cc

@@ -2 +2 @@
 //# MBFITSreader.cc: ATNF single-dish RPFITS reader.
 //#---------------------------------------------------------------------------
-//# Copyright (C) 2000-2006
-//# Mark Calabretta, ATNF
+//# livedata - processing pipeline for single-dish, multibeam spectral data.
+//# Copyright (C) 2000-2009, Australia Telescope National Facility, CSIRO
 //#
-//# 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 file is part of livedata.
 //#
-//# This library is distributed in the hope that it will be useful, but WITHOUT
+//# livedata is free software: you can redistribute it and/or modify it under
+//# the terms of the GNU General Public License as published by the Free
+//# Software Foundation, either version 3 of the License, or (at your option)
+//# any later version.
+//#
+//# livedata 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.
+//# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 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.
+//# You should have received a copy of the GNU General Public License along
+//# with livedata.  If not, see <http://www.gnu.org/licenses/>.
 //#
-//# 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,
+//# Correspondence concerning livedata may be directed to:
+//#        Internet email: mcalabre@atnf.csiro.au
+//#        Postal address: Dr. Mark Calabretta
+//#                        Australia Telescope National Facility, CSIRO
+//#                        PO Box 76
+//#                        Epping NSW 1710
 //#                        AUSTRALIA
 //#
-//# $Id$
+//# http://www.atnf.csiro.au/computing/software/livedata.html
+//# $Id: MBFITSreader.cc,v 19.57 2009-10-30 06:34:36 cal103 Exp $
 //#---------------------------------------------------------------------------
 //# The MBFITSreader class reads single dish RPFITS files (such as Parkes
@@ -35 +37 @@
 //#---------------------------------------------------------------------------
 
+#include <atnf/pks/pks_maths.h>
 #include <atnf/PKSIO/MBFITSreader.h>
-#include <atnf/PKSIO/PKSMBrecord.h>
-
-#include <RPFITS.h>
+#include <atnf/PKSIO/MBrecord.h>
+
+#include <casa/Logging/LogIO.h>
 
 #include <casa/math.h>
@@ -47 +50 @@
 #include <unistd.h>
 
+#include <RPFITS.h>
+
 using namespace std;
 
@@ -52 +57 @@
 const double PI = 3.141592653589793238462643;
 const double TWOPI = 2.0 * PI;
+const double HALFPI = PI / 2.0;
+const double R2D = 180.0 / PI;
+
+// Class name
+const string className = "MBFITSreader" ;
 
 //------------------------------------------------- MBFITSreader::MBFITSreader
@@ -81 +91 @@
   cRefChan   = 0x0;
 
-  cVis = new float[2*4*8163];
+  cVis = 0x0;
+  cWgt = 0x0;
 
   cBeamSel   = 0x0;
@@ -91 +102 @@
 
   cMBopen = 0;
-  jstat = -3;
+
+  // Tell RPFITSIN not to report errors directly.
+  //iostat_.errlun = -1;
 }
 
@@ -120 +133 @@
         int  &extraSysCal)
 {
+  const string methodName = "open()" ;
+  LogIO os( LogOrigin( className, methodName, WHERE ) ) ;
+
   if (cMBopen) {
     close();
@@ -127 +143 @@
 
   // Open the RPFITS file.
-  rpfitsin_(&jstat, cVis, weight, &baseline, &ut, &u, &v, &w, &flag, &bin,
-            &if_no, &sourceno);
-
-  if (jstat) {
-    fprintf(stderr, "Failed to open MBFITS file: %s\n", rpname);
+  int jstat = -3;
+  if (rpfitsin(jstat)) {
+    sprintf(cMsg, "Failed to open MBFITS file\n%s", rpname);
+    os << LogIO::SEVERE << cMsg << LogIO::POST ;
     return 1;
   }
@@ -147 +162 @@
   // Read the first header.
   jstat = -1;
-  rpfitsin_(&jstat, cVis, weight, &baseline, &ut, &u, &v, &w, &flag, &bin,
-            &if_no, &sourceno);
-
-  if (jstat) {
-    fprintf(stderr, "Failed to read MBFITS header: %s\n", rpname);
+  if (rpfitsin(jstat)) {
+    sprintf(cMsg, "Failed to read MBFITS header in file\n"
+                  "%s", rpname);
+    os << LogIO::SEVERE << cMsg << LogIO::POST ;
     close();
     return 1;
@@ -159 +173 @@
   cMopra = strncmp(names_.instrument, "ATMOPRA", 7) == 0;
 
-  // Tidbinbilla data has some more.
-  cTid = strncmp(names_.sta, "tid", 3) == 0;
-  if (cTid) {
-    // Telescope position is stored in the source table.
+  // Non-ATNF data may not store the position in (u,v,w).
+  if (strncmp(names_.sta, "tid", 3) == 0) {
+    sprintf(cMsg, "Found Tidbinbilla data");
+    cSUpos = 1;
+  } else if (strncmp(names_.sta, "HOB", 3) == 0) {
+    sprintf(cMsg, "Found Hobart data");
+    cSUpos = 1;
+  } else if (strncmp(names_.sta, "CED", 3) == 0) {
+    sprintf(cMsg, "Found Ceduna data");
+    cSUpos = 1;
+  } else {
+    cSUpos = 0;
+  }
+
+  if (cSUpos) {
+    strcat(cMsg, ", using telescope position\n         from SU table.");
+    os << LogIO::WARN << cMsg << LogIO::POST ;
     cInterp = 0;
   }
+
+  // Mean scan rate (for timestamp repairs).
+  cNRate = 0;
+  cAvRate[0] = 0.0;
+  cAvRate[1] = 0.0;
+  cCode5 = 0;
 
 
@@ -176 +209 @@
 
   if (cNBeam <= 0) {
-    fprintf(stderr, "Couldn't determine number of beams.\n");
+    os << LogIO::SEVERE << "Couldn't determine number of beams." << LogIO::POST ;
     close();
     return 1;
@@ -189 +222 @@
   // ...beams present in the data.
   for (int iBeam = 0; iBeam < anten_.nant; iBeam++) {
-    cBeams[anten_.ant_num[iBeam] - 1] = 1;
+    // 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,
+        "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, ".");
+      }
+
+      os << LogIO::WARN << cMsg << LogIO::POST ;
+    }
+
+    if (0 < beamNo && beamNo <= cNBeam) {
+      cBeams[beamNo-1] = 1;
+    }
   }
 
@@ -237 +300 @@
   }
 
-  // Is the vis array declared by RPFITS.h large enough?
-  if (8*8193 < maxProd) {
-    // Need to allocate more memory for RPFITSIN.
-    cVis = new float[2*maxProd];
-  }
+  // 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;
@@ -278 +341 @@
   // Read the first syscal record.
   if (rpget(1, cEOS)) {
-    fprintf(stderr, "Error: Failed to read first syscal record.\n");
+    os << LogIO::SEVERE << "Failed to read first syscal record." << LogIO::POST ;
     close();
     return 1;
@@ -304 +367 @@
         double antPos[3],
         char   obsType[32],
+        char   bunit[32],
         float  &equinox,
         char   radecsys[32],
@@ -312 +376 @@
         double &bandwidth)
 {
+  const string methodName = "getHeader()" ;
+  LogIO os( LogOrigin( className, methodName, WHERE ) ) ;
+
   if (!cMBopen) {
-    fprintf(stderr, "An MBFITS file has not been opened.\n");
+    os << LogIO::SEVERE << "An MBFITS file has not been opened." << LogIO::POST ;
     return 1;
   }
@@ -333 +400 @@
     antPos[1] =  2816759.046;
     antPos[2] = -3454035.950;
+
   } else if (strncmp(names_.sta, "HOH", 3) == 0) {
     // Parkes HOH receiver.
@@ -339 +407 @@
     antPos[1] =  2816759.046;
     antPos[2] = -3454035.950;
+
   } else if (strncmp(names_.sta, "CA0", 3) == 0) {
     // An ATCA antenna, use the array centre position.
@@ -345 +414 @@
     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.
+    // Mopra.  Updated position at epoch 2007/06/24 from Chris Phillips.
     sprintf(telescope, "%-16.16s", "ATMOPRA");
-    antPos[0] = -4682768.630;
-    antPos[1] =  2802619.060;
-    antPos[2] = -3291759.900;
+    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.
@@ -357 +433 @@
     antPos[1] =  2522347.567;
     antPos[2] = -4311562.569;
+
   } else if (strncmp(names_.sta, "CED", 3) == 0) {
-    // Ceduna.
+    // Ceduna.  Updated position at epoch 2007/06/24 from Chris Phillips.
     sprintf(telescope, "%-16.16s", "CEDUNA");
-    antPos[0] = -3749943.657;
-    antPos[1] =  3909017.709;
-    antPos[2] = -3367518.309;
+    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.
@@ -379 +457 @@
   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;
@@ -386 +472 @@
   // Time at start of observation.
   sprintf(datobs, "%-10.10s", names_.datobs);
-  utc = ut;
+  utc = cUTC;
 
   // Spectral parameters.
@@ -425 +511 @@
 //--------------------------------------------------------- MBFITSreader::read
 
-// Read the next data record.
+// Read the next data record (if you're feeling lucky).
 
 int MBFITSreader::read(
-        PKSMBrecord &MBrec)
+        MBrecord &MBrec)
 {
+  const string methodName = "read()" ;
+  LogIO os( LogOrigin( className, methodName, WHERE ) ) ;
+
   int beamNo = -1;
-  int haveData, status;
-  PKSMBrecord *iMBuff = 0x0;
+  int haveData, pCode = 0, status;
+  double raRate = 0.0, decRate = 0.0, paRate = 0.0;
+  MBrecord *iMBuff = 0x0;
 
   if (!cMBopen) {
-    fprintf(stderr, "An MBFITS file has not been opened.\n");
+    os << LogIO::SEVERE << "An MBFITS file has not been opened." << LogIO::POST ;
     return 1;
   }
 
-  // Positions recorded in the input records do not coincide with the midpoint
-  // of the integration and hence the input must be buffered so that true
-  // positions may be interpolated.
+  // 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
@@ -471 +561 @@
 
       // Read the next record.
+      pCode = 0;
       if ((status = rpget(0, cEOS)) == -1) {
         // EOF.
@@ -479 +570 @@
 
 #ifdef PKSIO_DEBUG
-        printf("End-of-file detected, flushing last scan.\n");
+        os << LogIO::DEBUGGING << "\nEnd-of-file detected, flushing last cycle.\n" << LogIO::POST ;
 #endif
 
@@ -507 +598 @@
           cXpolOff = new int[cNIF];
 
-          int simulIF = 0;
           int maxChan = 0;
           int maxXpol = 0;
 
+          cSimulIF = 0;
           for (int iIF = 0; iIF < cNIF; iIF++) {
             if (cIFs[iIF]) {
@@ -536 +627 @@
 
               // Maximum number of selected IFs in any simultaneous set.
-              simulIF = max(simulIF, cIFSel[iIF]+1);
+              cSimulIF = max(cSimulIF, cIFSel[iIF]+1);
 
               // Maximum memory required for any simultaneous set.
@@ -563 +654 @@
 
           if (cNBin > 1 && cNBeamSel > 1) {
-            fprintf(stderr, "Cannot handle binning mode for multiple "
-                            "beams.\n");
+            os << LogIO::SEVERE << "Cannot handle binning mode for multiple beams.\nSelect a single beam for input." << LogIO::POST ;
             close();
             return 1;
           }
 
-          // Allocate buffer data storage.
+          // 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 PKSMBrecord[nBuff];
+          cBuffer = new MBrecord[nBuff];
 
           // Allocate memory for spectral arrays.
           for (int ibuff = 0; ibuff < nBuff; ibuff++) {
-            cBuffer[ibuff].setNIFs(simulIF);
+            cBuffer[ibuff].setNIFs(cSimulIF);
             cBuffer[ibuff].allocate(0, maxChan, maxXpol);
+
+            // Signal that this IF in this buffer has been flushed.
+            for (int iIF = 0; iIF < cSimulIF; iIF++) {
+              cBuffer[ibuff].IFno[iIF] = 0;
+            }
           }
 
@@ -584 +681 @@
           cScanNo  = 1;
           cCycleNo = 0;
-          cUTC = 0.0;
-          cStaleness = new int[cNBeamSel];
-          for (int iBeamSel = 0; iBeamSel < cNBeamSel; iBeamSel++) {
-            cStaleness[iBeamSel] = 0;
-          }
+          cPrevUTC = -1.0;
         }
 
@@ -595 +688 @@
           cScanNo++;
           cCycleNo = 0;
-          cUTC = 0.0;
-        }
-
-        // Apply beam selection.
-        beamNo = int(baseline / 256.0);
+          cPrevUTC = -1.0;
+        }
+
+        // Apply beam and IF selection before the change-of-day test to allow
+        // a single selected beam and IF to be handled in binning-mode.
+        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 (if_no > cNIF) continue;
-
-        // Apply IF selection.
-        iIFSel = cIFSel[if_no - 1];
+        if (cIFno > cNIF) continue;
+
+        // Apply IF selection; iIFSel == 0 for the first selected IF, == 1
+        // for the second, etc.
+        iIFSel = cIFSel[cIFno - 1];
         if (iIFSel < 0) continue;
 
-        sprintf(cDateObs, "%-10.10s", names_.datobs);
-
-        // Change-of-day; note that the ut variable from RPFITS.h is global
-        // and will be preserved between calls to this function.
-        if (ut < cUTC - 85800.0) {
-          ut += 86400.0;
-        }
-
-        // New integration cycle?
-        if (ut > cUTC) {
-          cCycleNo++;
-          cUTC = ut + 0.0001;
-        }
 
         if (cNBin > 1) {
           // Binning mode: correct the time.
-          ut += param_.intbase * (bin - (cNBin + 1)/2.0);
-        }
+          cUTC += param_.intbase * (cBin - (cNBin + 1)/2.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
+          char buf[256] ;
+          sprintf(buf, "Change-of-day on cUTC: %.1f -> %.1f\n", cPrevUTC, cUTC);
+          os << LogIO::DEBUGGING << buf << LogIO::POST ;
+#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;
+
+        } 
+
+        if ((cUTC+cod) < cPrevUTC - 1.0) {
+          if (cBin == 1 && iIFSel) {
+            // Multiple-IF, binning-mode data is only partially time ordered.
+#ifdef PKSIO_DEBUG
+            fprintf(stderr, "New IF in multiple-IF, binning-mode data.\n");
+#endif
+            cCycleNo -= cNBin;
+            cPrevUTC = -1.0;
+
+          } else {
+            // All other data should be fully time ordered.
+            sprintf(cMsg,
+              "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);
+            os << LogIO::WARN << cMsg << LogIO::POST ;
+            cUTC += 86400.0;
+          }
+        }
+
+        // New integration cycle?
+        if ((cUTC+cod) > cPrevUTC) {
+          cCycleNo++;
+          cPrevUTC = cUTC + 0.0001;
+        }
+
+        sprintf(cDateObs, "%-10.10s", names_.datobs);
+        cDateObs[10] = '\0';
 
         // Compute buffer number.
         iMBuff = cBuffer + iBeamSel;
-        if (cNBin > 1) iMBuff += cNBeamSel*(bin-1);
+        if (cNBin > 1) iMBuff += cNBeamSel*(cBin-1);
 
         if (cCycleNo < iMBuff->cycleNo) {
@@ -640 +777 @@
 
         // Begin flush cycle?
-        if (cEOS || (iMBuff->nIF && ut > iMBuff->utc + 0.0001)) {
+        if (cEOS || (iMBuff->nIF && (cUTC+cod) > (iMBuff->utc+0.0001))) {
           cFlushing = 1;
           cFlushBin = 0;
@@ -647 +784 @@
 
 #ifdef PKSIO_DEBUG
-        printf(" In:%4d%4d%3d%3d\n", cScanNo, cCycleNo, beamNo, if_no);
-        if (cEOS) printf("Start of new scan, flushing previous scan.\n");
+        char rel = '=';
+        double dt = utcDiff(cUTC, cW);
+        if (dt < 0.0) {
+          rel = '<';
+        } else if (dt > 0.0) {
+          rel = '>';
+        }
+
+        sprintf(buf, "\n In:%4d%4d%3d%3d  %.3f %c %.3f (%+.3fs) - "
+          "%sflushing\n", cScanNo, cCycleNo, beamNo, cIFno, cUTC, rel, cW, dt,
+          cFlushing ? "" : "not ");
+        os << LogIO::DEBUGGING << buf << LogIO::POST ;
+        if (cEOS) {
+          sprintf(buf, "Start of new scan, flushing previous scan.\n");
+          os << LogIO::DEBUGGING << buf << LogIO::POST ;
+        }
 #endif
       }
@@ -663 +814 @@
           iMBuff = cBuffer + iBeamSel + cNBeamSel*cFlushBin;
 
-          // iMBuff->nIF is set to zero (below) to signal that all IFs in
-          // an integration have been flushed.
+          // iMBuff->nIF is decremented (below) and if zero signals that all
+          // IFs in an integration have been flushed.
           if (iMBuff->nIF) {
             if (cycleNo == 0 || iMBuff->cycleNo < cycleNo) {
@@ -677 +828 @@
           break;
         }
+
+        // Start with the first IF in the next bin.
+        cFlushIF = 0;
       }
 
@@ -684 +838 @@
 
         // Find the IF to flush.
-        for (; cFlushIF < iMBuff->nIF; cFlushIF++) {
+        for (; cFlushIF < cSimulIF; cFlushIF++) {
           if (iMBuff->IFno[cFlushIF]) break;
         }
@@ -696 +850 @@
 
         // The last record read must have been the first of a new cycle.
-        beamNo = int(baseline / 256.0);
+        beamNo = int(cBaseline / 256.0);
         iBeamSel = cBeamSel[beamNo-1];
 
         // Compute buffer number.
         iMBuff = cBuffer + iBeamSel;
-        if (cNBin > 1) iMBuff += cNBeamSel*(bin-1);
+        if (cNBin > 1) iMBuff += cNBeamSel*(cBin-1);
       }
     }
 
 
-    if (cFlushing && cFlushBin == 0 && cFlushIF == 0 && cInterp) {
-      // Interpolate the beam position at the start of the flush cycle.
+    if (cInterp && cFlushing == 1) {
+      // 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
-      printf("Doing position interpolation for beam %d.\n", iMBuff->beamNo);
+      sprintf(buf, "This (%d) ra, dec, UTC: %9.4f %9.4f %10.3f %9.4f\n",
+        iMBuff->cycleNo, thisRA*R2D, thisDec*R2D, thisUTC, thisPA*R2D);
+      os << LogIO::DEBUGGING << buf << LogIO::POST ;
 #endif
 
-      double prevRA  = iMBuff->ra;
-      double prevDec = iMBuff->dec;
-      double prevUTC = cPosUTC[iBeamSel];
-
-      if (!cEOF && !cEOS) {
-        // 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
-        // 'ut', this usually being a second or two later.
-        //
-        // Note that the time recorded as the 'w' visibility coordinate
-        // cycles through 86400 back to 0 at midnight, whereas that in 'ut'
-        // continues to increase past 86400.
-
-        double thisRA  = u;
-        double thisDec = v;
-        double thisUTC = w;
-
-        if (thisUTC < prevUTC) {
-          // Must have cycled through midnight.
-          thisUTC += 86400.0;
-        }
-
-        // Guard against RA cycling through 24h in either direction.
-        if (fabs(thisRA - prevRA) > PI) {
-          if (thisRA < prevRA) {
-            thisRA += TWOPI;
+      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 {
-            thisRA -= TWOPI;
+            // (cU,cV) are probably rubbish (not yet seen in practice).
+            pCode = 2;
+            cU = thisRA;
+            cV = thisDec;
           }
-        }
-
-        // 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.
-
-        double dUTC = thisUTC - prevUTC;
-
-        // Scan rate for this beam.
-        if (dUTC > 0.0) {
-          iMBuff->raRate  = (thisRA  - prevRA)  / dUTC;
-          iMBuff->decRate = (thisDec - prevDec) / dUTC;
-
-          if (cInterp == 2) {
-            // Use the same interpolation scheme as the original pksmbfits
-            // client.  This incorrectly assumed that (thisUTC - prevUTC) 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;
-            if (utc - prevUTC > 100.0) {
-              // Must have cycled through midnight.
-              utc -= 86400.0;
+
+#ifdef PKSIO_DEBUG
+          sprintf(buf, "Next (%d) ra, dec, UTC: %9.4f %9.4f %10.3f "
+            "(0.000s)\n", cCycleNo, cU*R2D, cV*R2D, cW);
+          os << LogIO::DEBUGGING << buf << LogIO::POST ;
+#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
+          sprintf(buf, "Next (%d) ra, dec, UTC: %9.4f %9.4f %10.3f "
+            "(%+.3fs)\n", cCycleNo, nextRA*R2D, nextDec*R2D, nextUTC,
+            utcDiff(nextUTC, thisUTC));
+          os << LogIO::DEBUGGING << buf << LogIO::POST ;
+#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;
             }
 
-            double tw1 = 1.0 - (utc - prevUTC) / iMBuff->exposure;
-            double tw2 = 1.0 - (thisUTC - utc) / iMBuff->exposure;
-            double gamma = (tw2 / (tw1 + tw2)) * dUTC / (utc - prevUTC);
-
-            iMBuff->raRate  *= gamma;
-            iMBuff->decRate *= gamma;
-          }
-
-          cStaleness[iBeamSel] = 0;
-
-        } else {
-          // Issue warnings.
-          int nch = 0;
-          fprintf(stderr, "WARNING, scan %d,%n cycle %d: Position ",
-            iMBuff->scanNo, &nch, iMBuff->cycleNo);
-
-          if (dUTC < 0.0) {
-            fprintf(stderr, "timestamp went backwards!\n");
+            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 (thisRA != prevRA || thisDec != prevDec) {
-              fprintf(stderr, "changed but timestamp unchanged!\n");
+            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 {
-              fprintf(stderr, "and timestamp unchanged!\n");
+              // 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;
             }
           }
-
-          cStaleness[iBeamSel]++;
-          fprintf(stderr, "%-*s Using stale scan rate, staleness = %d "
-            "cycle%s.\n", nch, "WARNING,", cStaleness[iBeamSel],
-            (cStaleness[iBeamSel] == 1) ? "" : "s");
-
-          if (thisRA != prevRA || thisDec != prevDec) {
-            if (iMBuff->raRate == 0.0 && iMBuff->decRate == 0.0) {
-              fprintf(stderr, "%-*s But the previous rate was zero!  "
-                "Position will be inaccurate.\n", nch, "WARNING,");
+        }
+      }
+
+
+      // 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]);
+      sprintf(buf, "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);
+      os << LogIO::DEBUGGING << buf << LogIO::POST ;
+#endif
+
 
       // Compute the position of this beam for all bins.
@@ -817 +1116 @@
         cBuffer[jbuff].raRate  = iMBuff->raRate;
         cBuffer[jbuff].decRate = iMBuff->decRate;
-
-        double dutc = cBuffer[jbuff].utc - prevUTC;
-        if (dutc > 100.0) {
+        cBuffer[jbuff].paRate  = iMBuff->paRate;
+
+        double dUTC = utcDiff(cBuffer[jbuff].utc, thisUTC);
+        if (dUTC > 100.0) {
           // Must have cycled through midnight.
-          dutc -= 86400.0;
-        }
-
-        cBuffer[jbuff].ra  = prevRA  + cBuffer[jbuff].raRate  * dutc;
-        cBuffer[jbuff].dec = prevDec + cBuffer[jbuff].decRate * dutc;
-        if (cBuffer[jbuff].ra < 0.0) {
-          cBuffer[jbuff].ra += TWOPI;
-        } else if (cBuffer[jbuff].ra > TWOPI) {
-          cBuffer[jbuff].ra -= TWOPI;
-        }
-      }
+          dUTC -= 86400.0;
+        }
+
+        applyRate(cRA0, cDec0, thisRA, thisDec,
+          cBuffer[jbuff].raRate, cBuffer[jbuff].decRate, dUTC,
+          cBuffer[jbuff].ra, cBuffer[jbuff].dec);
+
+#ifdef PKSIO_DEBUG
+        sprintf(buf, "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);
+        os << LogIO::DEBUGGING << buf << LogIO::POST ;
+#endif
+      }
+
+      cFlushing = 2;
     }
 
@@ -841 +1147 @@
 
 #ifdef PKSIO_DEBUG
-      printf("Out:%4d%4d%3d%3d\n", MBrec.scanNo, MBrec.cycleNo, MBrec.beamNo,
-        MBrec.IFno[0]);
+      sprintf(buf, "Out:%4d%4d%3d%3d\n", MBrec.scanNo, MBrec.cycleNo,
+        MBrec.beamNo, MBrec.IFno[0]);
+      os << LogIO::DEBUGGING << buf << LogIO::POST ;
 #endif
 
@@ -848 +1155 @@
       iMBuff->IFno[cFlushIF] = 0;
 
-      if (cFlushIF == iMBuff->nIF - 1) {
-        // Signal that all IFs in this buffer location have been flushed.
-        iMBuff->nIF = 0;
+      iMBuff->nIF--;
+      if (iMBuff->nIF == 0) {
+        // All IFs in this buffer location have been flushed.  Stop cEOS
+        // being set when the next integration is read.
+        iMBuff->cycleNo = 0;
+
       } else {
         // Carry on flushing the other IFs.
@@ -859 +1169 @@
       if (cFlushBin == cNBin - 1) {
         if (cEOS || cEOF) {
-          // Stop cEOS being set when the next integration is read.
-          iMBuff->cycleNo = 0;
-
           // Carry on flushing other buffers.
           cFlushIF = 0;
@@ -869 +1176 @@
         cFlushing = 0;
 
-        beamNo = int(baseline / 256.0);
+        beamNo = int(cBaseline / 256.0);
         iBeamSel = cBeamSel[beamNo-1];
 
         // Compute buffer number.
         iMBuff = cBuffer + iBeamSel;
-        if (cNBin > 1) iMBuff += cNBeamSel*(bin-1);
+        if (cNBin > 1) iMBuff += cNBeamSel*(cBin-1);
       }
     }
@@ -880 +1187 @@
     if (!cFlushing) {
       // Buffer this MBrec.
-      if (cCycleNo == 1 && iMBuff->IFno[0]) {
+      if ((cScanNo > iMBuff->scanNo) && iMBuff->IFno[0]) {
         // Sanity check on the number of IFs in the new scan.
         if (if_.n_if != cNIF) {
-          fprintf(stderr, "WARNING, scan %d has %d IFs instead of %d, "
-            "continuing.\n", cScanNo, if_.n_if, cNIF);
-        }
-      }
-
-      iMBuff->scanNo  = cScanNo;
-      iMBuff->cycleNo = cCycleNo;
-
-      // Times.
-      strncpy(iMBuff->datobs, cDateObs, 10);
-      iMBuff->utc = ut;
-      iMBuff->exposure = param_.intbase;
-
-      // Source identification.
-      sprintf(iMBuff->srcName, "%-16.16s",
-              names_.su_name + (sourceno-1)*16);
-      iMBuff->srcRA  = doubles_.su_ra[sourceno-1];
-      iMBuff->srcDec = doubles_.su_dec[sourceno-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.40575e6) < 100.0) {
-          iMBuff->restFreq = 1420.40575e6;
-        }
-      }
-
-      // 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 (cTid) {
-        // Tidbinbilla data.
-        iMBuff->ra  = doubles_.su_ra[sourceno-1];
-        iMBuff->dec = doubles_.su_dec[sourceno-1];
-      } else {
-        iMBuff->ra  = u;
-        iMBuff->dec = v;
-      }
-      cPosUTC[iBeamSel] = w;
+          sprintf(cMsg, "Scan %d has %d IFs instead of %d, "
+            "continuing.", cScanNo, if_.n_if, cNIF);
+          os << LogIO::WARN << cMsg << LogIO::POST ;
+        }
+      }
+
+      // 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
+      sprintf(buf, "Buf:%4d%4d%3d%3d\n", cScanNo, cCycleNo, beamNo, cIFno);
+      os << LogIO::DEBUGGING << buf << LogIO::POST ;
+#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 = if_no - 1;
+      int iIF = cIFno - 1;
       int startChan = cStartChan[iIF];
       int endChan   = cEndChan[iIF];
@@ -945 +1278 @@
 
       iIFSel = cIFSel[iIF];
-      iMBuff->nIF++;
-      iMBuff->IFno[iIFSel]  = if_no;
+      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, "Integration cycle %d:%d, beam %2d, \n"
+                      "IF %d was duplicated.", cScanNo, cCycleNo-1,
+                      beamNo, cIFno);
+        os << LogIO::WARN << cMsg << LogIO::POST ;
+      }
       iMBuff->nChan[iIFSel] = nChan;
       iMBuff->nPol[iIFSel]  = cNPol[iIF];
@@ -1005 +1347 @@
 
       // The baseline flag may be set independently.
-      if (rpflag == 0) rpflag = flag;
+      if (rpflag == 0) rpflag = cFlag;
 
       // Copy and scale data.
@@ -1046 +1388 @@
 
 
-      // Parallactic angle.
-      iMBuff->parAngle = sc_.sc_cal[scq*iBeam + 11];
-
       // Calibration factor applied to the data by the correlator.
       if (scq > 14) {
@@ -1059 +1398 @@
       }
 
-      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]));
+      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;
-        }
+          // 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;
+          // 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;
+          }
         }
       }
@@ -1128 +1473 @@
 int MBFITSreader::rpget(int syscalonly, int &EOS)
 {
+  const string methodName = "rpget()" ;
+  LogIO os( LogOrigin( className, methodName, WHERE ) ) ;
+
   EOS = 0;
 
@@ -1135 +1483 @@
   int numErr = 0;
 
-  jstat = 0;
+  int jstat = 0;
   while (numErr < 10) {
     int lastjstat = jstat;
-    rpfitsin_(&jstat, cVis, weight, &baseline, &ut, &u, &v, &w, &flag, &bin,
-              &if_no, &sourceno);
-
-    switch(jstat) {
+
+    switch(rpfitsin(jstat)) {
     case -1:
       // Read failed; retry.
       numErr++;
-      fprintf(stderr, "RPFITS read failed - retrying.\n");
+      os << LogIO::WARN << "RPFITS read failed - retrying." << LogIO::POST ;
       jstat = 0;
       break;
@@ -1152 +1498 @@
       // Successful read.
       if (lastjstat == 0) {
-        if (baseline == -1) {
+        if (cBaseline == -1) {
           // Syscal data.
           if (syscalonly) {
@@ -1201 +1547 @@
     default:
       // Shouldn't reach here.
-      fprintf(stderr, "Unrecognized RPFITSIN return code: %d (retrying)\n",
-              jstat);
+      sprintf(cMsg, "Unrecognized RPFITSIN return code: %d "
+                    "(retrying).", jstat);
+      os << LogIO::WARN << cMsg << LogIO::POST ;
       jstat = 0;
       break;
@@ -1208 +1555 @@
   }
 
-  fprintf(stderr, "RPFITS read failed too many times.\n");
+  os << LogIO::SEVERE << "RPFITS read failed too many times." << LogIO::POST ;
   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;
 }
 
@@ -1219 +1730 @@
 {
   if (cMBopen) {
-    jstat = 1;
-    rpfitsin_(&jstat, cVis, weight, &baseline, &ut, &u, &v, &w, &flag, &bin,
-              &if_no, &sourceno);
+    int jstat = 1;
+    rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
+              &cBin, &cIFno, &cSrcNo);
 
     if (cBeams)     delete [] cBeams;
@@ -1232 +1743 @@
     if (cRefChan)   delete [] cRefChan;
 
-    if (cVis)       delete [] cVis;
+    if (cVis) delete [] cVis;
+    if (cWgt) delete [] cWgt;
 
     if (cBeamSel)   delete [] cBeamSel;
@@ -1244 +1756 @@
   }
 }
+
+//-------------------------------------------------------------------- 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);
+}