Changeset 1452 for trunk/external/atnf/PKSIO/MBFITSreader.cc

Timestamp:

11/19/08 20:41:16 (16 years ago)

Author:

Malte Marquarding

Message:

update from livedata CVS

File:

• trunk/external/atnf/PKSIO/MBFITSreader.cc (modified) (43 diffs)

trunk/external/atnf/PKSIO/MBFITSreader.cc

@@ -27 +27 @@
 //#                        AUSTRALIA
 //#
-//# $Id: MBFITSreader.cc,v 19.38 2008-06-26 02:24:22 cal103 Exp $
+//# $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
@@ -35 +35 @@
 //#---------------------------------------------------------------------------
 
+#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/math.h>
@@ -47 +46 @@
 #include <unistd.h>
 
+#include <RPFITS.h>
+
 using namespace std;
 
@@ -52 +53 @@
 const double PI = 3.141592653589793238462643;
 const double TWOPI = 2.0 * PI;
+const double HALFPI = PI / 2.0;
 const double R2D = 180.0 / PI;
 
@@ -93 +95 @@
 
   cMBopen = 0;
+
+  // Tell RPFITSIN not to report errors directly.
+  iostat_.errlun = -1;
+
+  // By default, messages are written to stderr.
+  initMsg();
 }
 
@@ -121 +129 @@
         int  &extraSysCal)
 {
+  // Clear the message stack.
+  clearMsg();
+
   if (cMBopen) {
     close();
@@ -129 +140 @@
   // Open the RPFITS file.
   int jstat = -3;
-  rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
-            &cBin, &cIFno, &cSrcNo);
-
-  if (jstat) {
-    fprintf(stderr, "ERROR, failed to open MBFITS file: %s\n", rpname);
+  if (rpfitsin(jstat)) {
+    sprintf(cMsg, "ERROR: failed to open MBFITS file\n       %s", rpname);
+    logMsg(cMsg);
     return 1;
   }
@@ -149 +158 @@
   // Read the first header.
   jstat = -1;
-  rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
-            &cBin, &cIFno, &cSrcNo);
-
-  if (jstat) {
-    fprintf(stderr, "ERROR, failed to read MBFITS header: %s\n", rpname);
+  if (rpfitsin(jstat)) {
+    sprintf(cMsg, "ERROR: failed to read MBFITS header in file\n"
+                  "       %s", rpname);
+    logMsg(cMsg);
     close();
     return 1;
@@ -163 +171 @@
   // Non-ATNF data may not store the position in (u,v,w).
   if (strncmp(names_.sta, "tid", 3) == 0) {
-    fprintf(stderr, "WARNING, found Tidbinbilla data");
+    sprintf(cMsg, "WARNING: Found Tidbinbilla data");
     cSUpos = 1;
   } else if (strncmp(names_.sta, "HOB", 3) == 0) {
-    fprintf(stderr, "WARNING, found Hobart data");
+    sprintf(cMsg, "WARNING: Found Hobart data");
     cSUpos = 1;
   } else if (strncmp(names_.sta, "CED", 3) == 0) {
-    fprintf(stderr, "WARNING, found Ceduna data");
+    sprintf(cMsg, "WARNING: Found Ceduna data");
     cSUpos = 1;
   } else {
@@ -176 +184 @@
 
   if (cSUpos) {
-    fprintf(stderr, ", using telescope position from SU table.\n");
+    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;
 
 
@@ -190 +205 @@
 
   if (cNBeam <= 0) {
-    fprintf(stderr, "ERROR, couldn't determine number of beams.\n");
+    logMsg("ERROR, couldn't determine number of beams.");
     close();
     return 1;
@@ -203 +218 @@
   // ...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,
+        "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;
+    }
   }
 
@@ -292 +337 @@
   // Read the first syscal record.
   if (rpget(1, cEOS)) {
-    fprintf(stderr, "ERROR, failed to read first syscal record.\n");
+    logMsg("ERROR, failed to read first syscal record.");
     close();
     return 1;
@@ -328 +373 @@
 {
   if (!cMBopen) {
-    fprintf(stderr, "ERROR, an MBFITS file has not been opened.\n");
+    logMsg("ERROR, an MBFITS file has not been opened.");
     return 1;
   }
@@ -348 +393 @@
     antPos[1] =  2816759.046;
     antPos[2] = -3454035.950;
+
   } else if (strncmp(names_.sta, "HOH", 3) == 0) {
     // Parkes HOH receiver.
@@ -354 +400 @@
     antPos[1] =  2816759.046;
     antPos[2] = -3454035.950;
+
   } else if (strncmp(names_.sta, "CA0", 3) == 0) {
     // An ATCA antenna, use the array centre position.
@@ -360 +407 @@
     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.
@@ -372 +426 @@
     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.
@@ -448 +504 @@
 //--------------------------------------------------------- MBFITSreader::read
 
-// Read the next data record.
+// Read the next data record (if you're feeling lucky).
 
 int MBFITSreader::read(
-        PKSMBrecord &MBrec)
+        MBrecord &MBrec)
 {
   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, "ERROR, an MBFITS file has not been opened.\n");
+    logMsg("ERROR, an MBFITS file has not been opened.");
     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
@@ -494 +551 @@
 
       // Read the next record.
+      pCode = 0;
       if ((status = rpget(0, cEOS)) == -1) {
         // EOF.
@@ -502 +560 @@
 
 #ifdef PKSIO_DEBUG
-        printf("\nEnd-of-file detected, flushing last scan.\n");
+        fprintf(stderr, "\nEnd-of-file detected, flushing last cycle.\n");
 #endif
 
@@ -586 +644 @@
 
           if (cNBin > 1 && cNBeamSel > 1) {
-            fprintf(stderr, "ERROR, cannot handle binning mode for multiple "
-                            "beams.\n");
+            logMsg("ERROR, cannot handle binning mode for multiple beams.");
             close();
             return 1;
           }
 
-          // Allocate buffer data storage; the PKSMBrecord constructor zeroes
+          // 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.
@@ -602 +659 @@
             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;
+            }
           }
 
@@ -609 +671 @@
           cScanNo  = 1;
           cCycleNo = 0;
-          cPrevUTC = 0.0;
+          cPrevUTC = -1.0;
         }
 
@@ -616 +678 @@
           cScanNo++;
           cCycleNo = 0;
-          cPrevUTC = 0.0;
+          cPrevUTC = -1.0;
         }
 
         // Check for change-of-day.
-        if (cUTC < cPrevUTC - 85800.0) {
+        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;
         }
@@ -630 +713 @@
 
         // New integration cycle?
-        if (cUTC > cPrevUTC) {
+        if ((cUTC+cod) > cPrevUTC) {
           cCycleNo++;
           cPrevUTC = cUTC + 0.0001;
@@ -637 +720 @@
         // 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;
@@ -648 +736 @@
 
         sprintf(cDateObs, "%-10.10s", names_.datobs);
+        cDateObs[10] = '\0';
 
         // Compute buffer number.
@@ -660 +749 @@
 
         // Begin flush cycle?
-        if (cEOS || (iMBuff->nIF && cUTC > iMBuff->utc + 0.0001)) {
+        if (cEOS || (iMBuff->nIF && (cUTC+cod) > (iMBuff->utc+0.0001))) {
           cFlushing = 1;
           cFlushBin = 0;
@@ -667 +756 @@
 
 #ifdef PKSIO_DEBUG
-        printf("\n In:%4d%4d%3d%3d\n", cScanNo, cCycleNo, beamNo, cIFno);
-        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 = '>';
+        }
+
+        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
       }
@@ -733 +834 @@
       // 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.
+      // 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 will be
-      // interpolated to the midpoint of that integration using the position
-      // recorded in the "next" integration which is currently sitting in the
-      // RPFITS commons.  The interpolation method used here is based on the
-      // scan rate.  At the end of a scan, or if the next position has not
-      // been updated, the most recent determination of the scan rate will be
-      // used for extrapolation.
+      // 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.
       //
-      // The rate "age" is the offset from "this" integration (in iMBuff) of
-      // the earliest integration in the pair used to compute the rate.  A
-      // rate "age" of 0 thus refers to the normal situation where the rate
-      // is determined from "this" integration and the "next" one.  An age
-      // of 1 cycle means that it is determined from "this" integration and
-      // the one preceding it, which should be equally reliable.  An age
-      // of 2 cycles means that the rate is determined from the previous
-      // integration and the one before that, so the extrapolation spans one
-      // integration cycle.  Thus it has a "staleness" of 1.
-
+      // 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) {
-        iMBuff->rateAge++;
-        iMBuff->rateson = 0;
-        
+        // Use rates from the last cycle.
+        raRate  = iMBuff->raRate;
+        decRate = iMBuff->decRate;
+        paRate  = iMBuff->paRate;
+
       } else {
-        // Note that the time recorded as the 'w' visibility coordinate
-        // cycles through 86400 back to 0 at midnight, whereas that in 'cUTC'
-        // continues to increase past 86400.
-
-        double nextRA  = cU;
-        double nextDec = cV;
-        double nextUTC = cW;
-
-        if (nextUTC < thisUTC) {
-          // Must have cycled through midnight.
-          nextUTC += 86400.0;
-        }
-
-        // Guard against RA cycling through 24h in either direction.
-        if (fabs(nextRA - thisRA) > PI) {
-          if (nextRA < thisRA) {
-            nextRA += TWOPI;
+        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 {
-            nextRA -= TWOPI;
+            // (cU,cV) are probably rubbish (not yet seen in practice).
+            pCode = 2;
+            cU = thisRA;
+            cV = thisDec;
           }
-        }
 
 #ifdef PKSIO_DEBUG
-        printf("Previous ra, dec, UTC:  %8.4f  %8.4f  %7.1f\n", thisRA*R2D,
-          thisDec*R2D, thisUTC);
-        printf("Current  ra, dec, UTC:  %8.4f  %8.4f  %7.1f\n", nextRA*R2D,
-          nextDec*R2D, nextUTC);
+          fprintf(stderr, "Next (%d) ra, dec, UTC: %9.4f %9.4f %10.3f "
+            "(0.000s)\n", cCycleNo, cU*R2D, cV*R2D, cW);
 #endif
 
-        // 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 = nextUTC - thisUTC;
-
-        // Scan rate for this beam.
-        if (dUTC > 0.0) {
-          iMBuff->raRate  = (nextRA  - thisRA)  / dUTC;
-          iMBuff->decRate = (nextDec - thisDec) / dUTC;
-          iMBuff->rateAge = 0;
-          iMBuff->rateson = 0;
-
-          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;
-            if (utc - thisUTC > 100.0) {
-              // Must have cycled through midnight.
-              utc -= 86400.0;
+        } 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;
             }
 
-            double tw1 = 1.0 - (utc - thisUTC) / iMBuff->exposure;
-            double tw2 = 1.0 - (nextUTC - utc) / iMBuff->exposure;
-            double gamma = (tw2 / (tw1 + tw2)) * dUTC / (utc - thisUTC);
-
-            iMBuff->raRate  *= gamma;
-            iMBuff->decRate *= gamma;
-          }
-
-        } else {
-          iMBuff->rateAge++;
-
-          // Staleness codes.
-          if (dUTC < 0.0) {
-            iMBuff->rateson = 3;
+            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 (nextRA != thisRA || nextDec != thisDec) {
-              iMBuff->rateson = 2;
+            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 {
-              iMBuff->rateson = 1;
+              // 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;
             }
           }
@@ -844 +1004 @@
       }
 
+
+      // 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
-      printf("Doing position interpolation for beam %d.\n", iMBuff->beamNo);
-      printf("RA and Dec rates and age: %7.4f  %7.4f  %d\n",
-        iMBuff->raRate*R2D, iMBuff->decRate*R2D, iMBuff->rateAge);
+      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.
@@ -856 +1079 @@
         cBuffer[jbuff].raRate  = iMBuff->raRate;
         cBuffer[jbuff].decRate = iMBuff->decRate;
-
-        double dutc = cBuffer[jbuff].utc - thisUTC;
-        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  = thisRA  + cBuffer[jbuff].raRate  * dutc;
-        cBuffer[jbuff].dec = thisDec + 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
+        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
       }
     }
@@ -880 +1107 @@
 
 #ifdef PKSIO_DEBUG
-      printf("Out:%4d%4d%3d%3d\n", MBrec.scanNo, MBrec.cycleNo, MBrec.beamNo,
-        MBrec.IFno[0]);
+      fprintf(stderr, "Out:%4d%4d%3d%3d\n", MBrec.scanNo, MBrec.cycleNo,
+        MBrec.beamNo, MBrec.IFno[0]);
 #endif
 
@@ -923 +1150 @@
         // 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);
+          sprintf(cMsg, "WARNING: Scan %d has %d IFs instead of %d, "
+            "continuing.", cScanNo, if_.n_if, cNIF);
+          logMsg(cMsg);
         }
       }
@@ -935 +1163 @@
       }
 
-      iMBuff->scanNo  = cScanNo;
-      iMBuff->cycleNo = cCycleNo;
-
-      // Times.
-      strncpy(iMBuff->datobs, cDateObs, 10);
-      iMBuff->utc = cUTC;
-      iMBuff->exposure = param_.intbase;
-
-      // Source identification.
-      sprintf(iMBuff->srcName, "%-16.16s",
-              names_.su_name + (cSrcNo-1)*16);
-      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.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 (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;
+#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.
@@ -992 +1237 @@
 
       iIFSel = cIFSel[iIF];
-      iMBuff->nIF++;
-      iMBuff->IFno[iIFSel]  = cIFno;
+      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];
@@ -1093 +1347 @@
 
 
-      // Parallactic angle.
-      iMBuff->parAngle = sc_.sc_cal[scq*iBeam + 11];
-
       // Calibration factor applied to the data by the correlator.
       if (scq > 14) {
@@ -1106 +1357 @@
       }
 
-      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;
+          }
         }
       }
@@ -1185 +1442 @@
   while (numErr < 10) {
     int lastjstat = jstat;
-    rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
-              &cBin, &cIFno, &cSrcNo);
-
-    switch(jstat) {
+
+    switch(rpfitsin(jstat)) {
     case -1:
       // Read failed; retry.
       numErr++;
-      fprintf(stderr, "RPFITS read failed - retrying.\n");
+      logMsg("WARNING: RPFITS read failed - retrying.");
       jstat = 0;
       break;
@@ -1248 +1503 @@
     default:
       // Shouldn't reach here.
-      fprintf(stderr, "Unrecognized RPFITSIN return code: %d (retrying)\n",
-              jstat);
+      sprintf(cMsg, "WARNING: Unrecognized RPFITSIN return code: %d "
+                    "(retrying).", jstat);
+      logMsg(cMsg);
       jstat = 0;
       break;
@@ -1255 +1511 @@
   }
 
-  fprintf(stderr, "ERROR, RPFITS read failed too many times.\n");
+  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;
 }
 
@@ -1292 +1711 @@
   }
 }
+
+//-------------------------------------------------------------------- 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);
+}