Index: branches/mergetest/external/atnf/PKSIO/SDFITSreader.cc
===================================================================
--- branches/mergetest/external/atnf/PKSIO/SDFITSreader.cc	(revision 1777)
+++ branches/mergetest/external/atnf/PKSIO/SDFITSreader.cc	(revision 1779)
@@ -38,8 +38,9 @@
 
 #include <atnf/pks/pks_maths.h>
-#include <atnf/PKSIO/PKSmsg.h>
 #include <atnf/PKSIO/MBrecord.h>
 #include <atnf/PKSIO/SDFITSreader.h>
 
+#include <casa/Logging/LogIO.h>
+#include <casa/Quanta/MVTime.h>
 #include <casa/math.h>
 #include <casa/stdio.h>
@@ -67,4 +68,7 @@
 const double D2R = PI / 180.0;
 
+// Class name
+const string className = "SDFITSreader" ;
+
 //---------------------------------------------------- SDFITSreader::(statics)
 
@@ -97,7 +101,5 @@
   cEndChan   = 0x0;
   cRefChan   = 0x0;
-
-  // By default, messages are written to stderr.
-  initMsg();
+  cPols      = 0x0;
 }
 
@@ -128,6 +130,5 @@
         int    &extraSysCal)
 {
-  // Clear the message stack.
-  clearMsg();
+  const string methodName = "open()" ;
 
   if (cSDptr) {
@@ -139,5 +140,5 @@
   if (fits_open_file(&cSDptr, sdName, READONLY, &cStatus)) {
     sprintf(cMsg, "ERROR: Failed to open SDFITS file\n       %s", sdName);
-    logMsg(cMsg);
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, cMsg);
     return 1;
   }
@@ -162,5 +163,5 @@
 
       } else {
-        logMsg("ERROR: Failed to locate SDFITS binary table.");
+        log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Failed to locate SDFITS binary table.");
         close();
         return 1;
@@ -201,5 +202,5 @@
     cNAxes = 5;
     if (readDim(DATA, 1, &cNAxes, cNAxis)) {
-      logMsg();
+      log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
       close();
       return 1;
@@ -220,10 +221,10 @@
     if (cNAxes < 4) {
       // Need at least four axes (for now).
-      logMsg("ERROR: DATA array contains fewer than four axes.");
+      log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "DATA array contains fewer than four axes.");
       close();
       return 1;
     } else if (cNAxes > 5) {
       // We support up to five axes.
-      logMsg("ERROR: DATA array contains more than five axes.");
+      log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "DATA array contains more than five axes.");
       close();
       return 1;
@@ -236,5 +237,5 @@
     findData(DATAXED, "DATAXED", TSTRING);
     if (cData[DATAXED].colnum < 0) {
-      logMsg("ERROR: DATA array column absent from binary table.");
+      log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "DATA array column absent from binary table.");
       close();
       return 1;
@@ -266,5 +267,5 @@
 
   if (cStatus) {
-    logMsg();
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
     close();
     return 1;
@@ -280,4 +281,5 @@
   char *timeCRPIX = 0;
   char *beamCRVAL = 0;
+  char *polCRVAL = 0;
 
   cFreqAxis   = -1;
@@ -297,4 +299,5 @@
     } else if (strncmp(ctype[iaxis], "STOKES", 6) == 0) {
       cStokesAxis = iaxis;
+      polCRVAL = CRVAL[iaxis];
 
     } else if (strncmp(ctype[iaxis], "RA", 2) == 0) {
@@ -318,5 +321,6 @@
         sprintf(cMsg, "DATA array contains a TIME axis of length %ld.",
           cNAxisTime);
-        logMsg(cMsg);
+        //logMsg(cMsg);
+        log(LogOrigin( className, methodName, WHERE ), LogIO::NORMAL, cMsg);
       }
 
@@ -337,7 +341,8 @@
   }
 
+
   // Check that required axes are present.
   if (cFreqAxis < 0 || cStokesAxis < 0 || cRaAxis < 0 || cDecAxis < 0) {
-    logMsg("ERROR: Could not find required DATA array axes.");
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Could not find required DATA array axes.");
     close();
     return 1;
@@ -403,6 +408,14 @@
   findData(WINDDIRE, "WINDDIRE", TFLOAT);       // Shared.
 
+  findData(STOKES,    polCRVAL,  TINT);
+  findData(SIG,       "SIG",     TSTRING);
+  findData(CAL,       "CAL",     TSTRING);
+
+  findData(RVSYS,     "RVSYS",   TDOUBLE);
+  findData(VFRAME,    "VFRAME",  TDOUBLE);
+  findData(VELDEF,    "VELDEF",  TSTRING);
+
   if (cStatus) {
-    logMsg();
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
     close();
     return 1;
@@ -433,4 +446,9 @@
   cCycleNo = 0;
   cLastUTC = 0.0;
+  for ( int i = 0 ; i < 4 ; i++ ) {
+    cGLastUTC[i] = 0.0 ;
+    cGLastScan[i] = -1 ;
+    cGCycleNo[i] = 0 ;
+  }
 
   // Beam number, 1-relative by default.
@@ -536,5 +554,5 @@
   fits_get_num_rows(cSDptr, &cNRow, &cStatus);
   if (!cNRow) {
-    logMsg("ERROR: Table contains no entries.");
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Table contains no entries.");
     close();
     return 1;
@@ -550,5 +568,5 @@
                       &beamNul, beamCol, &anynul, &cStatus)) {
       delete [] beamCol;
-      logMsg();
+      log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
       close();
       return 1;
@@ -565,5 +583,5 @@
       if (beamCol[irow] < cBeam_1rel) {
         delete [] beamCol;
-        logMsg("ERROR: SDFITS file contains invalid beam number.");
+        log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "SDFITS file contains invalid beam number.");
         close();
         return 1;
@@ -606,5 +624,5 @@
                       &IFNul, IFCol, &anynul, &cStatus)) {
       delete [] IFCol;
-      logMsg();
+      log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
       close();
       return 1;
@@ -621,5 +639,5 @@
       if (IFCol[irow] < cIF_1rel) {
         delete [] IFCol;
-        logMsg("ERROR: SDFITS file contains invalid IF number.");
+        log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "SDFITS file contains invalid IF number.");
         close();
         return 1;
@@ -653,5 +671,5 @@
             // Variable dimension array.
             if (readDim(DATA, irow+1, &cNAxes, cNAxis)) {
-              logMsg();
+              log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
               close();
               return 1;
@@ -681,5 +699,5 @@
 
           if (readDim(XPOLDATA, irow+1, &nAxis, nAxes)) {
-            logMsg();
+            log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE );
             close();
             return 1;
@@ -719,4 +737,44 @@
     cNPol[0] = cNIF;
     cNIF = 1;
+  }
+
+  // For GBT data that stores spectra for each polarization in separate rows
+  if ( cData[STOKES].colnum > 0 ) {
+    int *stokesCol = new int[cNRow];
+    int stokesNul = 1;
+    int   anynul;
+    if (fits_read_col(cSDptr, TINT, cData[STOKES].colnum, 1, 1, cNRow,
+                      &stokesNul, stokesCol, &anynul, &cStatus)) {
+      delete [] stokesCol;
+      log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
+      close();
+      return 1;
+    }
+
+    vector<int> pols ;
+    pols.push_back( stokesCol[0] ) ;
+    for ( int i = 0 ; i < cNRow ; i++ ) {
+      bool pmatch = false ;
+      for ( uint j = 0 ; j < pols.size() ; j++ ) {
+        if ( stokesCol[i] == pols[j] ) {
+          pmatch = true ;
+          break ;
+        }
+      }
+      if ( !pmatch ) {
+        pols.push_back( stokesCol[i] ) ;
+      }
+    }
+
+    cPols = new int[pols.size()] ;
+    for ( uint i = 0 ; i < pols.size() ; i++ ) {
+      cPols[i] = pols[i] ;
+    }
+
+    for ( int i = 0 ; i < cNIF ; i++ ) {
+      cNPol[i] = pols.size() ;
+    }
+
+    delete [] stokesCol ;
   }
 
@@ -806,4 +864,6 @@
         double &bandwidth)
 {
+  const string methodName = "getHeader()" ;
+  
   // Has the file been opened?
   if (!cSDptr) {
@@ -879,4 +939,7 @@
 
     // Look for VELFRAME, written by earlier versions of Livedata.
+    //
+    // Added few more codes currently (as of 2009 Oct) used in the GBT
+    // SDFITS (based io_sdfits_define.pro of GBTIDL). - TT
     if (readParm("VELFRAME", TSTRING, dopplerFrame)) {	// Additional.
       // No, try digging it out of the CTYPE card (AIPS convention).
@@ -890,7 +953,28 @@
           // LSR unqualified usually means LSR (kinematic).
           strcpy(dopplerFrame, "LSRK");
+        } else if (strcmp(dopplerFrame, "LSD") == 0) {
+          // LSR as a dynamical defintion 
+          strcpy(dopplerFrame, "LSRD");
         } else if (strcmp(dopplerFrame, "HEL") == 0) {
           // Almost certainly barycentric.
           strcpy(dopplerFrame, "BARYCENT");
+        } else if (strcmp(dopplerFrame, "BAR") == 0) {
+          // barycentric.
+          strcpy(dopplerFrame, "BARYCENT");
+        } else if (strcmp(dopplerFrame, "OBS") == 0) {
+          // observed or topocentric.
+          strcpy(dopplerFrame, "TOPO");
+        } else if (strcmp(dopplerFrame, "GEO") == 0) {
+          // geocentric 
+          strcpy(dopplerFrame, "GEO");
+        } else if (strcmp(dopplerFrame, "GAL") == 0) {
+          // galactic 
+          strcpy(dopplerFrame, "GAL");
+        } else if (strcmp(dopplerFrame, "LGR") == 0) {
+          // Local group 
+          strcpy(dopplerFrame, "LGROUP");
+        } else if (strcmp(dopplerFrame, "CMB") == 0) {
+          // Cosimic Microwave Backgroup
+          strcpy(dopplerFrame, "CMB");
         }
       } else {
@@ -898,5 +982,4 @@
       }
     }
-
     // Translate to FITS standard names.
     if (strncmp(dopplerFrame, "TOP", 3) == 0) {
@@ -908,9 +991,15 @@
     } else if (strncmp(dopplerFrame, "BARY", 4) == 0) {
       strcpy(dopplerFrame, "BARYCENT");
-    }
-  }
-
+    } else if (strncmp(dopplerFrame, "GAL", 3) == 0) {
+      strcpy(dopplerFrame, "GALACTOC");
+    } else if (strncmp(dopplerFrame, "LGROUP", 6) == 0) {
+      strcpy(dopplerFrame, "LOCALGRP");
+    } else if (strncmp(dopplerFrame, "CMB", 3) == 0) {
+      strcpy(dopplerFrame, "CMBDIPOL");
+    }
+  }
+  
   if (cStatus) {
-    logMsg();
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
     return 1;
   }
@@ -922,5 +1011,5 @@
 
   if (cStatus) {
-    logMsg();
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
     return 1;
   }
@@ -938,4 +1027,6 @@
         double* &endFreq)
 {
+  const string methodName = "getFreqInfo()" ;
+
   float  fqRefPix;
   double fqDelt, fqRefVal;
@@ -952,5 +1043,5 @@
                       &IFNul, IFCol, &anynul, &cStatus)) {
       delete [] IFCol;
-      logMsg();
+      log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
       close();
       return 1;
@@ -1015,4 +1106,6 @@
         double* &positions)
 {
+  const string methodName = "findRange()" ;
+
   // Has the file been opened?
   if (!cSDptr) {
@@ -1036,5 +1129,5 @@
       delete [] beamCol;
       delete [] sel;
-      logMsg();
+      log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
       return 1;
     }
@@ -1056,5 +1149,5 @@
       delete [] IFCol;
       delete [] sel;
-      logMsg();
+      log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
       return 1;
     }
@@ -1088,5 +1181,5 @@
     if (cData[AZIMUTH].colnum  < 0 ||
         cData[ELEVATIO].colnum < 0) {
-      logMsg("WARNING: Azimuth/elevation information absent.");
+      log(LogOrigin( className, methodName, WHERE ), LogIO::WARN, "Azimuth/elevation information absent.");
       cStatus = -1;
 
@@ -1115,5 +1208,5 @@
         cData[FOCUSROT].colnum < 0 ||
         cData[ELEVATIO].colnum < 0) {
-      logMsg("WARNING: ZPA/elevation information absent.");
+      log(LogOrigin( className, methodName, WHERE ), LogIO::WARN, "ZPA/elevation information absent.");
       cStatus = -1;
 
@@ -1191,6 +1284,7 @@
           cData[PARANGLE].colnum < 0 ||
           cData[FOCUSROT].colnum < 0) {
-        logMsg("WARNING: Insufficient information to compute feed-plane\n"
-               "         coordinates.");
+        log( LogOrigin( className, methodName, WHERE ), LogIO::WARN, 
+             "Insufficient information to compute feed-plane\n"
+             "         coordinates.");
         cStatus = -1;
 
@@ -1241,5 +1335,5 @@
     nSel = 0;
     delete [] positions;
-    logMsg();
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
     cStatus = 0;
     return 1;
@@ -1257,11 +1351,13 @@
         MBrecord &mbrec)
 {
+  const string methodName = "read()" ;
+
   // Has the file been opened?
   if (!cSDptr) {
     return 1;
   }
-
   // Find the next selected beam and IF.
   short iBeam = 0, iIF = 0;
+  int iPol = -1 ;
   while (1) {
     if (++cTimeIdx > cNAxisTime) {
@@ -1321,9 +1417,20 @@
               sprintf(cMsg, "ALFA cal factors for beam %d: %.3e, %.3e",
                 iBeam+1, sALFAcal[iBeam][0], sALFAcal[iBeam][1]);
-              logMsg(cMsg);
+              log(LogOrigin( className, methodName, WHERE ), LogIO::NORMAL, cMsg);
+              //logMsg(cMsg);
             }
           }
         }
 
+        // for GBT SDFITS
+        if (cData[STOKES].colnum > 0 ) {
+          readData(STOKES, cRow, &iPol ) ;
+          for ( int i = 0 ; i < cNPol[iIF] ; i++ ) {
+            if ( cPols[i] == iPol ) {
+              iPol = i ;
+              break ;
+            }
+          }
+        }
         break;
       }
@@ -1370,8 +1477,64 @@
   }
 
+  if ( iPol != -1 ) {
+    if ( mbrec.scanNo != cGLastScan[iPol] ) {
+      cGLastScan[iPol] = mbrec.scanNo ;
+      cGCycleNo[iPol] = 0 ;
+      mbrec.cycleNo = ++cGCycleNo[iPol] ;
+    }
+    else {
+      mbrec.cycleNo = ++cGCycleNo[iPol] ;
+    }
+  }
+
   readData(EXPOSURE, cRow, &mbrec.exposure);
 
   // Source identification.
   readData(OBJECT, cRow, mbrec.srcName);
+
+  if ( iPol != -1 ) {
+    char obsmode[32] ;
+    readData( OBSMODE, cRow, obsmode ) ;
+    char sig[1] ;
+    char cal[1] ;
+    readData( SIG, cRow, sig ) ;
+    readData( CAL, cRow, cal ) ;
+    if ( strstr( obsmode, "PSWITCH" ) != NULL ) {
+      // position switch
+      strcat( mbrec.srcName, "_p" ) ;
+      if ( strstr( obsmode, "PSWITCHON" ) != NULL ) {
+        strcat( mbrec.srcName, "s" ) ;
+      }
+      else if ( strstr( obsmode, "PSWITCHOFF" ) != NULL ) {
+        strcat( mbrec.srcName, "r" ) ;
+      }
+    }
+    else if ( strstr( obsmode, "Nod" ) != NULL ) {
+      // nod
+      strcat( mbrec.srcName, "_n" ) ;
+      if ( sig[0] == 'T' ) {
+        strcat( mbrec.srcName, "s" ) ;
+      }
+      else {
+        strcat( mbrec.srcName, "r" ) ;
+      }
+    }
+    else if ( strstr( obsmode, "FSWITCH" ) != NULL ) {
+      // frequency switch
+      strcat( mbrec.srcName, "_f" ) ;
+      if ( sig[0] == 'T' ) {
+        strcat( mbrec.srcName, "s" ) ;
+      }
+      else {
+        strcat( mbrec.srcName, "r" ) ;
+      }
+    }
+    if ( cal[0] == 'T' ) {
+      strcat( mbrec.srcName, "c" ) ;
+    }
+    else {
+      strcat( mbrec.srcName, "o" ) ;
+    }
+  }
 
   readData(OBJ_RA,  cRow, &mbrec.srcRA);
@@ -1421,4 +1584,7 @@
   int nPol = cNPol[iIF];
 
+  if ( cData[STOKES].colnum > 0 )
+    nPol = 1 ;
+
   if (cGetSpectra || cGetXPol) {
     int nxpol = cGetXPol ? 2*nChan : 0;
@@ -1430,4 +1596,5 @@
   mbrec.nChan[0] = nChan;
   mbrec.nPol[0]  = nPol;
+  mbrec.polNo = iPol ;
 
   readData(FqRefPix, cRow, mbrec.fqRefPix);
@@ -1448,5 +1615,5 @@
 
   if (cStatus) {
-    logMsg();
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
     return 1;
   }
@@ -1485,5 +1652,5 @@
 
   if (cStatus) {
-    logMsg();
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
     return 1;
   }
@@ -1535,8 +1702,8 @@
 
         if ((status = readDim(DATA, cRow, &naxes, cNAxis))) {
-          logMsg();
+          log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
 
         } else if ((status = (naxes != cNAxes))) {
-          logMsg("ERROR: DATA array dimensions changed.");
+          log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "DATA array dimensions changed.");
         }
 
@@ -1552,5 +1719,5 @@
           blc, trc, inc, 0, mbrec.spectra[0] + iPol*nChan, &anynul,
           &cStatus)) {
-        logMsg();
+        log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
         delete [] blc;
         delete [] trc;
@@ -1613,5 +1780,5 @@
             cNAxis, blc, trc, inc, 0, mbrec.flagged[0] + iPol*nChan, &anynul,
             &cStatus)) {
-          logMsg();
+          log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
           delete [] blc;
           delete [] trc;
@@ -1664,5 +1831,5 @@
     if (fits_read_subset_flt(cSDptr, cData[XPOLDATA].colnum, nAxis, nAxes,
         blc, trc, inc, 0, mbrec.xpol[0], &anynul, &cStatus)) {
-      logMsg();
+      log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
       delete [] blc;
       delete [] trc;
@@ -1695,5 +1862,5 @@
 
   if (cStatus) {
-    logMsg();
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
     return 1;
   }
@@ -1729,6 +1896,16 @@
   mbrec.windAz    *= D2R;
 
+  // For GBT data, source velocity can be evaluated
+  if ( cData[RVSYS].colnum > 0 && cData[VFRAME].colnum > 0 ) {
+    float vframe;
+    readData(VFRAME, cRow, &vframe);
+    float rvsys;
+    readData(RVSYS,  cRow, &rvsys);
+    //mbrec.srcVelocity = rvsys - vframe ;
+    mbrec.srcVelocity = rvsys ;
+  }
+
   if (cStatus) {
-    logMsg();
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
     return 1;
   }
@@ -1756,21 +1933,24 @@
 }
 
-//------------------------------------------------------- SDFITSreader::logMsg
+//------------------------------------------------------- SDFITSreader::log
 
 // 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)
+void SDFITSreader::log(LogOrigin origin, LogIO::Command cmd, const char *msg)
 {
-  FITSreader::logMsg(msg);
+  LogIO os( origin ) ;
+
+  os << msg << endl ;
 
   if (cStatus > 0) {
     fits_get_errstatus(cStatus, cMsg);
-    FITSreader::logMsg(cMsg);
+    os << cMsg << endl ;
 
     while (fits_read_errmsg(cMsg)) {
-      FITSreader::logMsg(cMsg);
-    }
-  }
+      os << cMsg << endl ;
+    }
+  }
+  os << LogIO::POST ;
 }
 
@@ -2107,4 +2287,12 @@
   if (cData[TIME].colnum >= 0) {
     readData(TIME, iRow, &utc);
+  } else if (cGBT) {
+    Int yy, mm ;
+    Double dd, hour, min, sec ;
+    sscanf( datobs, "%d-%d-%lfT%lf:%lf:%lf", &yy, &mm, &dd, &hour, &min, &sec ) ;
+    dd = dd + ( hour * 3600.0 + min * 60.0 + sec ) / 86400.0 ;
+    MVTime mvt( yy, mm, dd ) ;
+    dd = mvt.day() ;
+    utc = fmod( dd, 1.0 ) * 86400.0 ;
   } else if (cNAxisTime > 1) {
     double timeDelt, timeRefPix, timeRefVal;
@@ -2153,4 +2341,6 @@
         short iPol)
 {
+  const string methodName = "alfaCal()" ;
+
   int  calOn;
   char chars[32];
@@ -2205,5 +2395,5 @@
   if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAxes, cNAxis,
       blc, trc, inc, 0, spectrum, &anynul, &cStatus)) {
-    logMsg();
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
     delete [] blc;
     delete [] trc;
