Changeset 1452 for trunk/external/atnf/PKSIO

Timestamp:

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

Author:

Malte Marquarding

Message:

update from livedata CVS

Location:

trunk/external/atnf/PKSIO

Files:

• FITSreader.cc (modified) (3 diffs)
• FITSreader.h (modified) (7 diffs)
• MBFITSreader.cc (modified) (43 diffs)
• MBFITSreader.h (modified) (6 diffs)
• MBrecord.cc (added)
• MBrecord.h (added)
• PKSFITSreader.cc (modified) (14 diffs)
• PKSFITSreader.h (modified) (7 diffs)
• PKSMBrecord.cc (deleted)
• PKSMBrecord.h (deleted)
• PKSMS2reader.cc (modified) (15 diffs)
• PKSMS2reader.h (modified) (5 diffs)
• PKSMS2writer.cc (modified) (17 diffs)
• PKSMS2writer.h (modified) (5 diffs)
• PKSSDwriter.cc (modified) (16 diffs)
• PKSSDwriter.h (modified) (5 diffs)
• PKSmsg.cc (added)
• PKSmsg.h (added)
• PKSreader.cc (modified) (7 diffs)
• PKSreader.h (modified) (9 diffs)
• PKSrecord.h (added)
• PKSwriter.h (modified) (3 diffs)
• SDFITSreader.cc (modified) (47 diffs)
• SDFITSreader.h (modified) (6 diffs)
• SDFITSwriter.cc (modified) (15 diffs)
• SDFITSwriter.h (modified) (7 diffs)
• makefile (modified) (1 diff)

trunk/external/atnf/PKSIO/FITSreader.cc

@@ -27 +27 @@
 //#                        AUSTRALIA
 //#
-//# $Id: FITSreader.cc,v 19.2 2006/05/19 02:14:50 mcalabre Exp $
+//# $Id: FITSreader.cc,v 19.3 2008-09-12 02:37:37 cal103 Exp $
 //#---------------------------------------------------------------------------
 //# The FITSreader class is an abstract base class for the Parkes Multibeam
@@ -53 +53 @@
         const int getSpectra,
         const int getXPol,
-        const int getFeedPos)
+        const int coordSys)
 {
   int maxNChan = 0;
@@ -83 +83 @@
   cGetSpectra = getSpectra && cHaveSpectra;
   cGetXPol    = getXPol    && cGetXPol;
-  cGetFeedPos = getFeedPos;
+  cCoordSys   = coordSys;
 
   return maxNChan;

trunk/external/atnf/PKSIO/FITSreader.h

@@ -27 +27 @@
 //#                        AUSTRALIA
 //#
-//# $Id: FITSreader.h,v 19.6 2007/11/12 03:37:56 cal103 Exp $
+//# $Id: FITSreader.h,v 19.9 2008-11-17 06:28:04 cal103 Exp $
 //#---------------------------------------------------------------------------
 //# The FITSreader class is an abstract base class for the Parkes Multibeam
@@ -38 +38 @@
 #define ATNF_FITSREADER_H
 
-#include <atnf/PKSIO/PKSMBrecord.h>
+#include <atnf/PKSIO/MBrecord.h>
+#include <atnf/PKSIO/PKSmsg.h>
+
+using namespace std;
 
 // <summary>
@@ -44 +47 @@
 // </summary>
 
-class FITSreader
+class FITSreader : public PKSmsg
 {
   public:
@@ -91 +94 @@
     // Set data selection criteria.  Channel numbering is 1-relative, zero or
     // negative channel numbers are taken to be offsets from the last channel.
+    // Coordinate systems are
+    //   0: equatorial (RA,Dec),
+    //   1: vertical (Az,El),
+    //   2: feed-plane.
     int select(
         const int startChan[],
@@ -96 +103 @@
         const int refChan[],
         const int getSpectra = 1,
-        const int getXPol = 0,
-        const int getFeedPos = 0);
+        const int getXPol  = 0,
+        const int coordSys = 0);
 
     // Find the range in time and position of the data selected.
@@ -109 +116 @@
     // Read the next data record.
     virtual int read(
-        PKSMBrecord &record) = 0;
+        MBrecord &record) = 0;
 
     // Close the RPFITS file.
@@ -115 +122 @@
 
   protected:
-    int    *cBeams, *cEndChan, cGetFeedPos, cGetSpectra, cGetXPol, cHaveBase,
-           cHaveSpectra, *cHaveXPol, *cIFs, cNBeam, *cNChan, cNIF, *cNPol,
-           *cRefChan, *cStartChan;
+    int  *cBeams, *cEndChan, cCoordSys, cGetSpectra, cGetXPol, cHaveBase,
+         cHaveSpectra, *cHaveXPol, *cIFs, cNBeam, *cNChan, cNIF, *cNPol,
+         *cRefChan, *cStartChan;
+
+    // For use in constructing messages.
+    char cMsg[256];
 };
 

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);
+}

trunk/external/atnf/PKSIO/MBFITSreader.h

@@ -27 +27 @@
 //#                        AUSTRALIA
 //#
-//# $Id: MBFITSreader.h,v 19.15 2008-06-26 02:14:36 cal103 Exp $
+//# $Id: MBFITSreader.h,v 19.21 2008-11-17 06:33:10 cal103 Exp $
 //#---------------------------------------------------------------------------
 //# The MBFITSreader class reads single dish RPFITS files (such as Parkes
@@ -39 +39 @@
 
 #include <atnf/PKSIO/FITSreader.h>
-#include <atnf/PKSIO/PKSMBrecord.h>
+#include <atnf/PKSIO/MBrecord.h>
+
+using namespace std;
 
 // <summary>
@@ -102 +104 @@
 
     // Read the next data record.
-    virtual int read(PKSMBrecord &record);
+    virtual int read(MBrecord &record);
 
     // Close the RPFITS file.
@@ -112 +114 @@
     float cUTC, cU, cV, cW, *cVis, *cWgt;
 
-    char   cDateObs[10];
+    char   cDateObs[12];
     int    *cBeamSel, *cChanOff, cFirst, *cIFSel, cInterp, cIntTime, cMBopen,
            cMopra, cNBeamSel, cNBin, cRetry, cSUpos, *cXpolOff;
@@ -119 +121 @@
     int    cEOF, cEOS, cFlushBin, cFlushIF, cFlushing;
     double *cPosUTC;
-    PKSMBrecord *cBuffer;
+    MBrecord *cBuffer;
 
+    // Scan and cycle number bookkeeping.
     int    cCycleNo, cScanNo;
     double cPrevUTC;
@@ -126 +129 @@
     // Read the next data record from the RPFITS file.
     int rpget(int syscalonly, int &EOS);
+    int rpfitsin(int &jstat);
+
+    // Check and, if necessary, repair a position timestamp.
+    int    cCode5, cNRate;
+    double cAvRate[2];
+    int fixw(const char *datobs, int cycleNo, int beamNo, double avRate[2],
+             double thisRA, double thisDec, double thisUTC,
+             double nextRA, double nextDec, float &nextUTC);
+
+    // Subtract two UTCs (s).
+    double utcDiff(double utc1, double utc2);
+
+    // Compute and apply the scan rate corrected for grid convergence.
+    double cRA0, cDec0;
+    void  scanRate(double ra0, double dec0,
+                   double ra1, double dec1,
+                   double ra2, double dec2, double dt,
+                   double &raRate, double &decRate);
+    void applyRate(double ra0, double dec0,
+                   double ra1, double dec1,
+                   double raRate, double decRate, double dt,
+                   double &ra2, double &dec2);
+    void eulerx(double lng0, double lat0, double phi0, double theta,
+                double phi, double &lng1, double &lat1);
 };
 

trunk/external/atnf/PKSIO/PKSFITSreader.cc

@@ -25 +25 @@
 //#                        Charlottesville, VA 22903-2475 USA
 //#
-//# $Id: PKSFITSreader.cc,v 19.14 2008-06-26 02:07:21 cal103 Exp $
+//# $Id: PKSFITSreader.cc,v 19.21 2008-11-17 06:54:25 cal103 Exp $
 //#---------------------------------------------------------------------------
 //# Original: 2000/08/02, Mark Calabretta, ATNF
 //#---------------------------------------------------------------------------
 
+#include <atnf/PKSIO/PKSmsg.h>
 #include <atnf/PKSIO/MBFITSreader.h>
 #include <atnf/PKSIO/SDFITSreader.h>
 #include <atnf/PKSIO/PKSFITSreader.h>
-#include <atnf/PKSIO/PKSMBrecord.h>
-
+#include <atnf/PKSIO/PKSrecord.h>
+
+#include <casa/stdio.h>
 #include <casa/Arrays/Array.h>
 #include <casa/BasicMath/Math.h>
 #include <casa/Quanta/MVTime.h>
-
-#include <casa/stdio.h>
 
 //----------------------------------------------- PKSFITSreader::PKSFITSreader
@@ -50 +50 @@
         const Bool   interpolate)
 {
-  cFITSMBrec.setNIFs(1);
+  cMBrec.setNIFs(1);
 
   if (fitsType == "SDFITS") {
@@ -57 +57 @@
     cReader = new MBFITSreader(retry, interpolate ? 1 : 0);
   }
+
+  // By default, messages are written to stderr.
+  initMsg();
 }
 
@@ -67 +70 @@
   close();
   delete cReader;
+}
+
+//------------------------------------------------------ PKSFITSreader::setMsg
+
+// Set message disposition.  If fd is non-zero messages will be written
+// to that file descriptor, else stored for retrieval by getMsg().
+
+Int PKSFITSreader::setMsg(FILE *fd)
+{
+  PKSmsg::setMsg(fd);
+  cReader->setMsg(fd);
+
+  return 0;
 }
 
@@ -83 +99 @@
         Bool   &haveSpectra)
 {
+  clearMsg();
+
   int    extraSysCal, haveBase_, *haveXPol_, haveSpectra_, nBeam, *nChan_,
          nIF, *nPol_, status;
-  if ((status = cReader->open((char *)fitsName.chars(), nBeam, cBeams, nIF,
-                              cIFs, nChan_, nPol_, haveXPol_, haveBase_,
-                              haveSpectra_, extraSysCal))) {
+  status = cReader->open((char *)fitsName.chars(), nBeam, cBeams, nIF, cIFs,
+                         nChan_, nPol_, haveXPol_, haveBase_, haveSpectra_,
+                         extraSysCal);
+  logMsg(cReader->getMsg());
+  cReader->clearMsg();
+  if (status) {
     return status;
   }
@@ -146 +167 @@
   char   bunit_[32], datobs[32], dopplerFrame_[32], observer_[32],
          obsType_[32], project_[32], radecsys[32], telescope[32];
+  int    status;
   float  equinox_;
   double antPos[3], utc;
 
-  if (cReader->getHeader(observer_, project_, telescope, antPos, obsType_,
-                         bunit_, equinox_, radecsys, dopplerFrame_,
-                         datobs, utc, refFreq, bandwidth)) {
+  status = cReader->getHeader(observer_, project_, telescope, antPos,
+                              obsType_, bunit_, equinox_, radecsys,
+                              dopplerFrame_, datobs, utc, refFreq, bandwidth);
+  logMsg(cReader->getMsg());
+  cReader->clearMsg();
+  if (status) {
     return 1;
   }
@@ -187 +212 @@
   double *startfreq, *endfreq;
 
-  Int status;
-  if (!(status = cReader->getFreqInfo(nIF, startfreq, endfreq))) {
+  Int status = cReader->getFreqInfo(nIF, startfreq, endfreq);
+
+  logMsg(cReader->getMsg());
+  cReader->clearMsg();
+  if (!status) {
     startFreq.takeStorage(IPosition(1,nIF), startfreq, TAKE_OVER);
     endFreq.takeStorage(IPosition(1,nIF), endfreq, TAKE_OVER);
@@ -208 +236 @@
         const Bool getSpectra,
         const Bool getXPol,
-        const Bool getFeedPos)
+        const Int  coordSys)
 {
   // Apply beam selection.
@@ -275 +303 @@
   cGetSpectra = getSpectra;
   cGetXPol    = getXPol;
-  cGetFeedPos = getFeedPos;
+  cCoordSys   = coordSys;
 
   uInt maxNChan = cReader->select(start, end, ref, cGetSpectra, cGetXPol,
-                                  cGetFeedPos);
+                                  cCoordSys);
+  logMsg(cReader->getMsg());
+  cReader->clearMsg();
 
   delete [] end;
@@ -301 +331 @@
   double* posns;
 
-  Int status;
-  if (!(status = cReader->findRange(nRow, nSel, dateSpan, utcSpan, posns))) {
+  Int status = cReader->findRange(nRow, nSel, dateSpan, utcSpan, posns);
+  logMsg(cReader->getMsg());
+  cReader->clearMsg();
+
+  if (!status) {
     timeSpan.resize(2);
 
@@ -321 +354 @@
 // Read the next data record.
 
-Int PKSFITSreader::read(MBrecord &MBrec)
-{
-  Int status;
-
-  if ((status = cReader->read(cFITSMBrec))) {
+Int PKSFITSreader::read(PKSrecord &pksrec)
+{
+  Int status = cReader->read(cMBrec);
+  logMsg(cReader->getMsg());
+  cReader->clearMsg();
+
+  if (status) {
     if (status != -1) {
       status = 1;
@@ -334 +369 @@
 
 
-  uInt nChan = cFITSMBrec.nChan[0];
-  uInt nPol  = cFITSMBrec.nPol[0];
-
-  MBrec.scanNo  = cFITSMBrec.scanNo;
-  MBrec.cycleNo = cFITSMBrec.cycleNo;
+  uInt nChan = cMBrec.nChan[0];
+  uInt nPol  = cMBrec.nPol[0];
+
+  pksrec.scanNo  = cMBrec.scanNo;
+  pksrec.cycleNo = cMBrec.cycleNo;
 
   // Extract MJD.
   Int day, month, year;
-  sscanf(cFITSMBrec.datobs, "%4d-%2d-%2d", &year, &month, &day);
-  MBrec.mjd = MVTime(year, month, Double(day)).day() + cFITSMBrec.utc/86400.0;
-
-  MBrec.interval  = cFITSMBrec.exposure;
-
-  MBrec.fieldName = trim(cFITSMBrec.srcName);
-  MBrec.srcName   = MBrec.fieldName;
-
-  MBrec.srcDir.resize(2);
-  MBrec.srcDir(0) = cFITSMBrec.srcRA;
-  MBrec.srcDir(1) = cFITSMBrec.srcDec;
-
-  MBrec.srcPM.resize(2);
-  MBrec.srcPM(0)  = 0.0;
-  MBrec.srcPM(1)  = 0.0;
-  MBrec.srcVel    = 0.0;
-  MBrec.obsType   = trim(cFITSMBrec.obsType);
-
-  MBrec.IFno = cFITSMBrec.IFno[0];
-  Double chanWidth = fabs(cFITSMBrec.fqDelt[0]);
-  MBrec.refFreq   = cFITSMBrec.fqRefVal[0];
-  MBrec.bandwidth = chanWidth * nChan;
-  MBrec.freqInc   = cFITSMBrec.fqDelt[0];
-  MBrec.restFreq  = cFITSMBrec.restFreq;
-
-  MBrec.tcal.resize(nPol);
+  sscanf(cMBrec.datobs, "%4d-%2d-%2d", &year, &month, &day);
+  pksrec.mjd = MVTime(year, month, Double(day)).day() + cMBrec.utc/86400.0;
+
+  pksrec.interval  = cMBrec.exposure;
+
+  pksrec.fieldName = trim(cMBrec.srcName);
+  pksrec.srcName   = pksrec.fieldName;
+
+  pksrec.srcDir.resize(2);
+  pksrec.srcDir(0) = cMBrec.srcRA;
+  pksrec.srcDir(1) = cMBrec.srcDec;
+
+  pksrec.srcPM.resize(2);
+  pksrec.srcPM(0)  = 0.0;
+  pksrec.srcPM(1)  = 0.0;
+  pksrec.srcVel    = 0.0;
+  pksrec.obsType   = trim(cMBrec.obsType);
+
+  pksrec.IFno = cMBrec.IFno[0];
+  Double chanWidth = fabs(cMBrec.fqDelt[0]);
+  pksrec.refFreq   = cMBrec.fqRefVal[0];
+  pksrec.bandwidth = chanWidth * nChan;
+  pksrec.freqInc   = cMBrec.fqDelt[0];
+  pksrec.restFreq  = cMBrec.restFreq;
+
+  pksrec.tcal.resize(nPol);
   for (uInt ipol = 0; ipol < nPol; ipol++) {
-    MBrec.tcal(ipol) = cFITSMBrec.tcal[0][ipol];
-  }
-  MBrec.tcalTime  = trim(cFITSMBrec.tcalTime);
-  MBrec.azimuth   = cFITSMBrec.azimuth;
-  MBrec.elevation = cFITSMBrec.elevation;
-  MBrec.parAngle  = cFITSMBrec.parAngle;
-  MBrec.focusAxi  = cFITSMBrec.focusAxi;
-  MBrec.focusTan  = cFITSMBrec.focusTan;
-  MBrec.focusRot  = cFITSMBrec.focusRot;
-
-  MBrec.temperature = cFITSMBrec.temp;
-  MBrec.pressure    = cFITSMBrec.pressure;
-  MBrec.humidity    = cFITSMBrec.humidity;
-  MBrec.windSpeed   = cFITSMBrec.windSpeed;
-  MBrec.windAz      = cFITSMBrec.windAz;
-
-  MBrec.refBeam = cFITSMBrec.refBeam;
-  MBrec.beamNo  = cFITSMBrec.beamNo;
-
-  MBrec.direction.resize(2);
-  MBrec.direction(0) = cFITSMBrec.ra;
-  MBrec.direction(1) = cFITSMBrec.dec;
-
-  MBrec.scanRate.resize(2);
-  MBrec.scanRate(0)  = cFITSMBrec.raRate;
-  MBrec.scanRate(1)  = cFITSMBrec.decRate;
-  MBrec.rateAge      = cFITSMBrec.rateAge;
-  MBrec.rateson      = cFITSMBrec.rateson;
-
-  MBrec.tsys.resize(nPol);
-  MBrec.sigma.resize(nPol);
-  MBrec.calFctr.resize(nPol);
+    pksrec.tcal(ipol) = cMBrec.tcal[0][ipol];
+  }
+  pksrec.tcalTime  = trim(cMBrec.tcalTime);
+  pksrec.azimuth   = cMBrec.azimuth;
+  pksrec.elevation = cMBrec.elevation;
+  pksrec.parAngle  = cMBrec.parAngle;
+
+  pksrec.focusAxi  = cMBrec.focusAxi;
+  pksrec.focusTan  = cMBrec.focusTan;
+  pksrec.focusRot  = cMBrec.focusRot;
+
+  pksrec.temperature = cMBrec.temp;
+  pksrec.pressure    = cMBrec.pressure;
+  pksrec.humidity    = cMBrec.humidity;
+  pksrec.windSpeed   = cMBrec.windSpeed;
+  pksrec.windAz      = cMBrec.windAz;
+
+  pksrec.refBeam = cMBrec.refBeam;
+  pksrec.beamNo  = cMBrec.beamNo;
+
+  pksrec.direction.resize(2);
+  pksrec.direction(0) = cMBrec.ra;
+  pksrec.direction(1) = cMBrec.dec;
+  pksrec.pCode        = cMBrec.pCode;
+  pksrec.rateAge      = cMBrec.rateAge;
+  pksrec.scanRate.resize(2);
+  pksrec.scanRate(0)  = cMBrec.raRate;
+  pksrec.scanRate(1)  = cMBrec.decRate;
+  pksrec.paRate       = cMBrec.paRate;
+
+  pksrec.tsys.resize(nPol);
+  pksrec.sigma.resize(nPol);
+  pksrec.calFctr.resize(nPol);
   for (uInt ipol = 0; ipol < nPol; ipol++) {
-    MBrec.tsys(ipol)  = cFITSMBrec.tsys[0][ipol];
-    MBrec.sigma(ipol) = (MBrec.tsys(ipol) / 0.81) /
-                         sqrt(MBrec.interval * chanWidth);
-    MBrec.calFctr(ipol) = cFITSMBrec.calfctr[0][ipol];
-  }
-
-  if (cFITSMBrec.haveBase) {
-    MBrec.baseLin.resize(2,nPol);
-    MBrec.baseSub.resize(9,nPol);
+    pksrec.tsys(ipol)  = cMBrec.tsys[0][ipol];
+    pksrec.sigma(ipol) = (pksrec.tsys(ipol) / 0.81) /
+                            sqrt(pksrec.interval * chanWidth);
+    pksrec.calFctr(ipol) = cMBrec.calfctr[0][ipol];
+  }
+
+  if (cMBrec.haveBase) {
+    pksrec.baseLin.resize(2,nPol);
+    pksrec.baseSub.resize(9,nPol);
 
     for (uInt ipol = 0; ipol < nPol; ipol++) {
-      MBrec.baseLin(0,ipol) = cFITSMBrec.baseLin[0][ipol][0];
-      MBrec.baseLin(1,ipol) = cFITSMBrec.baseLin[0][ipol][1];
+      pksrec.baseLin(0,ipol) = cMBrec.baseLin[0][ipol][0];
+      pksrec.baseLin(1,ipol) = cMBrec.baseLin[0][ipol][1];
 
       for (uInt j = 0; j < 9; j++) {
-        MBrec.baseSub(j,ipol) = cFITSMBrec.baseSub[0][ipol][j];
+        pksrec.baseSub(j,ipol) = cMBrec.baseSub[0][ipol][j];
       }
     }
 
   } else {
-    MBrec.baseLin.resize(0,0);
-    MBrec.baseSub.resize(0,0);
-  }
-
-  if (cGetSpectra && cFITSMBrec.haveSpectra) {
-    MBrec.spectra.resize(nChan,nPol);
-    MBrec.spectra.takeStorage(IPosition(2,nChan,nPol), cFITSMBrec.spectra[0],
+    pksrec.baseLin.resize(0,0);
+    pksrec.baseSub.resize(0,0);
+  }
+
+  if (cGetSpectra && cMBrec.haveSpectra) {
+    pksrec.spectra.resize(nChan,nPol);
+    pksrec.spectra.takeStorage(IPosition(2,nChan,nPol), cMBrec.spectra[0],
       SHARE);
 
-    MBrec.flagged.resize(nChan,nPol);
-    MBrec.flagged.takeStorage(IPosition(2,nChan,nPol), cFITSMBrec.flagged[0],
+    pksrec.flagged.resize(nChan,nPol);
+    pksrec.flagged.takeStorage(IPosition(2,nChan,nPol), cMBrec.flagged[0],
       SHARE);
 
   } else {
-    MBrec.spectra.resize(0,0);
-    MBrec.flagged.resize(0,0);
+    pksrec.spectra.resize(0,0);
+    pksrec.flagged.resize(0,0);
   }
 
   if (cGetXPol) {
-    MBrec.xCalFctr = Complex(cFITSMBrec.xcalfctr[0][0],
-                             cFITSMBrec.xcalfctr[0][1]);
-    MBrec.xPol.resize(nChan);
-    MBrec.xPol.takeStorage(IPosition(1,nChan), (Complex *)cFITSMBrec.xpol[0],
+    pksrec.xCalFctr = Complex(cMBrec.xcalfctr[0][0],
+                             cMBrec.xcalfctr[0][1]);
+    pksrec.xPol.resize(nChan);
+    pksrec.xPol.takeStorage(IPosition(1,nChan), (Complex *)cMBrec.xpol[0],
       SHARE);
   }
@@ -451 +487 @@
 }
 
-//-------------------------------------------------------- PKSFITSreader::read
-
-// Read the next data record, just the basics.
-
-Int PKSFITSreader::read(
-        Int           &IFno,
-        Vector<Float> &tsys,
-        Vector<Float> &calFctr,
-        Matrix<Float> &baseLin,
-        Matrix<Float> &baseSub,
-        Matrix<Float> &spectra,
-        Matrix<uChar> &flagged)
-{
-  Int status;
-
-  if ((status = cReader->read(cFITSMBrec))) {
-    if (status != -1) {
-      status = 1;
-    }
-
-    return status;
-  }
-
-  IFno = cFITSMBrec.IFno[0];
-
-  uInt nChan = cFITSMBrec.nChan[0];
-  uInt nPol  = cFITSMBrec.nPol[0];
-
-  tsys.resize(nPol);
-  calFctr.resize(nPol);
-  for (uInt ipol = 0; ipol < nPol; ipol++) {
-    tsys(ipol) = cFITSMBrec.tsys[0][ipol];
-    calFctr(ipol) = cFITSMBrec.calfctr[0][ipol];
-  }
-
-  if (cFITSMBrec.haveBase) {
-    baseLin.resize(2,nPol);
-    baseSub.resize(9,nPol);
-
-    for (uInt ipol = 0; ipol < nPol; ipol++) {
-      baseLin(0,ipol) = cFITSMBrec.baseLin[0][ipol][0];
-      baseLin(1,ipol) = cFITSMBrec.baseLin[0][ipol][1];
-
-      for (uInt j = 0; j < 9; j++) {
-        baseSub(j,ipol) = cFITSMBrec.baseSub[0][ipol][j];
-      }
-    }
-
-  } else {
-    baseLin.resize(0,0);
-    baseSub.resize(0,0);
-  }
-
-  if (cGetSpectra && cFITSMBrec.haveSpectra) {
-    spectra.resize(nChan,nPol);
-    spectra.takeStorage(IPosition(2,nChan,nPol), cFITSMBrec.spectra[0],
-      SHARE);
-
-    flagged.resize(nChan,nPol);
-    flagged.takeStorage(IPosition(2,nChan,nPol), cFITSMBrec.flagged[0],
-      SHARE);
-
-  } else {
-    spectra.resize(0,0);
-    flagged.resize(0,0);
-  }
-
-  return 0;
-}
-
 //------------------------------------------------------- PKSFITSreader::close
 
@@ -528 +494 @@
 {
   cReader->close();
+  logMsg(cReader->getMsg());
+  cReader->clearMsg();
 }
 

trunk/external/atnf/PKSIO/PKSFITSreader.h

@@ -26 +26 @@
 //#                        Charlottesville, VA 22903-2475 USA
 //#
-//# $Id: PKSFITSreader.h,v 19.13 2008-06-26 02:02:43 cal103 Exp $
+//# $Id: PKSFITSreader.h,v 19.17 2008-11-17 06:38:05 cal103 Exp $
 //#---------------------------------------------------------------------------
 //# This class is basically a wrapper class for reading data from either an
@@ -39 +39 @@
 
 #include <atnf/PKSIO/FITSreader.h>
+#include <atnf/PKSIO/PKSrecord.h>
 #include <atnf/PKSIO/PKSreader.h>
 
 #include <casa/aips.h>
+#include <casa/stdio.h>
 #include <casa/Arrays/Vector.h>
 #include <casa/Arrays/Matrix.h>
@@ -47 +49 @@
 #include <casa/BasicSL/String.h>
 
+#include <casa/namespace.h>
+
 // <summary>
 // Class to read Parkes Multibeam data from a FITS file.
 // </summary>
 
-#include <casa/namespace.h>
 class PKSFITSreader : public PKSreader
 {
@@ -63 +66 @@
     // Destructor.
     virtual ~PKSFITSreader();
+
+    // Set message disposition.
+    virtual Int setMsg(
+        FILE *fd = 0x0);
 
     // Open the FITS file for reading.
@@ -104 +111 @@
         const Bool getSpectra = True,
         const Bool getXPol    = False,
-        const Bool getFeedPos = False);
+        const Int  coordSys   = 0);
 
     // Find the range of the data selected in time and position.
@@ -114 +121 @@
 
     // Read the next data record.
-    virtual Int read(MBrecord &mbrec);
-
-    // Read the next data record, just the basics.
-    virtual Int read(
-        Int           &IFno,
-        Vector<Float> &tsys,
-        Vector<Float> &calFctr,
-        Matrix<Float> &baseLin,
-        Matrix<Float> &baseSub,
-        Matrix<Float> &spectra,
-        Matrix<uChar> &flagged);
+    virtual Int read(PKSrecord &pksrec);
 
     // Close the FITS file.
@@ -132 +129 @@
     Int    *cBeams, *cIFs;
     uInt   cNBeam, cNIF;
-    PKSMBrecord cFITSMBrec;
+    MBrecord cMBrec;
     FITSreader  *cReader;
 

trunk/external/atnf/PKSIO/PKSMS2reader.cc

@@ -26 +26 @@
 //#                        Charlottesville, VA 22903-2475 USA
 //#
-//# $Id: PKSMS2reader.cc,v 19.13 2008-06-26 02:08:02 cal103 Exp $
+//# $Id: PKSMS2reader.cc,v 19.21 2008-11-17 06:55:18 cal103 Exp $
 //#---------------------------------------------------------------------------
 //# Original: 2000/08/03, Mark Calabretta, ATNF
 //#---------------------------------------------------------------------------
 
-
-// AIPS++ includes.
+#include <atnf/pks/pks_maths.h>
+#include <atnf/PKSIO/PKSmsg.h>
+#include <atnf/PKSIO/PKSrecord.h>
+#include <atnf/PKSIO/PKSMS2reader.h>
+
 #include <casa/stdio.h>
 #include <casa/Arrays/ArrayMath.h>
@@ -39 +42 @@
 #include <tables/Tables.h>
 
-// Parkes includes.
-#include <atnf/pks/pks_maths.h>
-#include <atnf/PKSIO/PKSMS2reader.h>
-
-
 //------------------------------------------------- PKSMS2reader::PKSMS2reader
 
@@ -51 +49 @@
 {
   cMSopen = False;
+
+  // By default, messages are written to stderr.
+  initMsg();
 }
 
@@ -353 +354 @@
         const Bool getSpectra,
         const Bool getXPol,
-        const Bool getFeedPos)
+        const Int  coordSys)
 {
   if (!cMSopen) {
@@ -432 +433 @@
   cGetXPol = cGetXPol && getXPol;
 
-  // Get feed positions?  (Not available.)
-  cGetFeedPos = False;
+  // Coordinate system?  (Only equatorial available.)
+  cCoordSys = 0;
 
   return maxNChan;
@@ -486 +487 @@
 // Read the next data record.
 
-Int PKSMS2reader::read(MBrecord &MBrec)
+Int PKSMS2reader::read(PKSrecord &pksrec)
 {
   if (!cMSopen) {
@@ -514 +515 @@
 
   // Renumerate scan no. Here still is 1-based
-  MBrec.scanNo = cScanNoCol(cIdx) - cScanNoCol(0) + 1;
-
-  if (MBrec.scanNo != cScanNo) {
+  pksrec.scanNo = cScanNoCol(cIdx) - cScanNoCol(0) + 1;
+
+  if (pksrec.scanNo != cScanNo) {
     // Start of new scan.
-    cScanNo  = MBrec.scanNo;
+    cScanNo  = pksrec.scanNo;
     cCycleNo = 1;
     cTime    = cTimeCol(cIdx);
@@ -524 +525 @@
 
   Double time = cTimeCol(cIdx);
-  MBrec.mjd      = time/86400.0;
-  MBrec.interval = cIntervalCol(cIdx);
+  pksrec.mjd      = time/86400.0;
+  pksrec.interval = cIntervalCol(cIdx);
 
   // Reconstruct the integration cycle number; due to small latencies the
   // integration time is usually slightly less than the time between cycles,
   // resetting cTime will prevent the difference from accumulating.
-  cCycleNo += nint((time - cTime)/MBrec.interval);
-  MBrec.cycleNo = cCycleNo;
-  cTime   = time;
+  cCycleNo += nint((time - cTime)/pksrec.interval);
+  pksrec.cycleNo = cCycleNo;
+  cTime = time;
 
   Int fieldId = cFieldIdCol(cIdx);
-  MBrec.fieldName = cFieldNameCol(fieldId);
+  pksrec.fieldName = cFieldNameCol(fieldId);
 
   Int srcId = cSrcIdCol(fieldId);
-  MBrec.srcName = cSrcNameCol(srcId);
-  MBrec.srcDir  = cSrcDirCol(srcId);
-  MBrec.srcPM   = cSrcPMCol(srcId);
+  pksrec.srcName = cSrcNameCol(srcId);
+  pksrec.srcDir  = cSrcDirCol(srcId);
+  pksrec.srcPM   = cSrcPMCol(srcId);
 
   // Systemic velocity.
   if (!cHaveSrcVel) {
-    MBrec.srcVel = 0.0f;
+    pksrec.srcVel = 0.0f;
   } else {
-    MBrec.srcVel  = cSrcVelCol(srcId)(IPosition(1,0));
+    pksrec.srcVel = cSrcVelCol(srcId)(IPosition(1,0));
   }
 
   // Observation type.
   Int stateId = cStateIdCol(cIdx);
-  MBrec.obsType = cObsModeCol(stateId);
-
-  MBrec.IFno = iIF + 1;
+  pksrec.obsType = cObsModeCol(stateId);
+
+  pksrec.IFno = iIF + 1;
   Int nChan = abs(cEndChan(iIF) - cStartChan(iIF)) + 1;
 
@@ -560 +561 @@
   if (nChan == 1) {
     cout << "The input is continuum data. "<< endl;
-    MBrec.freqInc  = chanFreq(0);
-    MBrec.refFreq  = chanFreq(0);
-    MBrec.restFreq = 0.0f;
+    pksrec.freqInc  = chanFreq(0);
+    pksrec.refFreq  = chanFreq(0);
+    pksrec.restFreq = 0.0f;
   } else {
     if (cStartChan(iIF) <= cEndChan(iIF)) {
-      MBrec.freqInc = chanFreq(1) - chanFreq(0);
+      pksrec.freqInc = chanFreq(1) - chanFreq(0);
     } else {
-      MBrec.freqInc = chanFreq(0) - chanFreq(1);
-    }
-
-    MBrec.refFreq  = chanFreq(cRefChan(iIF)-1);
-    MBrec.restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0));
-  }
-  MBrec.bandwidth = abs(MBrec.freqInc * nChan);
-
-  MBrec.tcal.resize(cNPol(iIF));
-  MBrec.tcal      = 0.0f;
-  MBrec.tcalTime  = "";
-  MBrec.azimuth   = 0.0f;
-  MBrec.elevation = 0.0f;
-  MBrec.parAngle  = 0.0f;
-  MBrec.focusAxi  = 0.0f;
-  MBrec.focusTan  = 0.0f;
-  MBrec.focusRot  = 0.0f;
+      pksrec.freqInc = chanFreq(0) - chanFreq(1);
+    }
+
+    pksrec.refFreq  = chanFreq(cRefChan(iIF)-1);
+    pksrec.restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0));
+  }
+  pksrec.bandwidth = abs(pksrec.freqInc * nChan);
+
+  pksrec.tcal.resize(cNPol(iIF));
+  pksrec.tcal      = 0.0f;
+  pksrec.tcalTime  = "";
+  pksrec.azimuth   = 0.0f;
+  pksrec.elevation = 0.0f;
+  pksrec.parAngle  = 0.0f;
+
+  pksrec.focusAxi  = 0.0f;
+  pksrec.focusTan  = 0.0f;
+  pksrec.focusRot  = 0.0f;
 
   // Find the appropriate entry in the WEATHER subtable.
@@ -596 +598 @@
   if (weatherIdx < 0) {
     // No appropriate WEATHER entry.
-    MBrec.temperature = 0.0f;
-    MBrec.pressure    = 0.0f;
-    MBrec.humidity    = 0.0f;
+    pksrec.temperature = 0.0f;
+    pksrec.pressure    = 0.0f;
+    pksrec.humidity    = 0.0f;
   } else {
-    MBrec.temperature = cTemperatureCol(weatherIdx);
-    MBrec.pressure    = cPressureCol(weatherIdx);
-    MBrec.humidity    = cHumidityCol(weatherIdx);
-  }
-
-  MBrec.windSpeed = 0.0f;
-  MBrec.windAz    = 0.0f;
-
-  MBrec.refBeam = 0;
-  MBrec.beamNo  = ibeam + 1;
+    pksrec.temperature = cTemperatureCol(weatherIdx);
+    pksrec.pressure    = cPressureCol(weatherIdx);
+    pksrec.humidity    = cHumidityCol(weatherIdx);
+  }
+
+  pksrec.windSpeed = 0.0f;
+  pksrec.windAz    = 0.0f;
+
+  pksrec.refBeam = 0;
+  pksrec.beamNo  = ibeam + 1;
 
   Matrix<Double> pointingDir = cPointingCol(fieldId);
-  MBrec.direction = pointingDir.column(0);
+  pksrec.direction = pointingDir.column(0);
+  pksrec.pCode = 0;
+  pksrec.rateAge = 0.0f;
   uInt ncols = pointingDir.ncolumn();
   if (ncols == 1) {
-    MBrec.scanRate = 0.0f;
+    pksrec.scanRate = 0.0f;
   } else {
-    MBrec.scanRate  = pointingDir.column(1);
-  }
-  MBrec.rateAge = 0;
-  MBrec.rateson = 0;
+    pksrec.scanRate(0) = pointingDir.column(1)(0);
+    pksrec.scanRate(1) = pointingDir.column(1)(1);
+  }
+  pksrec.paRate = 0.0f;
 
   // Get Tsys assuming that entries in the SYSCAL table match the main table.
@@ -630 +634 @@
   }
   if (cHaveTsys) {
-    cTsysCol.get(cIdx, MBrec.tsys, True);
+    cTsysCol.get(cIdx, pksrec.tsys, True);
   } else {
     Int numReceptor;
     cNumReceptorCol.get(0, numReceptor);
-    MBrec.tsys.resize(numReceptor);
-    MBrec.tsys = 1.0f;
-  }
-  cSigmaCol.get(cIdx, MBrec.sigma, True);
+    pksrec.tsys.resize(numReceptor);
+    pksrec.tsys = 1.0f;
+  }
+  cSigmaCol.get(cIdx, pksrec.sigma, True);
 
   // Calibration factors (if available).
-  MBrec.calFctr.resize(cNPol(iIF));
+  pksrec.calFctr.resize(cNPol(iIF));
   if (cHaveCalFctr) {
-    cCalFctrCol.get(cIdx, MBrec.calFctr);
+    cCalFctrCol.get(cIdx, pksrec.calFctr);
   } else {
-    MBrec.calFctr = 0.0f;
+    pksrec.calFctr = 0.0f;
   }
 
   // Baseline parameters (if available).
   if (cHaveBaseLin) {
-    MBrec.baseLin.resize(2,cNPol(iIF));
-    cBaseLinCol.get(cIdx, MBrec.baseLin);
-
-    MBrec.baseSub.resize(9,cNPol(iIF));
-    cBaseSubCol.get(cIdx, MBrec.baseSub);
+    pksrec.baseLin.resize(2,cNPol(iIF));
+    cBaseLinCol.get(cIdx, pksrec.baseLin);
+
+    pksrec.baseSub.resize(9,cNPol(iIF));
+    cBaseSubCol.get(cIdx, pksrec.baseSub);
 
   } else {
-    MBrec.baseLin.resize(0,0);
-    MBrec.baseSub.resize(0,0);
+    pksrec.baseLin.resize(0,0);
+    pksrec.baseSub.resize(0,0);
   }
 
@@ -670 +674 @@
     // Transpose spectra.
     Int nPol = tmpData.nrow();
-    MBrec.spectra.resize(nChan, nPol);
-    MBrec.flagged.resize(nChan, nPol);
+    pksrec.spectra.resize(nChan, nPol);
+    pksrec.flagged.resize(nChan, nPol);
     if (cEndChan(iIF) >= cStartChan(iIF)) {
       // Simple transposition.
       for (Int ipol = 0; ipol < nPol; ipol++) {
         for (Int ichan = 0; ichan < nChan; ichan++) {
-          MBrec.spectra(ichan,ipol) = tmpData(ipol,ichan);
-          MBrec.flagged(ichan,ipol) = tmpFlag(ipol,ichan);
+          pksrec.spectra(ichan,ipol) = tmpData(ipol,ichan);
+          pksrec.flagged(ichan,ipol) = tmpFlag(ipol,ichan);
         }
       }
@@ -686 +690 @@
       for (Int ipol = 0; ipol < nPol; ipol++) {
         for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
-          MBrec.spectra(ichan,ipol) = tmpData(ipol,jchan);
-          MBrec.flagged(ichan,ipol) = tmpFlag(ipol,jchan);
+          pksrec.spectra(ichan,ipol) = tmpData(ipol,jchan);
+          pksrec.flagged(ichan,ipol) = tmpFlag(ipol,jchan);
         }
       }
@@ -696 +700 @@
   if (cGetXPol) {
     if (cHaveXCalFctr) {
-      cXCalFctrCol.get(cIdx, MBrec.xCalFctr);
+      cXCalFctrCol.get(cIdx, pksrec.xCalFctr);
     } else {
-      MBrec.xCalFctr = Complex(0.0f, 0.0f);
-    }
-
-    cDataCol.get(cIdx, MBrec.xPol, True);
+      pksrec.xCalFctr = Complex(0.0f, 0.0f);
+    }
+
+    cDataCol.get(cIdx, pksrec.xPol, True);
 
     if (cEndChan(iIF) < cStartChan(iIF)) {
@@ -707 +711 @@
       Int jchan = nChan - 1;
       for (Int ichan = 0; ichan < nChan/2; ichan++, jchan--) {
-        ctmp = MBrec.xPol(ichan);
-        MBrec.xPol(ichan) = MBrec.xPol(jchan);
-        MBrec.xPol(jchan) = ctmp;
-      }
-    }
-  }
-
-  cIdx++;
-
-  return 0;
-}
-
-//--------------------------------------------------------- PKSMS2reader::read
-
-// Read the next data record, just the basics.
-
-Int PKSMS2reader::read(
-        Int           &IFno,
-        Vector<Float> &tsys,
-        Vector<Float> &calFctr,
-        Matrix<Float> &baseLin,
-        Matrix<Float> &baseSub,
-        Matrix<Float> &spectra,
-        Matrix<uChar> &flagged)
-{
-  if (!cMSopen) {
-    return 1;
-  }
-
-  // Check for EOF.
-  if (cIdx >= cNRow) {
-    return -1;
-  }
-
-  // Find the next selected beam and IF.
-  Int ibeam;
-  Int iIF;
-  while (True) {
-    ibeam = cBeamNoCol(cIdx);
-    iIF   = cDataDescIdCol(cIdx);
-    if (cBeams(ibeam) && cIFs(iIF)) {
-      break;
-    }
-
-    // Check for EOF.
-    if (++cIdx >= cNRow) {
-      return -1;
-    }
-  }
-
-  IFno = iIF + 1;
-
-  // Get Tsys assuming that entries in the SYSCAL table match the main table.
-  cTsysCol.get(cIdx, tsys, True);
-
-  // Calibration factors (if available).
-  if (cHaveCalFctr) {
-    cCalFctrCol.get(cIdx, calFctr, True);
-  } else {
-    calFctr.resize(cNPol(iIF));
-    calFctr = 0.0f;
-  }
-
-  // Baseline parameters (if available).
-  if (cHaveBaseLin) {
-    baseLin.resize(2,cNPol(iIF));
-    cBaseLinCol.get(cIdx, baseLin);
-
-    baseSub.resize(9,cNPol(iIF));
-    cBaseSubCol.get(cIdx, baseSub);
-
-  } else {
-    baseLin.resize(0,0);
-    baseSub.resize(0,0);
-  }
-
-  if (cGetSpectra) {
-    // Get spectral data.
-    Matrix<Float> tmpData;
-    Matrix<Bool>  tmpFlag;
-    cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
-    cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
-
-    // Transpose spectra.
-    Int nChan = tmpData.ncolumn();
-    Int nPol  = tmpData.nrow();
-    spectra.resize(nChan, nPol);
-    flagged.resize(nChan, nPol);
-    if (cEndChan(iIF) >= cStartChan(iIF)) {
-      // Simple transposition.
-      for (Int ipol = 0; ipol < nPol; ipol++) {
-        for (Int ichan = 0; ichan < nChan; ichan++) {
-          spectra(ichan,ipol) = tmpData(ipol,ichan);
-          flagged(ichan,ipol) = tmpFlag(ipol,ichan);
-        }
-      }
-
-    } else {
-      // Transpose with inversion.
-      Int jchan = nChan - 1;
-      for (Int ipol = 0; ipol < nPol; ipol++) {
-        for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
-          spectra(ichan,ipol) = tmpData(ipol,jchan);
-          flagged(ichan,ipol) = tmpFlag(ipol,jchan);
-        }
+        ctmp = pksrec.xPol(ichan);
+        pksrec.xPol(ichan) = pksrec.xPol(jchan);
+        pksrec.xPol(jchan) = ctmp;
       }
     }

trunk/external/atnf/PKSIO/PKSMS2reader.h

@@ -26 +26 @@
 //#                        Charlottesville, VA 22903-2475 USA
 //#
-//# $Id: PKSMS2reader.h,v 19.14 2008-06-26 02:03:51 cal103 Exp $
+//# $Id: PKSMS2reader.h,v 19.17 2008-11-17 06:38:53 cal103 Exp $
 //#---------------------------------------------------------------------------
 //# Original: 2000/08/03, Mark Calabretta, ATNF
@@ -35 +35 @@
 
 #include <atnf/PKSIO/PKSreader.h>
+#include <atnf/PKSIO/PKSrecord.h>
 
 #include <casa/aips.h>
@@ -46 +47 @@
 #include <tables/Tables/ScalarColumn.h>
 
+#include <casa/namespace.h>
+
 // <summary>
 // Class to read Parkes Multibeam data from a v2 MS.
 // </summary>
-
-#include <casa/namespace.h>
 
 class PKSMS2reader : public PKSreader
@@ -101 +102 @@
         const Bool getSpectra = True,
         const Bool getXPol    = False,
-        const Bool getFeedPos = False);
+        const Int  coordSys   = 0);
 
     // Find the range of the data selected in time and position.
@@ -111 +112 @@
 
     // Read the next data record.
-    virtual Int read(MBrecord &mbrec);
-
-    // Read the next data record, just the basics.
-    virtual Int read(
-        Int           &IFno,
-        Vector<Float> &tsys,
-        Vector<Float> &calFctr,
-        Matrix<Float> &baseLin,
-        Matrix<Float> &baseSub,
-        Matrix<Float> &spectra,
-        Matrix<uChar> &flagged);
+    virtual Int read(PKSrecord &pksrec);
 
     // Close the MS.

trunk/external/atnf/PKSIO/PKSMS2writer.cc

@@ -2 +2 @@
 //# PKSMS2writer.cc: Class to write Parkes multibeam data to a measurementset.
 //#---------------------------------------------------------------------------
-//# Copyright (C) 2000-2007
+//# Copyright (C) 2000-2008
 //# Associated Universities, Inc. Washington DC, USA.
 //#
@@ -26 +26 @@
 //#                        Charlottesville, VA 22903-2475 USA
 //#
-//# $Id: PKSMS2writer.cc,v 19.12 2007/11/12 03:37:56 cal103 Exp $
+//# $Id: PKSMS2writer.cc,v 19.14 2008-11-17 06:56:13 cal103 Exp $
 //#---------------------------------------------------------------------------
 
+#include <atnf/PKSIO/PKSrecord.h>
 #include <atnf/PKSIO/PKSMS2writer.h>
 
@@ -54 +55 @@
 PKSMS2writer::PKSMS2writer()
 {
+  cPKSMS = 0x0;
+
+  // By default, messages are written to stderr.
+  initMsg();
 }
 
@@ -84 +89 @@
         const Bool   haveBase)
 {
+  if (cPKSMS) {
+    logMsg("ERROR: Output MS already open, close it first.");
+    return 1;
+  }
+
+  // Clear the message stack.
+  clearMsg();
+
   // Open a MS table.
   TableDesc pksDesc = MS::requiredTableDesc();
@@ -369 +382 @@
 
 Int PKSMS2writer::write(
-        const Int             scanNo,
-        const Int             cycleNo,
-        const Double          mjd,
-        const Double          interval,
-        const String          fieldName,
-        const String          srcName,
-        const Vector<Double>  srcDir,
-        const Vector<Double>  srcPM,
-        const Double          srcVel,
-        const String          obsMode,
-        const Int             IFno,
-        const Double          refFreq,
-        const Double          bandwidth,
-        const Double          freqInc,
-        const Double          restFreq,
-        const Vector<Float>   tcal,
-        const String          tcalTime,
-        const Float           azimuth,
-        const Float           elevation,
-        const Float           parAngle,
-        const Float           focusAxi,
-        const Float           focusTan,
-        const Float           focusRot,
-        const Float           temperature,
-        const Float           pressure,
-        const Float           humidity,
-        const Float           windSpeed,
-        const Float           windAz,
-        const Int             refBeam,
-        const Int             beamNo,
-        const Vector<Double>  direction,
-        const Vector<Double>  scanRate,
-        const Vector<Float>   tsys,
-        const Vector<Float>   sigma,
-        const Vector<Float>   calFctr,
-        const Matrix<Float>   baseLin,
-        const Matrix<Float>   baseSub,
-        const Matrix<Float>   &spectra,
-        const Matrix<uChar>   &flagged,
-        const Complex         xCalFctr,
-        const Vector<Complex> &xPol)
+        const PKSrecord &pksrec)
 {
   // Extend the time range in the OBSERVATION subtable.
   Vector<Double> timerange(2);
   cObservationCols->timeRange().get(0, timerange);
-  Double time = mjd*86400.0;
+  Double time = pksrec.mjd*86400.0;
   if (timerange(0) == 0.0) {
     timerange(0) = time;
@@ -421 +394 @@
   cObservationCols->timeRange().put(0, timerange);
 
-  Int iIF = IFno - 1;
+  Int iIF = pksrec.IFno - 1;
   Int nChan = cNChan(iIF);
   Int nPol  = cNPol(iIF);
@@ -427 +400 @@
   // IFno is the 1-relative row number in the DATA_DESCRIPTION,
   // SPECTRAL_WINDOW, and POLARIZATION subtables.
-  if (Int(cDataDescription.nrow()) < IFno) {
+  if (Int(cDataDescription.nrow()) < pksrec.IFno) {
     // Add a new entry to each subtable.
-    addDataDescriptionEntry(IFno);
-    addSpectralWindowEntry(IFno, nChan, refFreq, bandwidth, freqInc);
-    addPolarizationEntry(IFno, nPol);
+    addDataDescriptionEntry(pksrec.IFno);
+    addSpectralWindowEntry(pksrec.IFno, nChan, pksrec.refFreq,
+      pksrec.bandwidth, pksrec.freqInc);
+    addPolarizationEntry(pksrec.IFno, nPol);
   }
 
   // Find or add the source to the SOURCE subtable.
-  Int srcId = addSourceEntry(srcName, srcDir, srcPM, restFreq, srcVel);
-
-  // Find or add the obsMode to the STATE subtable.
-  Int stateId = addStateEntry(obsMode);
+  Int srcId = addSourceEntry(pksrec.srcName, pksrec.srcDir, pksrec.srcPM,
+    pksrec.restFreq, pksrec.srcVel);
+
+  // Find or add the obsType to the STATE subtable.
+  Int stateId = addStateEntry(pksrec.obsType);
 
   // FIELD subtable.
-  Int fieldId = addFieldEntry(fieldName, time, direction, scanRate, srcId);
+  Vector<Double> scanRate(2);
+  scanRate(0) = pksrec.scanRate(0);
+  scanRate(1) = pksrec.scanRate(1);
+  Int fieldId = addFieldEntry(pksrec.fieldName, time, pksrec.direction,
+    scanRate, srcId);
 
   // POINTING subtable.
-  addPointingEntry(time, interval, fieldName, direction, scanRate);
+  addPointingEntry(time, pksrec.interval, pksrec.fieldName, pksrec.direction,
+    scanRate);
 
   // SYSCAL subtable.
-  addSysCalEntry(beamNo, iIF, time, interval, tcal, tsys);
+  addSysCalEntry(pksrec.beamNo, iIF, time, pksrec.interval, pksrec.tcal,
+    pksrec.tsys);
 
   // Handle weather information.
@@ -453 +434 @@
   Int nWeather = wTime.nrow();
   if (nWeather == 0 || time > wTime(nWeather-1)) {
-    addWeatherEntry(time, interval, pressure, humidity, temperature);
+    addWeatherEntry(time, pksrec.interval, pksrec.pressure, pksrec.humidity,
+      pksrec.temperature);
   }
 
@@ -465 +447 @@
   cMSCols->antenna1().put(irow, 0);
   cMSCols->antenna2().put(irow, 0);
-  cMSCols->feed1().put(irow, beamNo-1);
-  cMSCols->feed2().put(irow, beamNo-1);
+  cMSCols->feed1().put(irow, pksrec.beamNo-1);
+  cMSCols->feed2().put(irow, pksrec.beamNo-1);
   cMSCols->dataDescId().put(irow, iIF);
   cMSCols->processorId().put(irow, 0);
@@ -472 +454 @@
 
   // Non-key attributes.
-  cMSCols->interval().put(irow, interval);
-  cMSCols->exposure().put(irow, interval);
+  cMSCols->interval().put(irow, pksrec.interval);
+  cMSCols->exposure().put(irow, pksrec.interval);
   cMSCols->timeCentroid().put(irow, time);
-  cMSCols->scanNumber().put(irow, scanNo);
+  cMSCols->scanNumber().put(irow, pksrec.scanNo);
   cMSCols->arrayId().put(irow, 0);
   cMSCols->observationId().put(irow, 0);
@@ -485 +467 @@
   // Baseline fit parameters.
   if (cHaveBase) {
-    cBaseLinCol->put(irow, baseLin);
-
-    if (baseSub.nrow() == 9) {
-      cBaseSubCol->put(irow, baseSub);
+    cBaseLinCol->put(irow, pksrec.baseLin);
+
+    if (pksrec.baseSub.nrow() == 9) {
+      cBaseSubCol->put(irow, pksrec.baseSub);
 
     } else {
       Matrix<Float> tmp(9, 2, 0.0f);
       for (Int ipol = 0; ipol < nPol; ipol++) {
-        for (uInt j = 0; j < baseSub.nrow(); j++) {
-          tmp(j,ipol) = baseSub(j,ipol);
+        for (uInt j = 0; j < pksrec.baseSub.nrow(); j++) {
+          tmp(j,ipol) = pksrec.baseSub(j,ipol);
         }
       }
@@ -506 +488 @@
   for (Int ipol = 0; ipol < nPol; ipol++) {
     for (Int ichan = 0; ichan < nChan; ichan++) {
-      tmpData(ipol,ichan) = spectra(ichan,ipol);
-      tmpFlag(ipol,ichan) = flagged(ichan,ipol);
+      tmpData(ipol,ichan) = pksrec.spectra(ichan,ipol);
+      tmpFlag(ipol,ichan) = pksrec.flagged(ichan,ipol);
     }
   }
 
-  cCalFctrCol->put(irow, calFctr);
+  cCalFctrCol->put(irow, pksrec.calFctr);
   cMSCols->floatData().put(irow, tmpData);
   cMSCols->flag().put(irow, tmpFlag);
@@ -517 +499 @@
   // Cross-polarization spectra.
   if (cHaveXPol(iIF)) {
-    cXCalFctrCol->put(irow, xCalFctr);
-    cMSCols->data().put(irow, xPol);
-  }
-
-  cMSCols->sigma().put(irow, sigma);
+    cXCalFctrCol->put(irow, pksrec.xCalFctr);
+    cMSCols->data().put(irow, pksrec.xPol);
+  }
+
+  cMSCols->sigma().put(irow, pksrec.sigma);
 
   Vector<Float> weight(1, 1.0f);
@@ -586 +568 @@
   cSysCal          = MSSysCal();
   cWeather         = MSWeather();
+
   // Release the main table.
-  delete cPKSMS; cPKSMS=0;
+  delete cPKSMS;
+  cPKSMS = 0x0;
 }
 
@@ -984 +968 @@
 
 Int PKSMS2writer::addStateEntry(
-        const String obsMode)
+        const String obsType)
 {
   // Look for an entry in the STATE subtable.
   for (uInt n = 0; n < cStateCols->nrow(); n++) {
-    if (cStateCols->obsMode()(n) == obsMode) {
+    if (cStateCols->obsMode()(n) == obsType) {
       return n;
     }
@@ -998 +982 @@
 
   // Data.
-  if (obsMode.contains("RF")) {
+  if (obsType.contains("RF")) {
     cStateCols->sig().put(n, False);
     cStateCols->ref().put(n, True);
-  } else if (!obsMode.contains("PA")) {
+  } else if (!obsType.contains("PA")) {
     // Signal and reference are both false for "paddle" data.
     cStateCols->sig().put(n, True);
@@ -1010 +994 @@
   cStateCols->cal().put(n, 0.0);
   cStateCols->subScan().put(n, 0);
-  cStateCols->obsMode().put(n, obsMode);
+  cStateCols->obsMode().put(n, obsType);
 
   // Flags.

trunk/external/atnf/PKSIO/PKSMS2writer.h

@@ -26 +26 @@
 //#                        Charlottesville, VA 22903-2475 USA
 //#
-//# $Id: PKSMS2writer.h,v 19.12 2007/11/12 03:37:56 cal103 Exp $
+//# $Id: PKSMS2writer.h,v 19.13 2008-11-17 06:40:25 cal103 Exp $
 //#---------------------------------------------------------------------------
 
@@ -32 +32 @@
 #define ATNF_PKSMS2WRITER_H
 
+#include <atnf/PKSIO/PKSrecord.h>
 #include <atnf/PKSIO/PKSwriter.h>
 
@@ -38 +39 @@
 #include <casa/Arrays/Vector.h>
 #include <casa/BasicSL/Complex.h>
+#include <casa/BasicSL/String.h>
 #include <ms/MeasurementSets/MeasurementSet.h>
 #include <ms/MeasurementSets/MSColumns.h>
-#include <casa/BasicSL/String.h>
+
+#include <casa/namespace.h>
 
 // <summary>
@@ -46 +49 @@
 // </summary>
 
-#include <casa/namespace.h>
 class PKSMS2writer : public PKSwriter
 {
@@ -74 +76 @@
     // Write the next data record.
     virtual Int write(
-        const Int             scanNo,
-        const Int             cycleNo,
-        const Double          mjd,
-        const Double          interval,
-        const String          fieldName,
-        const String          srcName,
-        const Vector<Double>  srcDir,
-        const Vector<Double>  srcPM,
-        const Double          srcVel,
-        const String          obsMode,
-        const Int             IFno,
-        const Double          refFreq,
-        const Double          bandwidth,
-        const Double          freqInc,
-        const Double          restFreq,
-        const Vector<Float>   tcal,
-        const String          tcalTime,
-        const Float           azimuth,
-        const Float           elevation,
-        const Float           parAngle,
-        const Float           focusAxi,
-        const Float           focusTan,
-        const Float           focusRot,
-        const Float           temperature,
-        const Float           pressure,
-        const Float           humidity,
-        const Float           windSpeed,
-        const Float           windAz,
-        const Int             refBeam,
-        const Int             beamNo,
-        const Vector<Double>  direction,
-        const Vector<Double>  scanRate,
-        const Vector<Float>   tsys,
-        const Vector<Float>   sigma,
-        const Vector<Float>   calFctr,
-        const Matrix<Float>   baseLin,
-        const Matrix<Float>   baseSub,
-        const Matrix<Float>   &spectra,
-        const Matrix<uChar>   &flagged,
-        const Complex         xCalFctr,
-        const Vector<Complex> &xPol);
+        const PKSrecord &pksrec);
 
     // Close the MS, flushing all associated Tables.

trunk/external/atnf/PKSIO/PKSSDwriter.cc

@@ -2 +2 @@
 //# PKSSDwriter.cc: Class to write Parkes multibeam data to an SDFITS file.
 //#---------------------------------------------------------------------------
-//# Copyright (C) 2000-2007
+//# Copyright (C) 2000-2008
 //# Associated Universities, Inc. Washington DC, USA.
 //#
@@ -26 +26 @@
 //#                        Charlottesville, VA 22903-2475 USA
 //#
-//# $Id: PKSSDwriter.cc,v 19.13 2007/11/12 03:37:56 cal103 Exp $
+//# $Id: PKSSDwriter.cc,v 19.15 2008-11-17 06:56:50 cal103 Exp $
 //#---------------------------------------------------------------------------
 
-#include <atnf/PKSIO/PKSMBrecord.h>
+#include <atnf/PKSIO/MBrecord.h>
 #include <atnf/PKSIO/PKSSDwriter.h>
 
+#include <casa/stdio.h>
 #include <casa/Quanta/MVTime.h>
 
-
 //--------------------------------------------------- PKSSDwriter::PKSSDwriter
 
@@ -41 +41 @@
 PKSSDwriter::PKSSDwriter()
 {
+  // By default, messages are written to stderr.
+  initMsg();
 }
 
@@ -50 +52 @@
 {
   close();
+}
+
+//-------------------------------------------------------- PKSSDwriter::setMsg
+
+// Set message disposition.  If fd is non-zero messages will be written
+// to that file descriptor, else stored for retrieval by getMsg().
+
+Int PKSSDwriter::setMsg(FILE *fd)
+{
+  PKSmsg::setMsg(fd);
+  cSDwriter.setMsg(fd);
+
+  return 0;
 }
 
@@ -71 +86 @@
         const Bool   haveBase)
 {
+  // Clear the message stack.
+  clearMsg();
+
   double antPos[3];
   antPos[0] = antPosition(0);
@@ -107 +125 @@
         (int *)cNPol.getStorage(deleteIt),
         (int *)cHaveXPol.getStorage(deleteIt), (int)cHaveBase, 1);
+  logMsg(cSDwriter.getMsg());
+  cSDwriter.clearMsg();
   if (status) {
-    cSDwriter.reportError();
     cSDwriter.deleteFile();
     close();
@@ -121 +140 @@
 
 Int PKSSDwriter::write(
-        const Int             scanNo,
-        const Int             cycleNo,
-        const Double          mjd,
-        const Double          interval,
-        const String          fieldName,
-        const String          srcName,
-        const Vector<Double>  srcDir,
-        const Vector<Double>  srcPM,
-        const Double          srcVel,
-        const String          obsMode,
-        const Int             IFno,
-        const Double          refFreq,
-        const Double          bandwidth,
-        const Double          freqInc,
-        const Double          restFreq,
-        const Vector<Float>   tcal,
-        const String          tcalTime,
-        const Float           azimuth,
-        const Float           elevation,
-        const Float           parAngle,
-        const Float           focusAxi,
-        const Float           focusTan,
-        const Float           focusRot,
-        const Float           temperature,
-        const Float           pressure,
-        const Float           humidity,
-        const Float           windSpeed,
-        const Float           windAz,
-        const Int             refBeam,
-        const Int             beamNo,
-        const Vector<Double>  direction,
-        const Vector<Double>  scanRate,
-        const Vector<Float>   tsys,
-        const Vector<Float>   sigma,
-        const Vector<Float>   calFctr,
-        const Matrix<Float>   baseLin,
-        const Matrix<Float>   baseSub,
-        const Matrix<Float>   &spectra,
-        const Matrix<uChar>   &flagged,
-        const Complex         xCalFctr,
-        const Vector<Complex> &xPol)
+        const PKSrecord &pksrec)
 {
   // Do basic checks.
+  Int IFno = pksrec.IFno;
   uInt iIF = IFno - 1;
   if (IFno < 1 || Int(cNIF) < IFno) {
@@ -172 +152 @@
   }
 
-  uInt nChan = spectra.nrow();
+  uInt nChan = pksrec.spectra.nrow();
   if (nChan != cNChan(iIF)) {
     cerr << "PKSDwriter::write: "
@@ -181 +161 @@
   }
 
-  uInt nPol = spectra.ncolumn();
+  uInt nPol = pksrec.spectra.ncolumn();
   if (nPol != cNPol(iIF)) {
     cerr << "PKSDwriter::write: "
@@ -191 +171 @@
 
   // Extract calendar information from mjd.
-  MVTime time(mjd);
+  MVTime time(pksrec.mjd);
   Int year  = time.year();
   Int month = time.month();
   Int day   = time.monthday();
 
-  // Transfer data to a single-IF PKSMBrecord.
-  PKSMBrecord mbrec(1);
+  // Transfer data to a single-IF MBrecord.
+  MBrecord mbrec(1);
 
   // Start with basic beam- and IF-independent bookkeeping information.
-  mbrec.scanNo  = scanNo;
-  mbrec.cycleNo = cycleNo;
+  mbrec.scanNo  = pksrec.scanNo;
+  mbrec.cycleNo = pksrec.cycleNo;
 
   sprintf(mbrec.datobs, "%4.4d-%2.2d-%2.2d", year, month, day);
-  mbrec.utc      = fmod(mjd, 1.0) * 86400.0;
-  mbrec.exposure = float(interval);
-
-  strncpy(mbrec.srcName, (char *)srcName.chars(), 17);
-  mbrec.srcRA    = srcDir(0);
-  mbrec.srcDec   = srcDir(1);
-
-  mbrec.restFreq = restFreq;
-
-  strncpy(mbrec.obsType, (char *)obsMode.chars(), 16);
+  mbrec.utc      = fmod(pksrec.mjd, 1.0) * 86400.0;
+  mbrec.exposure = float(pksrec.interval);
+
+  strncpy(mbrec.srcName, (char *)pksrec.srcName.chars(), 17);
+  mbrec.srcRA    = pksrec.srcDir(0);
+  mbrec.srcDec   = pksrec.srcDir(1);
+
+  mbrec.restFreq = pksrec.restFreq;
+
+  strncpy(mbrec.obsType, (char *)pksrec.obsType.chars(), 16);
 
   // Now beam-dependent parameters.
-  mbrec.beamNo   = beamNo;
-  mbrec.ra       = direction(0);
-  mbrec.dec      = direction(1);
-  mbrec.raRate   = scanRate(0);
-  mbrec.decRate  = scanRate(1);
+  mbrec.beamNo   = pksrec.beamNo;
+  mbrec.ra       = pksrec.direction(0);
+  mbrec.dec      = pksrec.direction(1);
+  mbrec.raRate   = pksrec.scanRate(0);
+  mbrec.decRate  = pksrec.scanRate(1);
 
   // Now IF-dependent parameters.
@@ -229 +209 @@
 
   mbrec.fqRefPix[0] = (nChan/2) + 1;
-  mbrec.fqRefVal[0] = refFreq;
-  mbrec.fqDelt[0]   = freqInc;
+  mbrec.fqRefVal[0] = pksrec.refFreq;
+  mbrec.fqDelt[0]   = pksrec.freqInc;
 
   // Now the data itself.
-  for (uInt i = 0; i < tsys.nelements(); i++) {
-    mbrec.tsys[0][i] = tsys(i);
+  for (uInt i = 0; i < pksrec.tsys.nelements(); i++) {
+    mbrec.tsys[0][i] = pksrec.tsys(i);
   }
 
   for (uInt ipol = 0; ipol < nPol; ipol++) {
-    mbrec.calfctr[0][ipol] = calFctr(ipol);
+    mbrec.calfctr[0][ipol] = pksrec.calFctr(ipol);
   }
 
   if (cHaveXPol(iIF)) {
-    mbrec.xcalfctr[0][0] = xCalFctr.real();
-    mbrec.xcalfctr[0][1] = xCalFctr.imag();
+    mbrec.xcalfctr[0][0] = pksrec.xCalFctr.real();
+    mbrec.xcalfctr[0][1] = pksrec.xCalFctr.imag();
   } else {
     mbrec.xcalfctr[0][0] = 0.0f;
@@ -253 +233 @@
 
     for (uInt ipol = 0; ipol < nPol; ipol++) {
-      mbrec.baseLin[0][ipol][0] = baseLin(0,ipol);
-      mbrec.baseLin[0][ipol][1] = baseLin(1,ipol);
-
-      for (uInt j = 0; j < baseSub.nrow(); j++) {
-        mbrec.baseSub[0][ipol][j] = baseSub(j,ipol);
+      mbrec.baseLin[0][ipol][0] = pksrec.baseLin(0,ipol);
+      mbrec.baseLin[0][ipol][1] = pksrec.baseLin(1,ipol);
+
+      for (uInt j = 0; j < pksrec.baseSub.nrow(); j++) {
+        mbrec.baseSub[0][ipol][j] = pksrec.baseSub(j,ipol);
       }
-      for (uInt j = baseSub.nrow(); j < 9; j++) {
+      for (uInt j = pksrec.baseSub.nrow(); j < 9; j++) {
         mbrec.baseSub[0][ipol][j] = 0.0f;
       }
@@ -269 +249 @@
 
   Bool delSpectra = False;
-  const Float *specstor = spectra.getStorage(delSpectra);
+  const Float *specstor = pksrec.spectra.getStorage(delSpectra);
   mbrec.spectra[0] = (float *)specstor;
 
   Bool delFlagged = False;
-  const uChar *flagstor = flagged.getStorage(delFlagged);
+  const uChar *flagstor = pksrec.flagged.getStorage(delFlagged);
   mbrec.flagged[0] = (unsigned char *)flagstor;
 
@@ -279 +259 @@
   const Complex *xpolstor;
   if (cHaveXPol(iIF)) {
-    xpolstor = xPol.getStorage(delXPol);
+    xpolstor = pksrec.xPol.getStorage(delXPol);
   } else {
     xpolstor = 0;
@@ -287 +267 @@
   // Finish off with system calibration parameters.
   mbrec.extraSysCal = 1;
-  mbrec.refBeam     = refBeam;
-  for (uInt i = 0; i < tcal.nelements(); i++) {
-    mbrec.tcal[0][i] = tcal(i);
-  }
-  strncpy(mbrec.tcalTime, (char *)tcalTime.chars(), 16);
-  mbrec.azimuth   = azimuth;
-  mbrec.elevation = elevation;
-  mbrec.parAngle  = parAngle;
-  mbrec.focusAxi  = focusAxi;
-  mbrec.focusTan  = focusTan;
-  mbrec.focusRot  = focusRot;
-  mbrec.temp      = temperature;
-  mbrec.pressure  = pressure;
-  mbrec.humidity  = humidity;
-  mbrec.windSpeed = windSpeed;
-  mbrec.windAz    = windAz;
+  mbrec.refBeam     = pksrec.refBeam;
+  for (uInt i = 0; i < pksrec.tcal.nelements(); i++) {
+    mbrec.tcal[0][i] = pksrec.tcal(i);
+  }
+  strncpy(mbrec.tcalTime, (char *)pksrec.tcalTime.chars(), 16);
+  mbrec.azimuth   = pksrec.azimuth;
+  mbrec.elevation = pksrec.elevation;
+  mbrec.parAngle  = pksrec.parAngle;
+  mbrec.focusAxi  = pksrec.focusAxi;
+  mbrec.focusTan  = pksrec.focusTan;
+  mbrec.focusRot  = pksrec.focusRot;
+  mbrec.temp      = pksrec.temperature;
+  mbrec.pressure  = pksrec.pressure;
+  mbrec.humidity  = pksrec.humidity;
+  mbrec.windSpeed = pksrec.windSpeed;
+  mbrec.windAz    = pksrec.windAz;
 
   Int status = cSDwriter.write(mbrec);
+  logMsg(cSDwriter.getMsg());
+  cSDwriter.clearMsg();
   if (status) {
-    cSDwriter.reportError();
     status = 1;
   }
 
-  spectra.freeStorage(specstor, delSpectra);
-  flagged.freeStorage(flagstor, delFlagged);
-  xPol.freeStorage(xpolstor, delXPol);
+  pksrec.spectra.freeStorage(specstor, delSpectra);
+  pksrec.flagged.freeStorage(flagstor, delFlagged);
+  pksrec.xPol.freeStorage(xpolstor, delXPol);
 
   return status;
@@ -338 +319 @@
 {
   cSDwriter.close();
-}
+  logMsg(cSDwriter.getMsg());
+  cSDwriter.clearMsg();
+}

trunk/external/atnf/PKSIO/PKSSDwriter.h

@@ -26 +26 @@
 //#                        Charlottesville, VA 22903-2475 USA
 //#
-//# $Id: PKSSDwriter.h,v 19.14 2007/11/12 03:37:56 cal103 Exp $
+//# $Id: PKSSDwriter.h,v 19.16 2008-11-17 06:41:48 cal103 Exp $
 //# Original: 2000/07/21, Mark Calabretta, ATNF
 //#---------------------------------------------------------------------------
@@ -34 +34 @@
 
 #include <atnf/PKSIO/PKSwriter.h>
+#include <atnf/PKSIO/PKSrecord.h>
 #include <atnf/PKSIO/SDFITSwriter.h>
 
 #include <casa/aips.h>
+#include <casa/stdio.h>
 #include <casa/Arrays/Vector.h>
 #include <casa/Arrays/Matrix.h>
@@ -42 +44 @@
 #include <casa/BasicSL/String.h>
 
+#include <casa/namespace.h>
+
 // <summary>
 // Class to write Parkes multibeam data to an SDFITS file.
 // </summary>
 
-#include <casa/namespace.h>
 class PKSSDwriter : public PKSwriter
 {
@@ -55 +58 @@
     // Destructor.
     virtual ~PKSSDwriter();
+
+    // Set message disposition.
+    virtual Int setMsg(
+        FILE *fd = 0x0);
 
     // Create the SDFITS file and write static data.
@@ -74 +81 @@
     // Write the next data record.
     virtual Int write(
-        const Int             scanNo,
-        const Int             cycleNo,
-        const Double          mjd,
-        const Double          interval,
-        const String          fieldName,
-        const String          srcName,
-        const Vector<Double>  srcDir,
-        const Vector<Double>  srcPM,
-        const Double          srcVel,
-        const String          obsMode,
-        const Int             IFno,
-        const Double          refFreq,
-        const Double          bandwidth,
-        const Double          freqInc,
-        const Double          restFreq,
-        const Vector<Float>   tcal,
-        const String          tcalTime,
-        const Float           azimuth,
-        const Float           elevation,
-        const Float           parAngle,
-        const Float           focusAxi,
-        const Float           focusTan,
-        const Float           focusRot,
-        const Float           temperature,
-        const Float           pressure,
-        const Float           humidity,
-        const Float           windSpeed,
-        const Float           windAz,
-        const Int             refBeam,
-        const Int             beamNo,
-        const Vector<Double>  direction,
-        const Vector<Double>  scanRate,
-        const Vector<Float>   tsys,
-        const Vector<Float>   sigma,
-        const Vector<Float>   calFctr,
-        const Matrix<Float>   baselin,
-        const Matrix<Float>   basesub,
-        const Matrix<Float>   &spectra,
-        const Matrix<uChar>   &flagged,
-        const Complex         xCalFctr,
-        const Vector<Complex> &xPol);
+        const PKSrecord &pksrec);
 
     // Write a history record.

trunk/external/atnf/PKSIO/PKSreader.cc

@@ -26 +26 @@
 //#                        Charlottesville, VA 22903-2475 USA
 //#
-//# $Id: PKSreader.cc,v 19.6 2008-06-26 01:54:08 cal103 Exp $
+//# $Id: PKSreader.cc,v 19.12 2008-11-17 06:58:07 cal103 Exp $
 //#---------------------------------------------------------------------------
 //# Original: 2000/08/23, Mark Calabretta, ATNF
@@ -41 +41 @@
 #include <casa/OS/File.h>
 
-
 //--------------------------------------------------------------- getPKSreader
 
@@ -50 +49 @@
         const Int retry,
         const Int interpolate,
-        String &format,
-        Vector<Bool> &beams,
-        Vector<Bool> &IFs,
-        Vector<uInt> &nChan,
-        Vector<uInt> &nPol,
-        Vector<Bool> &haveXPol,
-        Bool   &haveBase,
-        Bool   &haveSpectra)
+        String &format)
 {
   // Check accessibility of the input.
@@ -63 +55 @@
   if (!inFile.exists()) {
     format = "DATASET NOT FOUND";
-    return 0;
+    return 0x0;
   }
 
   if (!inFile.isReadable()) {
     format = "DATASET UNREADABLE";
-    return 0;
+    return 0x0;
   }
 
   // Determine the type of input.
-  PKSreader *reader = 0;
+  PKSreader *reader = 0x0;
   if (inFile.isRegular()) {
     // Is it MBFITS or SDFITS?
@@ -112 +104 @@
   }
 
+  return reader;
+}
+
+//--------------------------------------------------------------- getPKSreader
+
+// Search a list of directories for a Parkes Multibeam dataset and return an
+// appropriate PKSreader for it.
+
+PKSreader* getPKSreader(
+        const String name,
+        const Vector<String> directories,
+        const Int retry,
+        const Int interpolate,
+        Int    &iDir,
+        String &format)
+{
+  PKSreader *reader = 0x0;
+
+  iDir = -1;
+  Int nDir = directories.nelements();
+  for (Int i = 0; i < nDir; i++) {
+    String inName = directories(i) + "/" + name;
+    reader = getPKSreader(inName, retry, interpolate, format);
+    if (reader) {
+      iDir = i;
+      break;
+    }
+  }
+
+  return reader;
+}
+
+//--------------------------------------------------------------- getPKSreader
+
+// Open an appropriate PKSreader for a Parkes Multibeam dataset.
+
+PKSreader* getPKSreader(
+        const String name,
+        const Int retry,
+        const Int interpolate,
+        String &format,
+        Vector<Bool> &beams,
+        Vector<Bool> &IFs,
+        Vector<uInt> &nChan,
+        Vector<uInt> &nPol,
+        Vector<Bool> &haveXPol,
+        Bool   &haveBase,
+        Bool   &haveSpectra)
+{
+  PKSreader *reader = getPKSreader(name, retry, interpolate, format);
 
   // Try to open it.
@@ -119 +161 @@
       format += " OPEN ERROR";
       delete reader;
-    } else {
-      return reader;
-    }
-  }
-
-  return 0;
-}
-
-
-//--------------------------------------------------------------- getPKSreader
-
-// Search a list of directories for a Parkes Multibeam dataset and return an
+      reader = 0x0;
+    }
+  }
+
+  return reader;
+}
+
+//--------------------------------------------------------------- getPKSreader
+
+// Search a list of directories for a Parkes Multibeam dataset and open an
 // appropriate PKSreader for it.
 
@@ -148 +188 @@
         Bool   &haveSpectra)
 {
-  Int nDir = directories.nelements();
-  for (iDir = 0; iDir < nDir; iDir++) {
-    String inName = directories(iDir) + "/" + name;
-    PKSreader *reader = getPKSreader(inName, retry, interpolate, format,
-                                     beams, IFs, nChan, nPol, haveXPol,
-                                     haveBase, haveSpectra);
-    if (reader != 0) {
-      return reader;
-    }
-  }
-
-  iDir = -1;
-  return 0;
-}
-
-
-//-------------------------------------------------------- PKSFITSreader::read
-
-// Read the next data record.
-
-Int PKSreader::read(
-        Int             &scanNo,
-        Int             &cycleNo,
-        Double          &mjd,
-        Double          &interval,
-        String          &fieldName,
-        String          &srcName,
-        Vector<Double>  &srcDir,
-        Vector<Double>  &srcPM,
-        Double          &srcVel,
-        String          &obsType,
-        Int             &IFno,
-        Double          &refFreq,
-        Double          &bandwidth,
-        Double          &freqInc,
-        Double          &restFreq,
-        Vector<Float>   &tcal,
-        String          &tcalTime,
-        Float           &azimuth,
-        Float           &elevation,
-        Float           &parAngle,
-        Float           &focusAxi,
-        Float           &focusTan,
-        Float           &focusRot,
-        Float           &temperature,
-        Float           &pressure,
-        Float           &humidity,
-        Float           &windSpeed,
-        Float           &windAz,
-        Int             &refBeam,
-        Int             &beamNo,
-        Vector<Double>  &direction,
-        Vector<Double>  &scanRate,
-        Vector<Float>   &tsys,
-        Vector<Float>   &sigma,
-        Vector<Float>   &calFctr,
-        Matrix<Float>   &baseLin,
-        Matrix<Float>   &baseSub,
-        Matrix<Float>   &spectra,
-        Matrix<uChar>   &flagged,
-        Complex         &xCalFctr,
-        Vector<Complex> &xPol)
-{
-  Int status;
-  MBrecord MBrec;
-
-  if ((status = read(MBrec))) {
-    if (status != -1) {
-      status = 1;
-    }
-
-    return status;
-  }
-
-  scanNo      = MBrec.scanNo;
-  cycleNo     = MBrec.cycleNo;
-  mjd         = MBrec.mjd;
-  interval    = MBrec.interval;
-  fieldName   = MBrec.fieldName;
-  srcName     = MBrec.srcName;
-  srcDir      = MBrec.srcDir;
-  srcPM       = MBrec.srcPM;
-  srcVel      = MBrec.srcVel;
-  obsType     = MBrec.obsType;
-  IFno        = MBrec.IFno;
-  refFreq     = MBrec.refFreq;
-  bandwidth   = MBrec.bandwidth;
-  freqInc     = MBrec.freqInc;
-  restFreq    = MBrec.restFreq;
-  tcal        = MBrec.tcal;
-  tcalTime    = MBrec.tcalTime;
-  azimuth     = MBrec.azimuth;
-  elevation   = MBrec.elevation;
-  parAngle    = MBrec.parAngle;
-  focusAxi    = MBrec.focusAxi;
-  focusTan    = MBrec.focusTan;
-  focusRot    = MBrec.focusRot;
-  temperature = MBrec.temperature;
-  pressure    = MBrec.pressure;
-  humidity    = MBrec.humidity;
-  windSpeed   = MBrec.windSpeed;
-  windAz      = MBrec.windAz;
-  refBeam     = MBrec.refBeam;
-  beamNo      = MBrec.beamNo;
-  direction   = MBrec.direction;
-  scanRate    = MBrec.scanRate;
-  tsys        = MBrec.tsys;
-  sigma       = MBrec.sigma;
-  calFctr     = MBrec.calFctr;
-  baseLin     = MBrec.baseLin;
-  baseSub     = MBrec.baseSub;
-  spectra     = MBrec.spectra;
-  flagged     = MBrec.flagged;
-  xCalFctr    = MBrec.xCalFctr;
-  xPol        = MBrec.xPol;
-
-  return 0;
-}
+  PKSreader *reader = getPKSreader(name, directories, retry, interpolate,
+                                   iDir, format);
+
+  // Try to open it.
+  if (reader) {
+    if (reader->open(name, beams, IFs, nChan, nPol, haveXPol, haveBase,
+                     haveSpectra)) {
+      format += " OPEN ERROR";
+      delete reader;
+      reader = 0x0;
+    }
+  }
+
+  return reader;
+}

trunk/external/atnf/PKSIO/PKSreader.h

@@ -26 +26 @@
 //#                        Charlottesville, VA 22903-2475 USA
 //#
-//# $Id: PKSreader.h,v 19.13 2008-06-26 01:50:24 cal103 Exp $
+//# $Id: PKSreader.h,v 19.22 2008-11-17 06:44:34 cal103 Exp $
 //#---------------------------------------------------------------------------
 //# Original: 2000/08/02, Mark Calabretta, ATNF
@@ -34 +34 @@
 #define ATNF_PKSREADER_H
 
+#include <atnf/PKSIO/PKSmsg.h>
+#include <atnf/PKSIO/PKSrecord.h>
+
 #include <casa/aips.h>
 #include <casa/Arrays/Matrix.h>
 #include <casa/Arrays/Vector.h>
-#include <casa/BasicSL/Complex.h>
 #include <casa/BasicSL/String.h>
 
@@ -45 +47 @@
 // Class to read Parkes multibeam data.
 // </summary>
+
+// Return an appropriate PKSreader for a Parkes Multibeam dataset.
+class PKSreader* getPKSreader(
+        const String name,
+        const Int retry,
+        const Int interpolate,
+        String &format);
+
+// As above, but search a list of directories for it.
+class PKSreader* getPKSreader(
+        const String name,
+        const Vector<String> directories,
+        const Int retry,
+        const Int interpolate,
+        Int    &iDir,
+        String &format);
 
 // Open an appropriate PKSreader for a Parkes Multibeam dataset.
@@ -76 +94 @@
         Bool   &haveSpectra);
 
-class MBrecord;
-
-class PKSreader
+class PKSreader : public PKSmsg
 {
   public:
@@ -116 +132 @@
     // Set data selection criteria.  Channel numbering is 1-relative, zero or
     // negative channel numbers are taken to be offsets from the last channel.
+    // Coordinate system selection (only supported for SDFITS input):
+    //   0: equatorial (RA,Dec),
+    //   1: vertical (Az,El),
+    //   2: feed-plane.
     virtual uInt select(
         const Vector<Bool> beamSel,
@@ -124 +144 @@
         const Bool getSpectra = True,
         const Bool getXPol    = False,
-        const Bool getFeedPos = False) = 0;
+        const Int  coordSys   = 0) = 0;
 
     // Find the range of the data selected in time and position.
@@ -133 +153 @@
         Matrix<Double> &positions) = 0;
 
-    // Read the next data record (MBrecord is defined below).
-    virtual Int read(MBrecord &mbrec) = 0;
-
-    // Read the next data record (for backwards compatibility, do not use).
-    virtual Int read(
-        Int             &scanNo,
-        Int             &cycleNo,
-        Double          &mjd,
-        Double          &interval,
-        String          &fieldName,
-        String          &srcName,
-        Vector<Double>  &srcDir,
-        Vector<Double>  &srcPM,
-        Double          &srcVel,
-        String          &obsType,
-        Int             &IFno,
-        Double          &refFreq,
-        Double          &bandwidth,
-        Double          &freqInc,
-        Double          &restFreq,
-        Vector<Float>   &tcal,
-        String          &tcalTime,
-        Float           &azimuth,
-        Float           &elevation,
-        Float           &parAngle,
-        Float           &focusAxi,
-        Float           &focusTan,
-        Float           &focusRot,
-        Float           &temperature,
-        Float           &pressure,
-        Float           &humidity,
-        Float           &windSpeed,
-        Float           &windAz,
-        Int             &refBeam,
-        Int             &beamNo,
-        Vector<Double>  &direction,
-        Vector<Double>  &scanRate,
-        Vector<Float>   &tsys,
-        Vector<Float>   &sigma,
-        Vector<Float>   &calFctr,
-        Matrix<Float>   &baseLin,
-        Matrix<Float>   &baseSub,
-        Matrix<Float>   &spectra,
-        Matrix<uChar>   &flagged,
-        Complex         &xCalFctr,
-        Vector<Complex> &xPol);
-
-    // Read the next data record, just the basics.
-    virtual Int read(
-        Int           &IFno,
-        Vector<Float> &tsys,
-        Vector<Float> &calFctr,
-        Matrix<Float> &baseLin,
-        Matrix<Float> &baseSub,
-        Matrix<Float> &spectra,
-        Matrix<uChar> &flagged) = 0;
+    // Read the next data record.
+    virtual Int read(PKSrecord &pksrec) = 0;
 
     // Close the input file.
@@ -194 +160 @@
 
   protected:
-    Bool   cGetFeedPos, cGetSpectra, cGetXPol;
+    Bool  cGetSpectra, cGetXPol;
+    Int   cCoordSys;
 
     Vector<uInt> cNChan, cNPol;
@@ -200 +167 @@
 };
 
-
-// Essentially just a struct used as a function argument.
-class MBrecord
-{
-  public:
-    Int             scanNo;
-    Int             cycleNo;
-    Double          mjd;
-    Double          interval;
-    String          fieldName;
-    String          srcName;
-    Vector<Double>  srcDir;
-    Vector<Double>  srcPM;
-    Double          srcVel;
-    String          obsType;
-    Int             IFno;
-    Double          refFreq;
-    Double          bandwidth;
-    Double          freqInc;
-    Double          restFreq;
-    Vector<Float>   tcal;
-    String          tcalTime;
-    Float           azimuth;
-    Float           elevation;
-    Float           parAngle;
-    Float           focusAxi;
-    Float           focusTan;
-    Float           focusRot;
-    Float           temperature;
-    Float           pressure;
-    Float           humidity;
-    Float           windSpeed;
-    Float           windAz;
-    Int             refBeam;
-    Int             beamNo;
-    Vector<Double>  direction;
-    Vector<Double>  scanRate;
-    Int             rateAge;
-    Int             rateson;
-    Vector<Float>   tsys;
-    Vector<Float>   sigma;
-    Vector<Float>   calFctr;
-    Matrix<Float>   baseLin;
-    Matrix<Float>   baseSub;
-    Matrix<Float>   spectra;
-    Matrix<uChar>   flagged;
-    Complex         xCalFctr;
-    Vector<Complex> xPol;
-};
-
 #endif

trunk/external/atnf/PKSIO/PKSwriter.h

@@ -26 +26 @@
 //#                        Charlottesville, VA 22903-2475 USA
 //#
-//# $Id: PKSwriter.h,v 19.14 2007/11/12 03:37:56 cal103 Exp $
+//# $Id: PKSwriter.h,v 19.16 2008-11-17 06:46:36 cal103 Exp $
 //#---------------------------------------------------------------------------
 
 #ifndef ATNF_PKSWRITER_H
 #define ATNF_PKSWRITER_H
+
+#include <atnf/PKSIO/PKSmsg.h>
+#include <atnf/PKSIO/PKSrecord.h>
 
 #include <casa/aips.h>
@@ -38 +41 @@
 #include <casa/BasicSL/String.h>
 
+#include <casa/namespace.h>
+
 // <summary>
 // Class to write out Parkes multibeam data.
 // </summary>
 
-#include <casa/namespace.h>
-class PKSwriter
+class PKSwriter : public PKSmsg
 {
   public:
@@ -67 +71 @@
     // Write the next data record.
     virtual Int write (
-        const Int             scanNo,
-        const Int             cycleNo,
-        const Double          mjd,
-        const Double          interval,
-        const String          fieldName,
-        const String          srcName,
-        const Vector<Double>  srcDir,
-        const Vector<Double>  srcPM,
-        const Double          srcVel,
-        const String          obsMode,
-        const Int             IFno,
-        const Double          refFreq,
-        const Double          bandwidth,
-        const Double          freqInc,
-        const Double          restFreq,
-        const Vector<Float>   tcal,
-        const String          tcalTime,
-        const Float           azimuth,
-        const Float           elevation,
-        const Float           parAngle,
-        const Float           focusAxi,
-        const Float           focusTan,
-        const Float           focusRot,
-        const Float           temperature,
-        const Float           pressure,
-        const Float           humidity,
-        const Float           windSpeed,
-        const Float           windAz,
-        const Int             refBeam,
-        const Int             beamNo,
-        const Vector<Double>  direction,
-        const Vector<Double>  scanRate,
-        const Vector<Float>   tsys,
-        const Vector<Float>   sigma,
-        const Vector<Float>   calFctr,
-        const Matrix<Float>   baseLin,
-        const Matrix<Float>   baseSub,
-        const Matrix<Float>   &spectra,
-        const Matrix<uChar>   &flagged,
-        const Complex         xCalFctr,
-        const Vector<Complex> &xPol) = 0;
+        const PKSrecord &pksrec) = 0;
 
     // Write a history record.

trunk/external/atnf/PKSIO/SDFITSreader.cc

@@ -26 +26 @@
 //#                        Charlottesville, VA 22903-2475 USA
 //#
-//# $Id: SDFITSreader.cc,v 19.24 2008-06-26 02:13:11 cal103 Exp $
+//# $Id: SDFITSreader.cc,v 19.33 2008-11-17 06:58:34 cal103 Exp $
 //#---------------------------------------------------------------------------
 //# The SDFITSreader class reads single dish FITS files such as those written
@@ -34 +34 @@
 //#---------------------------------------------------------------------------
 
+#include <atnf/pks/pks_maths.h>
+#include <atnf/PKSIO/PKSmsg.h>
+#include <atnf/PKSIO/MBrecord.h>
+#include <atnf/PKSIO/SDFITSreader.h>
+
+#include <casa/math.h>
+#include <casa/stdio.h>
+
 #include <algorithm>
 #include <strings.h>
-
-// AIPS++ includes.
-#include <casa/iostream.h>
-#include <casa/math.h>
-#include <casa/stdio.h>
-
-// ATNF includes.
-#include <atnf/pks/pks_maths.h>
-#include <atnf/PKSIO/PKSMBrecord.h>
-#include <atnf/PKSIO/SDFITSreader.h>
-
 
 class FITSparm
@@ -86 +83 @@
   cEndChan   = 0x0;
   cRefChan   = 0x0;
+
+  // By default, messages are written to stderr.
+  initMsg();
 }
 
@@ -114 +114 @@
         int    &extraSysCal)
 {
+  // Clear the message stack.
+  clearMsg();
+
   if (cSDptr) {
     close();
@@ -121 +124 @@
   cStatus = 0;
   if (fits_open_file(&cSDptr, sdName, READONLY, &cStatus)) {
-    cerr << "Failed to open SDFITS file: " << sdName << endl;
-    reportError();
+    sprintf(cMsg, "ERROR: Failed to open SDFITS file\n       %s", sdName);
+    logMsg(cMsg);
     return 1;
   }
@@ -139 +142 @@
         cALFA_CIMA = 1;
 
+        // Check for later versions of CIMAFITS.
+        float version;
+        readParm("VERSION", TFLOAT, &version);
+        if (version >= 2.0f) cALFA_CIMA = int(version);
+
       } else {
-        cerr << "Failed to locate SDFITS binary table." << endl;
-        reportError();
+        logMsg("ERROR: Failed to locate SDFITS binary table.");
         close();
         return 1;
@@ -177 +184 @@
     cNAxis = 5;
     if (readDim(DATA, 1, &cNAxis, cNAxes)) {
-      reportError();
+      logMsg();
       close();
       return 1;
@@ -196 +203 @@
     if (cNAxis < 4) {
       // Need at least four axes (for now).
-      cerr << "DATA array contains fewer than four axes." << endl;
+      logMsg("ERROR: DATA array contains fewer than four axes.");
       close();
       return 1;
     } else if (cNAxis > 5) {
       // We support up to five axes.
-      cerr << "DATA array contains more than five axes." << endl;
+      logMsg("ERROR: DATA array contains more than five axes.");
       close();
       return 1;
@@ -212 +219 @@
     findData(DATAXED, "DATAXED", TSTRING);
     if (cData[DATAXED].colnum < 0) {
-      cerr << "DATA array column absent from binary table." << endl;
+      logMsg("ERROR: DATA array column absent from binary table.");
       close();
       return 1;
@@ -242 +249 @@
 
   if (cStatus) {
-    reportError();
+    logMsg();
     close();
     return 1;
@@ -295 +302 @@
   for (int iaxis = 0; iaxis < 4; iaxis++) {
     if (cReqax[iaxis] < 0) {
-      cerr << "Could not find required DATA array axes." << endl;
+      logMsg("ERROR: Could not find required DATA array axes.");
       close();
       return 1;
@@ -345 +352 @@
 
   if (cStatus) {
-    reportError();
+    logMsg();
     close();
     return 1;
@@ -356 +363 @@
     cALFAscan = 0;
     cScanNo = 0;
-    if (cALFA_BD) {
+    if (cALFA_CIMA) {
+      findData(SCAN,  "SCAN_ID", TINT);
+      if (cALFA_CIMA > 1) {
+        findData(CYCLE, "RECNUM", TINT);
+      } else {
+        findData(CYCLE, "SUBSCAN", TINT);
+      }
+    } else if (cALFA_BD) {
       findData(SCAN,  "SCAN_NUMBER", TINT);
       findData(CYCLE, "PATTERN_NUMBER", TINT);
-    } else if (cALFA_CIMA) {
-      findData(SCAN,  "SCAN_ID", TINT);
-      findData(CYCLE, "SUBSCAN", TINT);
     }
   } else {
@@ -372 +383 @@
   // Beam number, 1-relative by default.
   cBeam_1rel = 1;
-  if (cData[BEAM].colnum < 0) {
+  if (cALFA) {
+    // ALFA INPUT_ID, 0-relative (overrides BEAM column if present).
+    findData(BEAM, "INPUT_ID", TSHORT);
+    cBeam_1rel = 0;
+
+  } else if (cData[BEAM].colnum < 0) {
     if (beamCRVAL) {
       // There is a BEAM axis.
       findData(BEAM, beamCRVAL, TDOUBLE);
-
     } else {
-      if (cALFA) {
-        // ALFA data, 0-relative.
-        findData(BEAM, "INPUT_ID", TSHORT);
-      } else {
-        // ms2sdfits output, 0-relative "feed" number.
-        findData(BEAM, "MAIN_FEED1", TSHORT);
-      }
-
+      // ms2sdfits output, 0-relative "feed" number.
+      findData(BEAM, "MAIN_FEED1", TSHORT);
       cBeam_1rel = 0;
     }
@@ -394 +403 @@
   if (cALFA && cData[IF].colnum < 0) {
     // ALFA data, 0-relative.
-    findData(IF, "IFVAL", TSHORT);
+    if (cALFA_CIMA > 1) {
+      findData(IF, "IFN", TSHORT);
+    } else {
+      findData(IF, "IFVAL", TSHORT);
+    }
     cIF_1rel = 0;
   }
@@ -477 +490 @@
   fits_get_num_rows(cSDptr, &cNRow, &cStatus);
   if (!cNRow) {
-    cerr << "Table contains no entries." << endl;
+    logMsg("ERROR: Table contains no entries.");
     close();
     return 1;
@@ -490 +503 @@
     if (fits_read_col(cSDptr, TSHORT, cData[BEAM].colnum, 1, 1, cNRow,
                       &beamNul, beamCol, &anynul, &cStatus)) {
-      reportError();
       delete [] beamCol;
+      logMsg();
       close();
       return 1;
@@ -505 +518 @@
       // Check validity.
       if (beamCol[irow] < cBeam_1rel) {
-        cerr << "SDFITS file contains invalid beam number." << endl;
         delete [] beamCol;
+        logMsg("ERROR: SDFITS file contains invalid beam number.");
         close();
         return 1;
@@ -546 +559 @@
     if (fits_read_col(cSDptr, TSHORT, cData[IF].colnum, 1, 1, cNRow,
                       &IFNul, IFCol, &anynul, &cStatus)) {
-      reportError();
       delete [] IFCol;
+      logMsg();
       close();
       return 1;
@@ -561 +574 @@
       // Check validity.
       if (IFCol[irow] < cIF_1rel) {
-        cerr << "SDFITS file contains invalid IF number." << endl;
         delete [] IFCol;
+        logMsg("ERROR: SDFITS file contains invalid IF number.");
         close();
         return 1;
@@ -594 +607 @@
             // Variable dimension array.
             if (readDim(DATA, irow+1, &cNAxis, cNAxes)) {
-              reportError();
+              logMsg();
               close();
               return 1;
@@ -622 +635 @@
 
           if (readDim(XPOLDATA, irow+1, &nAxis, nAxes)) {
-            reportError();
+            logMsg();
             close();
             return 1;
@@ -656 +669 @@
   }
 
-  if (cALFA) {
-    // ALFA labels each polarization as a separate IF.
+  if (cALFA && cALFA_CIMA < 2) {
+    // Older ALFA data labels each polarization as a separate IF.
     cNPol[0] = cNIF;
     cNIF = 1;
@@ -776 +789 @@
 
   // Brightness unit.
-  strcpy(bunit, cData[DATA].units);
+  if (cData[DATAXED].colnum >= 0) {
+    strcpy(bunit, "Jy");
+  } else {
+    strcpy(bunit, cData[DATA].units);
+  }
+
   if (strcmp(bunit, "JY") == 0) {
     bunit[1] = 'y';
@@ -836 +854 @@
 
   if (cStatus) {
-    reportError();
+    logMsg();
     return 1;
   }
@@ -849 +867 @@
 
   if (cStatus) {
-    reportError();
+    logMsg();
     return 1;
   }
@@ -902 +920 @@
     if (fits_read_col(cSDptr, TSHORT, cData[IF].colnum, 1, 1, cNRow,
                       &IFNul, IFCol, &anynul, &cStatus)) {
-      reportError();
       delete [] IFCol;
+      logMsg();
       close();
       return 1;
@@ -966 +984 @@
         double* &positions)
 {
-  int anynul;
-
   // Has the file been opened?
   if (!cSDptr) {
@@ -981 +997 @@
   }
 
+  int anynul;
   if (cData[BEAM].colnum > 0) {
     short *beamCol = new short[cNRow];
@@ -986 +1003 @@
     if (fits_read_col(cSDptr, TSHORT, cData[BEAM].colnum, 1, 1, cNRow,
                       &beamNul, beamCol, &anynul, &cStatus)) {
-      reportError();
       delete [] beamCol;
       delete [] sel;
+      logMsg();
       return 1;
     }
@@ -1006 +1023 @@
     if (fits_read_col(cSDptr, TSHORT, cData[IF].colnum, 1, 1, cNRow,
                       &IFNul, IFCol, &anynul, &cStatus)) {
-      reportError();
       delete [] IFCol;
       delete [] sel;
+      logMsg();
       return 1;
     }
@@ -1066 +1083 @@
 
   // Retrieve positions for selected data.
-  double *ra  = new double[cNRow];
-  double *dec = new double[cNRow];
-  fits_read_col(cSDptr, TDOUBLE, cData[RA].colnum,  1, 1, nRow, 0, ra,
-                &anynul, &cStatus);
-  fits_read_col(cSDptr, TDOUBLE, cData[DEC].colnum, 1, 1, nRow, 0, dec,
-                &anynul, &cStatus);
-
-  if (cALFA_BD) {
-    for (int irow = 0; irow < nRow; irow++) {
-      // Convert hours to degrees.
-      ra[irow] *= 15.0;
-    }
-  }
-
   int isel = 0;
   positions = new double[2*nSel];
 
-  // Parameters needed to compute feed-plane coordinates.
-  double *srcRA, *srcDec;
-  float  *par, *rot;
-  if (cGetFeedPos) {
-    srcRA  = new double[cNRow];
-    srcDec = new double[cNRow];
-    par    = new float[cNRow];
-    rot    = new float[cNRow];
-    fits_read_col(cSDptr, TDOUBLE, cData[OBJ_RA].colnum,   1, 1, nRow, 0,
-                  srcRA,  &anynul, &cStatus);
-    fits_read_col(cSDptr, TDOUBLE, cData[OBJ_DEC].colnum,  1, 1, nRow, 0,
-                  srcDec, &anynul, &cStatus);
-    fits_read_col(cSDptr, TFLOAT,  cData[PARANGLE].colnum, 1, 1, nRow, 0,
-                  par,    &anynul, &cStatus);
-    fits_read_col(cSDptr, TFLOAT,  cData[FOCUSROT].colnum, 1, 1, nRow, 0,
-                  rot,    &anynul, &cStatus);
-
-    for (int irow = 0; irow < nRow; irow++) {
-      if (sel[irow]) {
-        // Convert to feed-plane coordinates.
-        Double dist, pa;
-        distPA(ra[irow]*D2R, dec[irow]*D2R, srcRA[irow]*D2R, srcDec[irow]*D2R,
-               dist, pa);
-
-        Double spin = (par[irow] + rot[irow])*D2R - pa + PI;
-        if (spin > 2.0*PI) spin -= 2.0*PI;
-        Double squint = PI/2.0 - dist;
-
-        positions[isel++] = spin;
-        positions[isel++] = squint;
-      }
-    }
-
-    delete [] srcRA;
-    delete [] srcDec;
-    delete [] par;
-    delete [] rot;
+  if (cCoordSys == 1) {
+    // Vertical (Az,El).
+    if (cData[AZIMUTH].colnum  < 1 ||
+        cData[ELEVATIO].colnum < 1) {
+      logMsg("WARNING: Azimuth/elevation information absent.");
+      cStatus = -1;
+
+    } else {
+      float *az = new float[cNRow];
+      float *el = new float[cNRow];
+      fits_read_col(cSDptr, TFLOAT, cData[AZIMUTH].colnum,  1, 1, nRow, 0, az,
+                    &anynul, &cStatus);
+      fits_read_col(cSDptr, TFLOAT, cData[ELEVATIO].colnum, 1, 1, nRow, 0, el,
+                    &anynul, &cStatus);
+
+      if (!cStatus) {
+        for (int irow = 0; irow < nRow; irow++) {
+          if (sel[irow]) {
+            positions[isel++] = az[irow] * D2R;
+            positions[isel++] = el[irow] * D2R;
+          }
+        }
+      }
+
+      delete [] az;
+      delete [] el;
+    }
 
   } else {
-    for (int irow = 0; irow < nRow; irow++) {
-      if (sel[irow]) {
-        positions[isel++] =  ra[irow] * D2R;
-        positions[isel++] = dec[irow] * D2R;
-      }
-    }
-  }
-
+    double *ra  = new double[cNRow];
+    double *dec = new double[cNRow];
+    fits_read_col(cSDptr, TDOUBLE, cData[RA].colnum,  1, 1, nRow, 0, ra,
+                  &anynul, &cStatus);
+    fits_read_col(cSDptr, TDOUBLE, cData[DEC].colnum, 1, 1, nRow, 0, dec,
+                  &anynul, &cStatus);
+    if (cStatus) {
+      delete [] ra;
+      delete [] dec;
+      goto cleanup;
+    }
+
+    if (cALFA_BD) {
+      for (int irow = 0; irow < nRow; irow++) {
+        // Convert hours to degrees.
+        ra[irow] *= 15.0;
+      }
+    }
+
+    if (cCoordSys == 0) {
+      // Equatorial (RA,Dec).
+      for (int irow = 0; irow < nRow; irow++) {
+        if (sel[irow]) {
+          positions[isel++] =  ra[irow] * D2R;
+          positions[isel++] = dec[irow] * D2R;
+        }
+      }
+
+    } else if (cCoordSys == 2) {
+      // Feed-plane.
+      if (cData[OBJ_RA].colnum   < 0 ||
+          cData[OBJ_DEC].colnum  < 0 ||
+          cData[PARANGLE].colnum < 1 ||
+          cData[FOCUSROT].colnum < 1) {
+        logMsg("WARNING: Insufficient information to compute feed-plane\n"
+               "         coordinates.");
+        cStatus = -1;
+
+      } else {
+        double *srcRA  = new double[cNRow];
+        double *srcDec = new double[cNRow];
+        float  *par = new float[cNRow];
+        float  *rot = new float[cNRow];
+
+        if (cData[OBJ_RA].colnum == 0) {
+          // Header keyword.
+          readData(OBJ_RA, 0, srcRA);
+          for (int irow = 1; irow < nRow; irow++) {
+            srcRA[irow] = *srcRA;
+          }
+        } else {
+          // Table column.
+          fits_read_col(cSDptr, TDOUBLE, cData[OBJ_RA].colnum,   1, 1, nRow,
+                        0, srcRA,  &anynul, &cStatus);
+        }
+
+        if (cData[OBJ_DEC].colnum == 0) {
+          // Header keyword.
+          readData(OBJ_DEC, 0, srcDec);
+          for (int irow = 1; irow < nRow; irow++) {
+            srcDec[irow] = *srcDec;
+          }
+        } else {
+          // Table column.
+          fits_read_col(cSDptr, TDOUBLE, cData[OBJ_DEC].colnum,  1, 1, nRow,
+                        0, srcDec, &anynul, &cStatus);
+        }
+
+        fits_read_col(cSDptr, TFLOAT,  cData[PARANGLE].colnum, 1, 1, nRow, 0,
+                      par,    &anynul, &cStatus);
+        fits_read_col(cSDptr, TFLOAT,  cData[FOCUSROT].colnum, 1, 1, nRow, 0,
+                      rot,    &anynul, &cStatus);
+
+        if (!cStatus) {
+          for (int irow = 0; irow < nRow; irow++) {
+            if (sel[irow]) {
+              // Convert to feed-plane coordinates.
+              Double dist, pa;
+              distPA(ra[irow]*D2R, dec[irow]*D2R, srcRA[irow]*D2R,
+                     srcDec[irow]*D2R, dist, pa);
+
+              Double spin = (par[irow] + rot[irow])*D2R - pa + PI;
+              if (spin > 2.0*PI) spin -= 2.0*PI;
+              Double squint = PI/2.0 - dist;
+
+              positions[isel++] = spin;
+              positions[isel++] = squint;
+            }
+          }
+        }
+
+        delete [] srcRA;
+        delete [] srcDec;
+        delete [] par;
+        delete [] rot;
+      }
+    }
+
+    delete [] ra;
+    delete [] dec;
+  }
+
+cleanup:
   delete [] sel;
-  delete [] ra;
-  delete [] dec;
-
-  return cStatus;
+
+  if (cStatus) {
+    nSel = 0;
+    delete [] positions;
+    logMsg();
+    cStatus = 0;
+    return 1;
+  }
+
+  return 0;
 }
 
@@ -1143 +1237 @@
 
 int SDFITSreader::read(
-        PKSMBrecord &mbrec)
+        MBrecord &mbrec)
 {
   // Has the file been opened?
@@ -1175 +1269 @@
           readData(OBSMODE, cRow, chars);
           if (strcmp(chars, "CAL") == 0) {
-            // iIF is really the polarization in ALFA data.
-            alfaCal(iBeam, iIF);
-            continue;
+            if (cALFA_CIMA > 1) {
+              for (short iPol = 0; iPol < cNPol[iIF]; iPol++) {
+                alfaCal(iBeam, iIF, iPol);
+              }
+              continue;
+            } else {
+              // iIF is really the polarization in older ALFA data.
+              alfaCal(iBeam, 0, iIF);
+              continue;
+            }
           }
         }
@@ -1289 +1390 @@
     mbrec.decRate = scanrate[1] * D2R;
   }
-  mbrec.rateAge = 0;
-  mbrec.rateson = 0;
+  mbrec.paRate = 0.0f;
 
   // IF-dependent parameters.
@@ -1328 +1428 @@
 
   if (cStatus) {
-    reportError();
+    logMsg();
     return 1;
   }
@@ -1365 +1465 @@
 
   if (cStatus) {
-    reportError();
+    logMsg();
     return 1;
   }
@@ -1392 +1492 @@
       trc[cReqax[1]] = ipol+1;
 
-      if (cALFA) {
+      if (cALFA && cALFA_CIMA < 2) {
         // ALFA data: polarizations are stored in successive rows.
         blc[cReqax[1]] = 1;
@@ -1411 +1511 @@
 
         if ((status = readDim(DATA, cRow, &naxis, cNAxes))) {
-          reportError();
+          logMsg();
 
         } else if ((status = (naxis != cNAxis))) {
-          cerr << "DATA array dimensions changed." << endl;
+          logMsg("ERROR: DATA array dimensions changed.");
         }
 
@@ -1428 +1528 @@
           blc, trc, inc, 0, mbrec.spectra[0] + ipol*nChan, &anynul,
           &cStatus)) {
-        reportError();
+        logMsg();
         delete [] blc;
         delete [] trc;
@@ -1474 +1574 @@
             cNAxes, blc, trc, inc, 0, mbrec.flagged[0] + ipol*nChan, &anynul,
             &cStatus)) {
-          reportError();
+          logMsg();
           delete [] blc;
           delete [] trc;
@@ -1525 +1625 @@
     if (fits_read_subset_flt(cSDptr, cData[XPOLDATA].colnum, nAxis, nAxes,
         blc, trc, inc, 0, mbrec.xpol[0], &anynul, &cStatus)) {
-      reportError();
+      logMsg();
       delete [] blc;
       delete [] trc;
@@ -1556 +1656 @@
 
   if (cStatus) {
-    reportError();
+    logMsg();
     return 1;
   }
@@ -1564 +1664 @@
   readData(TCAL,     cRow, &mbrec.tcal[0]);
   readData(TCALTIME, cRow,  mbrec.tcalTime);
+
   readData(AZIMUTH,  cRow, &mbrec.azimuth);
   readData(ELEVATIO, cRow, &mbrec.elevation);
   readData(PARANGLE, cRow, &mbrec.parAngle);
+
   readData(FOCUSAXI, cRow, &mbrec.focusAxi);
   readData(FOCUSTAN, cRow, &mbrec.focusTan);
   readData(FOCUSROT, cRow, &mbrec.focusRot);
+
   readData(TAMBIENT, cRow, &mbrec.temp);
   readData(PRESSURE, cRow, &mbrec.pressure);
@@ -1588 +1691 @@
 
   if (cStatus) {
-    reportError();
+    logMsg();
     return 1;
   }
 
   return 0;
-}
-
-
-//------------------------------------------------------ SDFITSreader::alfaCal
-
-// Process ALFA calibration data.
-
-int SDFITSreader::alfaCal(
-        short iBeam,
-        short iPol)
-{
-  int  calOn;
-  char chars[32];
-  if (cALFA_BD) {
-    readData("OBS_NAME", TSTRING, cRow, chars);
-  } else {
-    readData("SCANTYPE", TSTRING, cRow, chars);
-  }
-
-  if (strcmp(chars, "ON") == 0) {
-    calOn = 1;
-  } else if (strcmp(chars, "OFF") == 0) {
-    calOn = 0;
-  } else {
-    return 1;
-  }
-
-  // Read cal data.
-  long *blc = new long[cNAxis+1];
-  long *trc = new long[cNAxis+1];
-  long *inc = new long[cNAxis+1];
-  for (int iaxis = 0; iaxis <= cNAxis; iaxis++) {
-    blc[iaxis] = 1;
-    trc[iaxis] = 1;
-    inc[iaxis] = 1;
-  }
-
-  // User channel selection.
-  int startChan = cStartChan[0];
-  int endChan   = cEndChan[0];
-
-  blc[cNAxis] = cRow;
-  trc[cNAxis] = cRow;
-  blc[cReqax[0]] = std::min(startChan, endChan);
-  trc[cReqax[0]] = std::max(startChan, endChan);
-  blc[cReqax[1]] = 1;
-  trc[cReqax[1]] = 1;
-
-  float spectrum[endChan];
-  int anynul;
-  if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAxis, cNAxes,
-      blc, trc, inc, 0, spectrum, &anynul, &cStatus)) {
-    reportError();
-    delete [] blc;
-    delete [] trc;
-    delete [] inc;
-    return 1;
-  }
-
-  // Average the spectrum.
-  float mean = 1e9f;
-  for (int k = 0; k < 2; k++) {
-    float discrim = 2.0f * mean;
-
-    int nChan = 0;
-    float sum = 0.0f;
-
-    float *chanN = spectrum + abs(endChan - startChan) + 1;
-    for (float *chan = spectrum; chan < chanN; chan++) {
-      // Simple discriminant that eliminates strong radar interference.
-      if (*chan < discrim) {
-        nChan++;
-        sum += *chan;
-      }
-    }
-
-    mean = sum / nChan;
-  }
-
-  if (calOn) {
-    cALFAcalOn[iBeam][iPol]  += mean;
-  } else {
-    cALFAcalOff[iBeam][iPol] += mean;
-  }
-
-  if (cALFAcalOn[iBeam][iPol] != 0.0f &&
-      cALFAcalOff[iBeam][iPol] != 0.0f) {
-    // Tcal should come from the TCAL table, it varies weakly with beam,
-    // polarization, and frequency.  However, TCAL is not written properly.
-    float Tcal = 12.0f;
-    cALFAcal[iBeam][iPol] = Tcal / (cALFAcalOn[iBeam][iPol] -
-                                    cALFAcalOff[iBeam][iPol]);
-
-    // Scale from K to Jy; the gain also varies weakly with beam,
-    // polarization, frequency, and zenith angle.
-    float fluxCal = 10.0f;
-    cALFAcal[iBeam][iPol] /= fluxCal;
-
-    cALFAcalOn[iBeam][iPol]  = 0.0f;
-    cALFAcalOff[iBeam][iPol] = 0.0f;
-  }
-
-  return 0;
-}
-
-
-//-------------------------------------------------- SDFITSreader::reportError
-
-// Print the error message corresponding to the input status value and all the
-// messages on the CFITSIO error stack to stderr.
-
-void SDFITSreader::reportError()
-{
-  fits_report_error(stderr, cStatus);
 }
 
@@ -1725 +1714 @@
     if (cEndChan)   delete [] cEndChan;
     if (cRefChan)   delete [] cRefChan;
+  }
+}
+
+//------------------------------------------------------- SDFITSreader::logMsg
+
+// Log a message.  If the current CFITSIO status value is non-zero, also log
+// the corresponding error message and the CFITSIO message stack.
+
+void SDFITSreader::logMsg(const char *msg)
+{
+  FITSreader::logMsg(msg);
+
+  if (cStatus > 0) {
+    fits_get_errstatus(cStatus, cMsg);
+    FITSreader::logMsg(cMsg);
+
+    while (fits_read_errmsg(cMsg)) {
+      FITSreader::logMsg(cMsg);
+    }
   }
 }
@@ -2011 +2019 @@
   }
 }
+
+//------------------------------------------------------ SDFITSreader::alfaCal
+
+// Process ALFA calibration data.
+
+int SDFITSreader::alfaCal(
+        short iBeam,
+        short iIF,
+        short iPol)
+{
+  int  calOn;
+  char chars[32];
+  if (cALFA_BD) {
+    readData("OBS_NAME", TSTRING, cRow, chars);
+  } else {
+    readData("SCANTYPE", TSTRING, cRow, chars);
+  }
+
+  if (strcmp(chars, "ON") == 0) {
+    calOn = 1;
+  } else if (strcmp(chars, "OFF") == 0) {
+    calOn = 0;
+  } else {
+    return 1;
+  }
+
+  // Read cal data.
+  long *blc = new long[cNAxis+1];
+  long *trc = new long[cNAxis+1];
+  long *inc = new long[cNAxis+1];
+  for (int iaxis = 0; iaxis <= cNAxis; iaxis++) {
+    blc[iaxis] = 1;
+    trc[iaxis] = 1;
+    inc[iaxis] = 1;
+  }
+
+  // User channel selection.
+  int startChan = cStartChan[iIF];
+  int endChan   = cEndChan[iIF];
+
+  blc[cNAxis] = cRow;
+  trc[cNAxis] = cRow;
+  blc[cReqax[0]] = std::min(startChan, endChan);
+  trc[cReqax[0]] = std::max(startChan, endChan);
+  if (cALFA_CIMA > 1) {
+    // CIMAFITS 2.x has a legitimate STOKES axis...
+    blc[cReqax[1]] = iPol+1;
+    trc[cReqax[1]] = iPol+1;
+  } else {
+    // ...older ALFA data does not.
+    blc[cReqax[1]] = 1;
+    trc[cReqax[1]] = 1;
+  }
+
+  float spectrum[endChan];
+  int anynul;
+  if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAxis, cNAxes,
+      blc, trc, inc, 0, spectrum, &anynul, &cStatus)) {
+    logMsg();
+    delete [] blc;
+    delete [] trc;
+    delete [] inc;
+    return 1;
+  }
+
+  // Average the spectrum.
+  float mean = 1e9f;
+  for (int k = 0; k < 2; k++) {
+    float discrim = 2.0f * mean;
+
+    int nChan = 0;
+    float sum = 0.0f;
+
+    float *chanN = spectrum + abs(endChan - startChan) + 1;
+    for (float *chan = spectrum; chan < chanN; chan++) {
+      // Simple discriminant that eliminates strong radar interference.
+      if (*chan < discrim) {
+        nChan++;
+        sum += *chan;
+      }
+    }
+
+    mean = sum / nChan;
+  }
+
+  if (calOn) {
+    cALFAcalOn[iBeam][iPol]  += mean;
+  } else {
+    cALFAcalOff[iBeam][iPol] += mean;
+  }
+
+  if (cALFAcalOn[iBeam][iPol] != 0.0f &&
+      cALFAcalOff[iBeam][iPol] != 0.0f) {
+    // Tcal should come from the TCAL table, it varies weakly with beam,
+    // polarization, and frequency.  However, TCAL is not written properly.
+    float Tcal = 12.0f;
+    cALFAcal[iBeam][iPol] = Tcal / (cALFAcalOn[iBeam][iPol] -
+                                    cALFAcalOff[iBeam][iPol]);
+
+    // Scale from K to Jy; the gain also varies weakly with beam,
+    // polarization, frequency, and zenith angle.
+    float fluxCal = 10.0f;
+    cALFAcal[iBeam][iPol] /= fluxCal;
+
+    cALFAcalOn[iBeam][iPol]  = 0.0f;
+    cALFAcalOff[iBeam][iPol] = 0.0f;
+  }
+
+  return 0;
+}

trunk/external/atnf/PKSIO/SDFITSreader.h

@@ -2 +2 @@
 //# SDFITSreader.h: ATNF CFITSIO interface class for SDFITS input.
 //#---------------------------------------------------------------------------
-//# Copyright (C) 2000-2007
+//# Copyright (C) 2000-2008
 //# Associated Universities, Inc. Washington DC, USA.
 //#
@@ -26 +26 @@
 //#                        Charlottesville, VA 22903-2475 USA
 //#
-//# $Id: SDFITSreader.h,v 19.13 2007/11/12 03:37:56 cal103 Exp $
+//# $Id: SDFITSreader.h,v 19.16 2008-11-17 06:47:05 cal103 Exp $
 //#---------------------------------------------------------------------------
 //# The SDFITSreader class reads single dish FITS files such as those written
@@ -38 +38 @@
 
 #include <atnf/PKSIO/FITSreader.h>
-#include <atnf/PKSIO/PKSMBrecord.h>
+#include <atnf/PKSIO/MBrecord.h>
 
 #include <fitsio.h>
+
+using namespace std;
 
 // <summary>
@@ -100 +102 @@
 
     // Read the next data record.
-    virtual int read(PKSMBrecord &record);
-
-    // Print out CFITSIO error messages.
-    void reportError(void);
+    virtual int read(MBrecord &record);
 
     // Close the SDFITS file.
@@ -125 +124 @@
           HUMIDITY, WINDSPEE, WINDDIRE, NDATA};
 
+    // Message handling.
+    virtual void logMsg(const char *msg = 0x0);
+
     void findData(int iData, char *name, int type);
     int   readDim(int iData, long iRow, int *naxis, long naxes[]);
@@ -132 +134 @@
     void  findCol(char *name, int *colnum);
 
-    // These are for ALFA data, "BDFITS" or "CIMAFITS".
+    // These are for ALFA data: "BDFITS" or "CIMAFITS".
     int   cALFA, cALFA_BD, cALFA_CIMA, cALFAscan, cScanNo;
     float cALFAcal[8][2], cALFAcalOn[8][2], cALFAcalOff[8][2];
-    int   alfaCal(short iBeam, short iIF);
+    int   alfaCal(short iBeam, short iIF, short iPol);
 
     // These are for GBT data.

trunk/external/atnf/PKSIO/SDFITSwriter.cc

@@ -2 +2 @@
 //# SDFITSwriter.cc: ATNF CFITSIO interface class for SDFITS output.
 //#---------------------------------------------------------------------------
-//# Copyright (C) 2000-2007
+//# Copyright (C) 2000-2008
 //# Mark Calabretta, ATNF
 //#
@@ -27 +27 @@
 //#                        AUSTRALIA
 //#
-//# $Id: SDFITSwriter.cc,v 19.12 2007/11/12 03:37:56 cal103 Exp $
+//# $Id: SDFITSwriter.cc,v 19.15 2008-11-17 06:58:58 cal103 Exp $
 //#---------------------------------------------------------------------------
 //# Original: 2000/07/24, Mark Calabretta, ATNF
 //#---------------------------------------------------------------------------
 
+#include <atnf/PKSIO/MBrecord.h>
+#include <atnf/PKSIO/PKSmsg.h>
+#include <atnf/PKSIO/SDFITSwriter.h>
+
+#include <casa/iostream.h>
+
 #include <algorithm>
 #include <math.h>
-
-// AIPS++ includes.
-#include <casa/iostream.h>
-
-// ATNF includes.
-#include <atnf/PKSIO/PKSMBrecord.h>
-#include <atnf/PKSIO/SDFITSwriter.h>
 
 using namespace std;
@@ -55 +54 @@
 {
   // Default constructor.
-  cSDptr = 0;
+  cSDptr = 0x0;
+
+  // By default, messages are written to stderr.
+  initMsg();
 }
 
@@ -86 +88 @@
         int    extraSysCal)
 {
+  if (cSDptr) {
+    logMsg("ERROR: Output file already open, close it first.");
+    return 1;
+  }
+
+  // Clear the message stack.
+  clearMsg();
+
   // Prepend an '!' to the output name to force it to be overwritten.
   char sdname[80];
@@ -94 +104 @@
   cStatus = 0;
   if (fits_create_file(&cSDptr, sdname, &cStatus)) {
+    sprintf(cMsg, "ERROR: Failed to create SDFITS file\n       %s", sdName);
+    logMsg(cMsg);
     return cStatus;
   }
@@ -141 +153 @@
   // Write required primary header keywords.
   if (fits_write_imghdr(cSDptr, 8, 0, 0, &cStatus)) {
+    logMsg("ERROR: Failed to write required primary header keywords.");
     return cStatus;
   }
@@ -160 +173 @@
   char version[7];
   char date[11];
-  sscanf("$Revision: 19.12 $", "%*s%s", version);
-  sscanf("$Date: 2007/11/12 03:37:56 $", "%*s%s", date);
+  sscanf("$Revision: 19.15 $", "%*s%s", version);
+  sscanf("$Date: 2008-11-17 06:58:58 $", "%*s%s", date);
   sprintf(text, "SDFITSwriter (v%s, %s)", version, date);
   fits_write_key_str(cSDptr, "ORIGIN", text, "output class", &cStatus);
@@ -171 +184 @@
   fits_write_comment(cSDptr, text, &cStatus);
 
+  if (cStatus) {
+    logMsg("ERROR: Failed in writing primary header.");
+    return cStatus;
+  }
+
+
   // Create an SDFITS extension.
   long nrow = 0;
@@ -176 +195 @@
   if (fits_create_tbl(cSDptr, BINARY_TBL, nrow, ncol, NULL, NULL, NULL,
       "SINGLE DISH", &cStatus)) {
+    logMsg("ERROR: Failed to create a binary table extension.");
     return 1;
   }
@@ -524 +544 @@
   }
 
+  if (cStatus) {
+    logMsg("ERROR: Failed in writing binary table header.");
+  }
+
   return cStatus;
 }
@@ -531 +555 @@
 // Write a record to the SDFITS file.
 
-int SDFITSwriter::write(PKSMBrecord &mbrec)
+int SDFITSwriter::write(MBrecord &mbrec)
 {
   char *cptr;
@@ -740 +764 @@
   }
 
+  if (cStatus) {
+    logMsg("ERROR: Failed in writing binary table entry.");
+  }
+
   return cStatus;
 }
@@ -751 +779 @@
 
 {
-  fits_write_history(cSDptr, text, &cStatus);
+  if (!cSDptr) {
+    return 1;
+  }
+
+  if (fits_write_history(cSDptr, text, &cStatus)) {
+    logMsg("ERROR: Failed in writing HISTORY records.");
+  }
 
   return cStatus;
-}
-
-//-------------------------------------------------- SDFITSwriter::reportError
-
-// Print the error message corresponding to the input status value and all the
-// messages on the CFITSIO error stack to stderr.
-
-void SDFITSwriter::reportError()
-{
-  fits_report_error(stderr, cStatus);
 }
 
@@ -774 +798 @@
   if (cSDptr) {
     cStatus = 0;
-    fits_close_file(cSDptr, &cStatus);
+    if (fits_close_file(cSDptr, &cStatus)) {
+      logMsg("ERROR: Failed to close file.");
+    }
+
     cSDptr = 0;
   }
@@ -787 +814 @@
   if (cSDptr) {
     cStatus = 0;
-    fits_delete_file(cSDptr, &cStatus);
+    if (fits_delete_file(cSDptr, &cStatus)) {
+      logMsg("ERROR: Failed to close and delete file.");
+    }
+
     cSDptr = 0;
   }
 }
+
+//------------------------------------------------------- SDFITSwriter::logMsg
+
+// Log a message.  If the current CFITSIO status value is non-zero, also log
+// the corresponding error message and dump the CFITSIO message stack.
+
+void SDFITSwriter::logMsg(const char *msg)
+{
+  PKSmsg::logMsg(msg);
+
+  if (cStatus) {
+    fits_get_errstatus(cStatus, cMsg);
+    PKSmsg::logMsg(cMsg);
+
+    while (fits_read_errmsg(cMsg)) {
+      PKSmsg::logMsg(cMsg);
+    }
+  }
+}

trunk/external/atnf/PKSIO/SDFITSwriter.h

@@ -2 +2 @@
 //# SDFITSwriter.h: ATNF CFITSIO interface class for SDFITS output.
 //#---------------------------------------------------------------------------
-//# Copyright (C) 2000-2007
+//# Copyright (C) 2000-2008
 //# Mark Calabretta, ATNF
 //#
@@ -27 +27 @@
 //#                        AUSTRALIA
 //#
-//# $Id: SDFITSwriter.h,v 19.7 2007/11/12 03:37:56 cal103 Exp $
+//# $Id: SDFITSwriter.h,v 19.9 2008-11-17 06:48:32 cal103 Exp $
 //#---------------------------------------------------------------------------
 //# Original: 2000/07/24, Mark Calabretta, ATNF
@@ -35 +35 @@
 #define ATNF_SDFITSWRITER_H
 
-#include <atnf/PKSIO/PKSMBrecord.h>
+#include <atnf/PKSIO/MBrecord.h>
+#include <atnf/PKSIO/PKSmsg.h>
 
 #include <fitsio.h>
+
+using namespace std;
 
 // <summary>
@@ -43 +46 @@
 // </summary>
 
-class SDFITSwriter
+class SDFITSwriter : public PKSmsg
 {
   public:
@@ -50 +53 @@
 
     // Destructor.
-    ~SDFITSwriter();
+    virtual ~SDFITSwriter();
 
     // Create a new SDFITSwriter and store static data.
@@ -71 +74 @@
 
     // Store time-variable data.
-    int write(PKSMBrecord &record);
+    int write(MBrecord &record);
 
     // Write a history record.
     int history(char* text);
-
-    // Print out CFITSIO error messages.
-    void reportError();
 
     // Close the SDFITS file.
@@ -90 +90 @@
          *cNChan, cNIF, *cNPol, cStatus;
     long cRow;
+
+    // Message handling.
+    char cMsg[256];
+    virtual void logMsg(const char *msg = 0x0);
 };
 

trunk/external/atnf/PKSIO/makefile

@@ -1 @@
-# $Id: makefile,v 19.0 2003/07/16 03:34:05 aips2adm Exp $
+# $Id: makefile,v 19.0 2003-07-16 03:34:05 aips2adm Exp $
 
 XLIBLIST := CFITSIO RPFITS