Changeset 1779 for branches/mergetest/external/atnf/PKSIO

branches/mergetest

Property svn:mergeinfo changed

merged: 1386-1389,1393,1400-1401,1445-1448,1451,1453-1461,1480,1482,1489-1496,1498,1514-1522,1530-1531,1533,1562,1566,1601-1616,1619-1634,1637,1639-1640,1649-1662,1666-1668,1671-1686,1693,1698-1707,1714,1717,1719,1723-1724,1729,1732,1757-1768,1770-1771

branches/mergetest/external/atnf/PKSIO/FITSreader.cc

-              r1720
+              r1779
         const int getSpectra,
         const int getXPol,
+        const int getFeedPos,
+        const int getPointing,
         const int coordSys)
+{
 …
   cGetSpectra = getSpectra && cHaveSpectra;
   cGetXPol    = getXPol    && cGetXPol;
+  cGetFeedPos = getFeedPos;
   cCoordSys   = coordSys;
   return maxNChan;

branches/mergetest/external/atnf/PKSIO/FITSreader.h

-              r1720
+              r1779
 #include <atnf/PKSIO/MBrecord.h>
-#include <atnf/PKSIO/PKSmsg.h>
 using namespace std;
 // <summary>
 …
 // </summary>
+class FITSreader : public PKSmsg
+//class FITSreader
+class FITSreader
+{
   public:
 …
         const int refChan[],
         const int getSpectra = 1,
+        const int getXPol  = 0,
+        const int getXPol = 0,
+        const int getFeedPos = 0,
+        const int getPointing = 0,
         const int coordSys = 0);
     // Find the range in time and position of the data selected.
 …
     // Read the next data record.
     virtual int read(
+//        PKSMBrecord &record) = 0;
         MBrecord &record) = 0;
 …
   protected:
     int  *cBeams, *cEndChan, cCoordSys, cGetSpectra, cGetXPol, cHaveBase,
          cHaveSpectra, *cHaveXPol, *cIFs, cNBeam, *cNChan, cNIF, *cNPol,
          *cRefChan, *cStartChan;
+    int  *cBeams, *cEndChan, cGetFeedPos, cCoordSys, cGetSpectra, cGetXPol,
+           cHaveBase, cHaveSpectra, *cHaveXPol, *cIFs, cNBeam, *cNChan, cNIF,
+           *cNPol, *cRefChan, *cStartChan;
     // For use in constructing messages.
     char cMsg[256];
 };

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

-              r1720
+              r1779
 #include <atnf/PKSIO/MBrecord.h>
+#include <casa/Logging/LogIO.h>
 #include <casa/math.h>
 #include <casa/iostream.h>
 …
 const double HALFPI = PI / 2.0;
 const double R2D = 180.0 / PI;
+// Class name
+const string className = "MBFITSreader" ;
 //------------------------------------------------- MBFITSreader::MBFITSreader
 …
   // Tell RPFITSIN not to report errors directly.
+  iostat_.errlun = -1;
+  // By default, messages are written to stderr.
+  initMsg();
+  //iostat_.errlun = -1;
+}
 …
         int  &extraSysCal)
+{
   // Clear the message stack.
   clearMsg();
+  const string methodName = "open()" ;
+  LogIO os( LogOrigin( className, methodName, WHERE ) ) ;
   if (cMBopen) {
 …
   int jstat = -3;
   if (rpfitsin(jstat)) {
     sprintf(cMsg, "ERROR: Failed to open MBFITS file\n       %s", rpname);
     logMsg(cMsg);
+    sprintf(cMsg, "Failed to open MBFITS file\n%s", rpname);
+    os << LogIO::SEVERE << cMsg << LogIO::POST ;
     return 1;
+  }
 …
   jstat = -1;
   if (rpfitsin(jstat)) {
     sprintf(cMsg, "ERROR: Failed to read MBFITS header in file\n"
                   "       %s", rpname);
     logMsg(cMsg);
+    sprintf(cMsg, "Failed to read MBFITS header in file\n"
+                  "%s", rpname);
+    os << LogIO::SEVERE << cMsg << LogIO::POST ;
     close();
     return 1;
 …
   // Non-ATNF data may not store the position in (u,v,w).
   if (strncmp(names_.sta, "tid", 3) == 0) {
     sprintf(cMsg, "WARNING: Found Tidbinbilla data");
+    sprintf(cMsg, "Found Tidbinbilla data");
     cSUpos = 1;
   } else if (strncmp(names_.sta, "HOB", 3) == 0) {
     sprintf(cMsg, "WARNING: Found Hobart data");
+    sprintf(cMsg, "Found Hobart data");
     cSUpos = 1;
   } else if (strncmp(names_.sta, "CED", 3) == 0) {
     sprintf(cMsg, "WARNING: Found Ceduna data");
+    sprintf(cMsg, "Found Ceduna data");
     cSUpos = 1;
   } else {
 …
   if (cSUpos) {
     strcat(cMsg, ", using telescope position\n         from SU table.");
     logMsg(cMsg);
+    os << LogIO::WARN << cMsg << LogIO::POST ;
     cInterp = 0;
+  }
 …
   if (cNBeam <= 0) {
     logMsg("ERROR: Couldn't determine number of beams.");
+    os << LogIO::SEVERE << "Couldn't determine number of beams." << LogIO::POST ;
     close();
     return 1;
 …
       sprintf(cMsg,
         "WARNING: RPFITSIN returned beam number %2d for AN table\n"
         "         entry %2d with name '%.8s'", beamNo, iBeam+1, sta);
+        "RPFITSIN returned beam number %2d for AN table\n"
+        "entry %2d with name '%.8s'", beamNo, iBeam+1, sta);
       char text[8];
 …
+      }
       logMsg(cMsg);
+      os << LogIO::WARN << cMsg << LogIO::POST ;
+    }
 …
   // Read the first syscal record.
   if (rpget(1, cEOS)) {
     logMsg("ERROR: Failed to read first syscal record.");
+    os << LogIO::SEVERE << "Failed to read first syscal record." << LogIO::POST ;
     close();
     return 1;
 …
         double &bandwidth)
+{
+  const string methodName = "getHeader()" ;
+  LogIO os( LogOrigin( className, methodName, WHERE ) ) ;
   if (!cMBopen) {
     logMsg("ERROR: An MBFITS file has not been opened.");
+    os << LogIO::SEVERE << "An MBFITS file has not been opened." << LogIO::POST ;
     return 1;
+  }
 …
         MBrecord &MBrec)
+{
+  const string methodName = "read()" ;
+  LogIO os( LogOrigin( className, methodName, WHERE ) ) ;
   int beamNo = -1;
   int haveData, pCode = 0, status;
 …
   if (!cMBopen) {
     logMsg("ERROR: An MBFITS file has not been opened.");
+    os << LogIO::SEVERE << "An MBFITS file has not been opened." << LogIO::POST ;
     return 1;
+  }
 …
 #ifdef PKSIO_DEBUG
         fprintf(stderr, "\nEnd-of-file detected, flushing last cycle.\n");
+        os << LogIO::DEBUGGING << "\nEnd-of-file detected, flushing last cycle.\n" << LogIO::POST ;
 #endif
 …
           if (cNBin > 1 && cNBeamSel > 1) {
+            logMsg("ERROR: Cannot handle binning mode for multiple beams.\n"
+                   "       Select a single beam for input.");
+            os << LogIO::SEVERE << "Cannot handle binning mode for multiple beams.\nSelect a single beam for input." << LogIO::POST ;
             close();
             return 1;
 …
           // the start of the next.
 #ifdef PKSIO_DEBUG
+          fprintf(stderr, "Change-of-day on cUTC: %.1f -> %.1f\n",
+            cPrevUTC, cUTC);
+          char buf[256] ;
+          sprintf(buf, "Change-of-day on cUTC: %.1f -> %.1f\n", cPrevUTC, cUTC);
+          os << LogIO::DEBUGGING << buf << LogIO::POST ;
 #endif
           // Can't change the recorded value of cUTC directly (without also
 …
           // an offset to be applied when comparing integration timestamps.
           cod = 86400.0;
+        }
+        }
         if ((cUTC+cod) < cPrevUTC - 1.0) {
 …
             // All other data should be fully time ordered.
             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,
+              "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);
+            //logMsg(cMsg);
+            os << LogIO::WARN << cMsg << LogIO::POST ;
             cUTC += 86400.0;
+          }
 …
+        }
         fprintf(stderr, "\n In:%4d%4d%3d%3d  %.3f %c %.3f (%+.3fs) - "
+        sprintf(buf, "\n In:%4d%4d%3d%3d  %.3f %c %.3f (%+.3fs) - "
           "%sflushing\n", cScanNo, cCycleNo, beamNo, cIFno, cUTC, rel, cW, dt,
           cFlushing ? "" : "not ");
+        os << LogIO::DEBUGGING << buf << LogIO::POST ;
         if (cEOS) {
+          fprintf(stderr, "Start of new scan, flushing previous scan.\n");
+          sprintf(buf, "Start of new scan, flushing previous scan.\n");
+          os << LogIO::DEBUGGING << buf << LogIO::POST ;
+        }
 #endif
 …
 #ifdef PKSIO_DEBUG
       fprintf(stderr, "This (%d) ra, dec, UTC: %9.4f %9.4f %10.3f %9.4f\n",
+      sprintf(buf, "This (%d) ra, dec, UTC: %9.4f %9.4f %10.3f %9.4f\n",
         iMBuff->cycleNo, thisRA*R2D, thisDec*R2D, thisUTC, thisPA*R2D);
+      os << LogIO::DEBUGGING << buf << LogIO::POST ;
 #endif
 …
 #ifdef PKSIO_DEBUG
           fprintf(stderr, "Next (%d) ra, dec, UTC: %9.4f %9.4f %10.3f "
+          sprintf(buf, "Next (%d) ra, dec, UTC: %9.4f %9.4f %10.3f "
             "(0.000s)\n", cCycleNo, cU*R2D, cV*R2D, cW);
+          os << LogIO::DEBUGGING << buf << LogIO::POST ;
 #endif
 …
 #ifdef PKSIO_DEBUG
           fprintf(stderr, "Next (%d) ra, dec, UTC: %9.4f %9.4f %10.3f "
+          sprintf(buf, "Next (%d) ra, dec, UTC: %9.4f %9.4f %10.3f "
             "(%+.3fs)\n", cCycleNo, nextRA*R2D, nextDec*R2D, nextUTC,
             utcDiff(nextUTC, thisUTC));
+          os << LogIO::DEBUGGING << buf << LogIO::POST ;
 #endif
 …
 #ifdef PKSIO_DEBUG
       double avRate = sqrt(cAvRate[0]*cAvRate[0] + cAvRate[1]*cAvRate[1]);
       fprintf(stderr, "RA, Dec, Av & PA rates: %8.4f %8.4f %8.4f %8.4f "
+      sprintf(buf, "RA, Dec, Av & PA rates: %8.4f %8.4f %8.4f %8.4f "
         "pCode %d\n", raRate*R2D, decRate*R2D, avRate*R2D, paRate*R2D, pCode);
+      os << LogIO::DEBUGGING << buf << LogIO::POST ;
 #endif
 …
 #ifdef PKSIO_DEBUG
         fprintf(stderr, "Intp (%d) ra, dec, UTC: %9.4f %9.4f %10.3f (pCode, "
+        sprintf(buf, "Intp (%d) ra, dec, UTC: %9.4f %9.4f %10.3f (pCode, "
           "age: %d %.1fs)\n", iMBuff->cycleNo, cBuffer[jbuff].ra*R2D,
           cBuffer[jbuff].dec*R2D, cBuffer[jbuff].utc, iMBuff->pCode,
           iMBuff->rateAge);
+        os << LogIO::DEBUGGING << buf << LogIO::POST ;
 #endif
+      }
 …
 #ifdef PKSIO_DEBUG
       fprintf(stderr, "Out:%4d%4d%3d%3d\n", MBrec.scanNo, MBrec.cycleNo,
+      sprintf(buf, "Out:%4d%4d%3d%3d\n", MBrec.scanNo, MBrec.cycleNo,
         MBrec.beamNo, MBrec.IFno[0]);
+      os << LogIO::DEBUGGING << buf << LogIO::POST ;
 #endif
 …
         // Sanity check on the number of IFs in the new scan.
         if (if_.n_if != cNIF) {
           sprintf(cMsg, "WARNING: Scan %d has %d IFs instead of %d, "
+          sprintf(cMsg, "Scan %d has %d IFs instead of %d, "
             "continuing.", cScanNo, if_.n_if, cNIF);
           logMsg(cMsg);
+          os << LogIO::WARN << cMsg << LogIO::POST ;
+        }
+      }
 …
 #ifdef PKSIO_DEBUG
+      fprintf(stderr, "Buf:%4d%4d%3d%3d\n", cScanNo, cCycleNo, beamNo, cIFno);
+      sprintf(buf, "Buf:%4d%4d%3d%3d\n", cScanNo, cCycleNo, beamNo, cIFno);
+      os << LogIO::DEBUGGING << buf << LogIO::POST ;
 #endif
 …
         // 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,
+        sprintf(cMsg, "Integration cycle %d:%d, beam %2d, \n"
+                      "IF %d was duplicated.", cScanNo, cCycleNo-1,
                       beamNo, cIFno);
         logMsg(cMsg);
+        os << LogIO::WARN << cMsg << LogIO::POST ;
+      }
       iMBuff->nChan[iIFSel] = nChan;
 …
 int MBFITSreader::rpget(int syscalonly, int &EOS)
+{
+  const string methodName = "rpget()" ;
+  LogIO os( LogOrigin( className, methodName, WHERE ) ) ;
   EOS = 0;
 …
       // Read failed; retry.
       numErr++;
       logMsg("WARNING: RPFITS read failed - retrying.");
+      os << LogIO::WARN << "RPFITS read failed - retrying." << LogIO::POST ;
       jstat = 0;
       break;
 …
     default:
       // Shouldn't reach here.
       sprintf(cMsg, "WARNING: Unrecognized RPFITSIN return code: %d "
+      sprintf(cMsg, "Unrecognized RPFITSIN return code: %d "
                     "(retrying).", jstat);
       logMsg(cMsg);
+      os << LogIO::WARN << cMsg << LogIO::POST ;
       jstat = 0;
       break;
 …
+  }
   logMsg("ERROR: RPFITS read failed too many times.");
+  os << LogIO::SEVERE << "RPFITS read failed too many times." << LogIO::POST ;
   return 2;
+}
 …
   // Handle messages from RPFITSIN.
+/**
   if (names_.errmsg[0] != ' ') {
     int i;
 …
     logMsg(cMsg);
+  }
+**/
   return jstat;
+}

branches/mergetest/external/atnf/PKSIO/MBrecord.cc

-              r1720
+              r1779
   refBeam = other.refBeam;
+  polNo = other.polNo ;
+  srcVelocity = other.srcVelocity ;
   return *this;
+}
 …
   refBeam = other.refBeam;
+  polNo = other.polNo ;
+  srcVelocity = other.srcVelocity ;
   return 0;
+}

branches/mergetest/external/atnf/PKSIO/MBrecord.h

r1720	r1779
166	166	short refBeam; // Reference beam, in beam-switching (MX)
167	167	// mode (added 1999/03/17).
	168	int polNo ; // polarization ID
	169	float srcVelocity ; // source velocity w.r.t. reference frame
168	170
169	171	private:

branches/mergetest/external/atnf/PKSIO/PKSFITSreader.cc

-              r1720
+              r1779
 //#---------------------------------------------------------------------------
-#include <atnf/PKSIO/PKSmsg.h>
 #include <atnf/PKSIO/MBFITSreader.h>
 #include <atnf/PKSIO/SDFITSreader.h>
 …
 #include <casa/BasicMath/Math.h>
 #include <casa/Quanta/MVTime.h>
+#include <casa/Logging/LogIO.h>
 //----------------------------------------------- PKSFITSreader::PKSFITSreader
 …
     cReader = new MBFITSreader(retry, interpolate ? 1 : 0);
+  }
-  // By default, messages are written to stderr.
-  initMsg();
+}
 …
+}
-//------------------------------------------------------ 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;
+}
 //-------------------------------------------------------- PKSFITSreader::open
 …
 Int PKSFITSreader::open(
         const String fitsName,
+        const String antenna,
         Vector<Bool> &beams,
         Vector<Bool> &IFs,
 …
         Bool   &haveSpectra)
+{
-  clearMsg();
   int    extraSysCal, haveBase_, *haveXPol_, haveSpectra_, nBeam, *nChan_,
          nIF, *nPol_, status;
 …
                          nChan_, nPol_, haveXPol_, haveBase_, haveSpectra_,
                          extraSysCal);
   logMsg(cReader->getMsg());
   cReader->clearMsg();
+  //logMsg(cReader->getMsg());
+  //cReader->clearMsg();
   if (status) {
     return status;
 …
                               obsType_, bunit_, equinox_, radecsys,
                               dopplerFrame_, datobs, utc, refFreq, bandwidth);
   logMsg(cReader->getMsg());
   cReader->clearMsg();
+  //logMsg(cReader->getMsg());
+  //cReader->clearMsg();
   if (status) {
     return 1;
 …
   Int status = cReader->getFreqInfo(nIF, startfreq, endfreq);
   logMsg(cReader->getMsg());
   cReader->clearMsg();
+  //logMsg(cReader->getMsg());
+  //cReader->clearMsg();
   if (!status) {
     startFreq.takeStorage(IPosition(1,nIF), startfreq, TAKE_OVER);
 …
         const Bool getSpectra,
         const Bool getXPol,
+        const Bool getFeedPos,
+        const Bool getPointing,
         const Int  coordSys)
+{
 …
   cGetSpectra = getSpectra;
   cGetXPol    = getXPol;
+  cGetFeedPos = getFeedPos;
+  cGetPointing = getPointing;
   cCoordSys   = coordSys;
   uInt maxNChan = cReader->select(start, end, ref, cGetSpectra, cGetXPol,
                                   cCoordSys);
   logMsg(cReader->getMsg());
   cReader->clearMsg();
+                                  cGetFeedPos, cGetPointing, cCoordSys);
+  //logMsg(cReader->getMsg());
+  //cReader->clearMsg();
   delete [] end;
 …
   Int status = cReader->findRange(nRow, nSel, dateSpan, utcSpan, posns);
   logMsg(cReader->getMsg());
   cReader->clearMsg();
+  //logMsg(cReader->getMsg());
+  //cReader->clearMsg();
   if (!status) {
 …
+{
   Int status = cReader->read(cMBrec);
   logMsg(cReader->getMsg());
   cReader->clearMsg();
+  //logMsg(cReader->getMsg());
+  //cReader->clearMsg();
   if (status) {
 …
   pksrec.scanNo  = cMBrec.scanNo;
   pksrec.cycleNo = cMBrec.cycleNo;
+  pksrec.polNo = cMBrec.polNo ;
   // Extract MJD.
   Int day, month, year;
+  sscanf(cMBrec.datobs, "%4d-%2d-%2d", &year, &month, &day);
+  pksrec.mjd = MVTime(year, month, Double(day)).day() + cMBrec.utc/86400.0;
+  if ( strstr( cMBrec.datobs, "T" ) == NULL ) {
+    sscanf(cMBrec.datobs, "%4d-%2d-%2d", &year, &month, &day);
+    pksrec.mjd = MVTime(year, month, Double(day)).day() + cMBrec.utc/86400.0;
+  }
+  else {
+    Double dd, hour, min, sec ;
+    sscanf( cMBrec.datobs, "%4d-%2d-%2lfT%lf:%lf:%lf", &year, &month, &dd, &hour, &min, &sec ) ;
+    dd = dd + ( hour * 3600.0 + min * 60.0 + sec ) / 86400.0 ;
+    pksrec.mjd = MVTime(year, month, dd).day() ;
+  }
   pksrec.interval  = cMBrec.exposure;
 …
   pksrec.fieldName = trim(cMBrec.srcName);
   pksrec.srcName   = pksrec.fieldName;
+  int namelen = pksrec.srcName.length() ;
+  if ( namelen > 4 ) {
+    String srcsub = pksrec.srcName.substr( namelen-4, 4 ) ;
+    if ( srcsub.find( "_psc" ) != string::npos ) {
+      pksrec.fieldName = pksrec.srcName.substr( 0, namelen-4 ) ;
+      pksrec.srcName = pksrec.fieldName + "_ps_calon" ;
+    }
+    else if ( srcsub.find( "_pso" ) != string::npos ) {
+      pksrec.fieldName = pksrec.srcName.substr( 0, namelen-4 ) ;
+      pksrec.srcName = pksrec.fieldName + "_ps" ;
+    }
+    else if ( srcsub.find( "_prc" ) != string::npos ) {
+      pksrec.fieldName = pksrec.srcName.substr( 0, namelen-4 ) ;
+      pksrec.srcName = pksrec.fieldName + "_psr_calon" ;
+    }
+    else if ( srcsub.find( "_pro" ) != string::npos ) {
+      pksrec.fieldName = pksrec.srcName.substr( 0, namelen-4 ) ;
+      pksrec.srcName = pksrec.fieldName + "_psr" ;
+    }
+    else if ( srcsub.find( "_fsc" ) != string::npos ) {
+      pksrec.fieldName = pksrec.srcName.substr( 0, namelen-4 ) ;
+      pksrec.srcName = pksrec.fieldName + "_fs_calon" ;
+    }
+    else if ( srcsub.find( "_fso" ) != string::npos ) {
+      pksrec.fieldName = pksrec.srcName.substr( 0, namelen-4 ) ;
+      pksrec.srcName = pksrec.fieldName + "_fs" ;
+    }
+    else if ( srcsub.find( "_frc" ) != string::npos ) {
+      pksrec.fieldName = pksrec.srcName.substr( 0, namelen-4 ) ;
+      pksrec.srcName = pksrec.fieldName + "_fsr_calon" ;
+    }
+    else if ( srcsub.find( "_fro" ) != string::npos ) {
+      pksrec.fieldName = pksrec.srcName.substr( 0, namelen-4 ) ;
+      pksrec.srcName = pksrec.fieldName + "_fsr" ;
+    }
+    else if ( srcsub.find( "_nsc" ) != string::npos || srcsub.find( "_nrc" ) != string::npos ) {
+      pksrec.fieldName = pksrec.srcName.substr( 0, namelen-4 ) ;
+      pksrec.srcName = pksrec.fieldName + "_nod_calon" ;
+    }
+    else if ( srcsub.find( "_nso" ) != string::npos || srcsub.find( "_nro" ) != string::npos ) {
+      pksrec.fieldName = pksrec.srcName.substr( 0, namelen-4 ) ;
+      pksrec.srcName = pksrec.fieldName + "_nod" ;
+    }
+  }
   pksrec.srcDir.resize(2);
 …
   pksrec.srcPM(0)  = 0.0;
   pksrec.srcPM(1)  = 0.0;
   pksrec.srcVel    = 0.0;
+  pksrec.srcVel    = cMBrec.srcVelocity;
   pksrec.obsType   = trim(cMBrec.obsType);
 …
   pksrec.bandwidth = chanWidth * nChan;
   pksrec.freqInc   = cMBrec.fqDelt[0];
+  pksrec.restFreq  = cMBrec.restFreq;
+  pksrec.restFreq.resize(1) ;
+  pksrec.restFreq(0)  = cMBrec.restFreq;
   pksrec.tcal.resize(nPol);
 …
+{
   cReader->close();
   logMsg(cReader->getMsg());
   cReader->clearMsg();
+  //logMsg(cReader->getMsg());
+  //cReader->clearMsg();
+}

branches/mergetest/external/atnf/PKSIO/PKSFITSreader.h

-              r1720
+              r1779
     virtual ~PKSFITSreader();
-    // Set message disposition.
-    virtual Int setMsg(
-        FILE *fd = 0x0);
     // Open the FITS file for reading.
     virtual Int open(
         const String fitsName,
+        const String antenna,
         Vector<Bool> &beams,
         Vector<Bool> &IFs,
 …
         const Bool getSpectra = True,
         const Bool getXPol    = False,
+        const Bool getFeedPos = False,
+        const Bool getPointing = False,
         const Int  coordSys   = 0);

branches/mergetest/external/atnf/PKSIO/PKSMS2reader.cc

-              r1720
+              r1779
 //#---------------------------------------------------------------------------
+#include <atnf/pks/pks_maths.h>
+#include <atnf/PKSIO/PKSmsg.h>
+#include <atnf/PKSIO/PKSrecord.h>
+#include <atnf/PKSIO/PKSMS2reader.h>
+// AIPS++ includes.
 #include <casa/stdio.h>
 #include <casa/Arrays/ArrayMath.h>
 …
 #include <ms/MeasurementSets/MSColumns.h>
 #include <tables/Tables.h>
+#include <casa/Quanta/MVTime.h>
+#include <casa/Quanta/MVAngle.h>
+#include <casa/BasicMath/Math.h>
+#include <casa/Logging/LogIO.h>
+#include <casa/Utilities/Sort.h>
+#include <measures/Measures/MeasConvert.h>
+#include <measures/Measures/MEpoch.h>
+#include <measures/Measures/MeasRef.h>
+// Parkes includes.
+#include <atnf/pks/pks_maths.h>
+#include <atnf/PKSIO/PKSrecord.h>
+#include <atnf/PKSIO/PKSMS2reader.h>
 //------------------------------------------------- PKSMS2reader::PKSMS2reader
 …
+{
   cMSopen = False;
-  // By default, messages are written to stderr.
-  initMsg();
+}
 …
 Int PKSMS2reader::open(
         const String msName,
+        const String antenna,
         Vector<Bool> &beams,
         Vector<Bool> &IFs,
 …
   cPKSMS  = MeasurementSet(msName);
+  // data selection by antenna
+  if ( antenna.length() == 0 ) {
+    cAntId.resize( 1 ) ;
+    cAntId[0] = 0 ;
+  }
+  else {
+    setupAntennaList( antenna ) ;
+    if ( cAntId.size() > 1 ) {
+      LogIO os( LogOrigin( "PKSMS2reader", "open()", WHERE ) ) ;
+      os << LogIO::WARN << "PKSMS2reader is not ready for multiple antenna selection. Use first antenna id " << cAntId[0] << "."<< LogIO::POST ;
+      Int tmp = cAntId[0] ;
+      cAntId.resize( 1 ) ;
+      cAntId[0] = tmp ;
+    }
+    stringstream ss ;
+    ss << "SELECT FROM $1 WHERE ANTENNA1 == ANTENNA2 && ANTENNA1 IN [" ;
+    for ( uInt i = 0 ; i < cAntId.size() ; i++ ) {
+      ss << cAntId[i] ;
+      if ( i == cAntId.size()-1 ) {
+        ss << "]" ;
+      }
+      else {
+        ss << "," ;
+      }
+    }
+    string taql = ss.str() ;
+    //cerr << "taql = " << taql << endl ;
+    cPKSMS = MeasurementSet( tableCommand( taql, cPKSMS ) ) ;
+  }
+  // taql access to the syscal table
+  cHaveSysCal = False;
+  if (cHaveSysCal=Table::isReadable(cPKSMS.sysCalTableName())) {
+    cSysCalTab = Table(cPKSMS.sysCalTableName());
+  }
+  // Lock the table for read access.
+  cPKSMS.lock(False);
   cIdx    = 0;
+  lastmjd = 0.0;
   cNRow   = cPKSMS.nrow();
   cMSopen = True;
-  // Lock the table for read access.
-  cPKSMS.lock(False);
   // Main MS table and subtable column access.
 …
   ROMSSysCalColumns       sysCalCols(cPKSMS.sysCal());
   ROMSWeatherColumns      weatherCols(cPKSMS.weather());
+  ROMSAntennaColumns      antennaCols(cPKSMS.antenna());
   // Column accessors for required columns.
 …
   cFieldIdCol.reference(msCols.fieldId());
   cFieldNameCol.reference(fieldCols.name());
+  cFieldDelayDirCol.reference(fieldCols.delayDir());
   cSrcIdCol.reference(fieldCols.sourceId());
+  cSrcId2Col.reference(sourceCols.sourceId());
   cSrcNameCol.reference(sourceCols.name());
   cSrcDirCol.reference(sourceCols.direction());
 …
   cStateIdCol.reference(msCols.stateId());
   cObsModeCol.reference(stateCols.obsMode());
+  cCalCol.reference(stateCols.cal());
+  cSigStateCol.reference(stateCols.sig());
+  cRefStateCol.reference(stateCols.ref());
   cDataDescIdCol.reference(msCols.dataDescId());
+  cSpWinIdCol.reference(dataDescCols.spectralWindowId());
   cChanFreqCol.reference(spWinCols.chanFreq());
+  cTotBWCol.reference(spWinCols.totalBandwidth());
   cWeatherTimeCol.reference(weatherCols.time());
 …
   cBeamNoCol.reference(msCols.feed1());
   cPointingCol.reference(pointingCols.direction());
+  cPointingTimeCol.reference(pointingCols.time());
   cSigmaCol.reference(msCols.sigma());
   cNumReceptorCol.reference(feedCols.numReceptors());
   // Optional columns.
+  cHaveTsys = False;
+  cHaveTcal = False;
   if ((cHaveSrcVel = cPKSMS.source().tableDesc().isColumn("SYSVEL"))) {
     cSrcVelCol.attach(cPKSMS.source(), "SYSVEL");
+  }
   if ((cHaveTsys = cPKSMS.sysCal().tableDesc().isColumn("TSYS"))) {
+  if (cHaveSysCal && (cHaveTsys = cPKSMS.sysCal().tableDesc().isColumn("TSYS"))) {
     cTsysCol.attach(cPKSMS.sysCal(), "TSYS");
+  }
+  if (cHaveSysCal && (cHaveTcal = cPKSMS.sysCal().tableDesc().isColumn("TCAL"))) {
+    cTcalCol.attach(cPKSMS.sysCal(), "TCAL");
+  }
 …
   // Spectral data should always be present.
   haveSpectra = True;
+  cFloatDataCol.reference(msCols.floatData());
+  cHaveDataCol = False;
+  cHaveCorrectedDataCol = False;
+  ROMSObservationColumns observationCols(cPKSMS.observation());
+  //String telName = observationCols.telescopeName()(0);
+  cTelName = observationCols.telescopeName()(0);
+  //cATF = cTelName.contains("ATF");
+  //cOSF = cTelName.contains("OSF");
+  //cALMA = cTelName.contains("ALMA");
+  cALMA = cTelName.contains("ATF")||cTelName.contains("OSF")||
+           cTelName.contains("ALMA");
+  if (cHaveDataCol = cPKSMS.isColumn(MSMainEnums::DATA)) {
+    if (cALMA) {
+      //try to read a single baseline interferometeric data
+      //and treat it as single dish data
+      //maybe extended for ALMA commissioning later
+      cDataCol.reference(msCols.data());
+      if (cHaveCorrectedDataCol = cPKSMS.isColumn(MSMainEnums::CORRECTED_DATA)) {
+        //cerr<<"Do have CORRECTED_DATA column"<<endl;
+        cCorrectedDataCol.reference(msCols.correctedData());
+      }
+    }
+  }
+  else {
+    cFloatDataCol.reference(msCols.floatData());
+  }
   cFlagCol.reference(msCols.flag());
+  if ((cGetXPol = cPKSMS.isColumn(MSMainEnums::DATA))) {
+  cFlagRowCol.reference(msCols.flagRow());
+  if (cGetXPol = (cPKSMS.isColumn(MSMainEnums::DATA) && (!cALMA))) {
     if ((cHaveXCalFctr = cPKSMS.tableDesc().isColumn("XCALFCTR"))) {
       cXCalFctrCol.attach(cPKSMS, "XCALFCTR");
 …
   // Number of IFs.
+  uInt nIF = dataDescCols.nrow();
+  //uInt nIF = dataDescCols.nrow();
+  uInt nIF =spWinCols.nrow();
+  Vector<Int> spWinIds = cSpWinIdCol.getColumn() ;
   IFs.resize(nIF);
   IFs = True;
+  for ( Int ispw = 0 ; ispw < nIF ; ispw++ ) {
+    if ( allNE( ispw, spWinIds ) ) {
+      IFs(ispw) = False ;
+    }
+  }
   // Number of polarizations and channels in each IF.
-  ROScalarColumn<Int> spWinIdCol(dataDescCols.spectralWindowId());
   ROScalarColumn<Int> numChanCol(spWinCols.numChan());
 …
   nPol.resize(nIF);
   for (uInt iIF = 0; iIF < nIF; iIF++) {
+    nChan(iIF) = numChanCol(spWinIdCol(iIF));
+    nPol(iIF)  = numPolCol(polIdCol(iIF));
+    if ( IFs(iIF) ) {
+      nChan(iIF) = numChanCol(cSpWinIdCol(iIF)) ;
+      nPol(iIF) = numPolCol(polIdCol(iIF)) ;
+    }
+    else {
+      nChan(iIF) = 0 ;
+      nPol(iIF) = 0 ;
+    }
+  }
 …
   haveXPol = False;
   if (cGetXPol) {
+  if (cGetXPol && !(cALMA)) {
     for (Int irow = 0; irow < cNRow; irow++) {
       if (cDataCol.isDefined(irow)) {
 …
         String &antName,
         Vector<Double> &antPosition,
+        // before merge...
+        //String &obsMode,
         String &obsType,
         String &bunit,
 …
         Double &mjd,
         Double &refFreq,
         Double &bandwidth)
+        Double &bandwidth)
+{
   if (!cMSopen) {
 …
   // Antenna name and ITRF coordinates.
   ROMSAntennaColumns antennaCols(cPKSMS.antenna());
+  antName = antennaCols.name()(0);
+  antPosition = antennaCols.position()(0);
+  //antName = antennaCols.name()(0);
+  antName = antennaCols.name()(cAntId[0]);
+  if (cALMA) {
+     antName = cTelName + "-" + antName;
+  }
+  //antPosition = antennaCols.position()(0);
+  antPosition = antennaCols.position()(cAntId[0]);
   // Observation type.
 …
+  }
+  // Brightness units.
+  bunit = cPKSMS.unit(MSMainEnums::FLOAT_DATA);
+  bunit = "";
+  if (cHaveDataCol) {
+    const TableRecord& keywordSet2
+       = cDataCol.columnDesc().keywordSet();
+    if(keywordSet2.isDefined("UNIT")) {
+      bunit = keywordSet2.asString("UNIT");
+    }
+  } else {
+    const TableRecord& keywordSet
+       = cFloatDataCol.columnDesc().keywordSet();
+    if(keywordSet.isDefined("UNIT")) {
+      bunit = keywordSet.asString("UNIT");
+    }
+  }
+/***
+  const TableRecord& keywordSet
+       = cFloatDataCol.columnDesc().keywordSet();
+  if(keywordSet.isDefined("UNIT")) {
+    fluxunit = keywordSet.asString("UNIT");
+  }
+***/
   // Coordinate equinox.
   ROMSPointingColumns pointingCols(cPKSMS.pointing());
   String dirref = pointingCols.direction().keywordSet().asRecord("MEASINFO").
                     asString("Ref");
+  cDirRef = dirref;
+  if (dirref =="AZELGEO" || dirref == "AZEL") {
+     dirref = "J2000";
+  }
   sscanf(dirref.chars()+1, "%f", &equinox);
 …
         const Bool getSpectra,
         const Bool getXPol,
+        const Bool getFeedPos,
+        const Bool getPointing,
         const Int  coordSys)
+{
 …
   cGetXPol = cGetXPol && getXPol;
+  // Get feed positions?  (Not available.)
+  cGetFeedPos = False;
+  // Get Pointing data (for MS)
+  cGetPointing = getPointing;
   // Coordinate system?  (Only equatorial available.)
   cCoordSys = 0;
 …
 // Read the next data record.
+/**
+Int PKSMS2reader::read(
+        Int             &scanNo,
+        Int             &cycleNo,
+        Double          &mjd,
+        Double          &interval,
+        String          &fieldName,
+        String          &srcName,
+        Vector<Double>  &srcDir,
+        Vector<Double>  &srcPM,
+        Double          &srcVel,
+        String          &obsMode,
+        Int             &IFno,
+        Double          &refFreq,
+        Double          &bandwidth,
+        Double          &freqInc,
+        Vector<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,
+        uInt            &flagrow,
+        Complex         &xCalFctr,
+        Vector<Complex> &xPol)
+**/
 Int PKSMS2reader::read(PKSrecord &pksrec)
+{
+  LogIO os( LogOrigin( "PKSMS2reader", "read()", WHERE ) ) ;
   if (!cMSopen) {
     return 1;
 …
   Int ibeam;
   Int iIF;
+  Int iDataDesc;
   while (True) {
     ibeam = cBeamNoCol(cIdx);
+    iIF   = cDataDescIdCol(cIdx);
+    iDataDesc   = cDataDescIdCol(cIdx);
+    iIF   =cSpWinIdCol(iDataDesc);
     if (cBeams(ibeam) && cIFs(iIF)) {
       break;
 …
+    }
+  }
   // Renumerate scan no. Here still is 1-based
+  pksrec.scanNo = cScanNoCol(cIdx) - cScanNoCol(0) + 1;
+  //scanNo = cScanNoCol(cIdx) - cScanNoCol(0) + 1;
+  //scanNo = cScanNoCol(cIdx);
+  pksrec.scanNo = cScanNoCol(cIdx);
   if (pksrec.scanNo != cScanNo) {
 …
   Int srcId = cSrcIdCol(fieldId);
+  pksrec.srcName = cSrcNameCol(srcId);
+  //For source with multiple spectral window setting, this is not
+  // correct. Source name of srcId may not be at 'srcId'th row of SrcNameCol
+  //srcName = cSrcNameCol(srcId);
+  for (uInt irow = 0; irow < cSrcId2Col.nrow(); irow++) {
+    if (cSrcId2Col(irow) == srcId) {
+      //srcName = cSrcNameCol(irow);
+      pksrec.srcName = cSrcNameCol(irow);
+    }
+  }
   pksrec.srcDir  = cSrcDirCol(srcId);
   pksrec.srcPM   = cSrcPMCol(srcId);
   // Systemic velocity.
   if (!cHaveSrcVel) {
+  if (!cHaveSrcVel || cALMA) {
     pksrec.srcVel = 0.0f;
   } else {
 …
+  }
+  ROMSAntennaColumns antennaCols(cPKSMS.antenna());
+  //String telescope = antennaCols.name()(0);
+  String telescope = antennaCols.name()(cAntId[0]);
+  Bool cGBT = telescope.contains("GBT");
+  //Bool cPM = telescope.contains("PM"); // ACA TP antenna
+  //Bool cDV = telescope.contains("DV"); // VERTEX
+  //Bool cCM = telescope.contains("CM"); // ACA 7m antenna
+  //Bool cALMA = cPM || cDV || cCM ;
   // Observation type.
+  Int stateId = cStateIdCol(cIdx);
+  pksrec.obsType = cObsModeCol(stateId);
+  // check if State Table exist
+  //Bool cHaveStateTab=Table::isReadable(cPKSMS.stateTableName());
+  Int stateId = 0;
+  Int StateNRow = 0;
+  StateNRow=cObsModeCol.nrow();
+  if (Table::isReadable(cPKSMS.stateTableName())) {
+        pksrec.obsType = " ";
+    if (StateNRow > 0) {
+      stateId = cStateIdCol(cIdx);
+      if (stateId == -1) {
+        //pksrec.obsType = " ";
+      } else {
+        pksrec.obsType = cObsModeCol(stateId);
+        Bool sigState =cSigStateCol(stateId);
+        Bool refState =cRefStateCol(stateId);
+        //DEBUG
+        //cerr <<"stateid="<<stateId<<" obsmode="<<pksrec.obsType<<endl;
+        if (cGBT) {
+          // split the obsType string and append a proper label
+          // (these are GBT specific)
+          int epos = pksrec.obsType.find_first_of(':');
+          int nextpos = pksrec.obsType.find_first_of(':',epos+1);
+          string obsMode1 = pksrec.obsType.substr(0,epos);
+          string obsMode2 = pksrec.obsType.substr(epos+1,nextpos-epos-1);
+          //cerr <<"obsMode2= "<<obsMode2<<endl;
+          if (!pksrec.srcName.contains("_ps")
+              &&!pksrec.srcName.contains("_psr")
+              &&!pksrec.srcName.contains("_nod")
+              &&!pksrec.srcName.contains("_fs")
+              &&!pksrec.srcName.contains("_fsr")) {
+            // if Nod mode observation , append '_nod'
+            if (obsMode1 == "Nod") {
+              //pksrec.srcName.append("_nod");
+              pksrec.srcType = SrcType::NOD ;
+            } else if (obsMode1 == "OffOn") {
+            // for GBT position switch observations (OffOn or OnOff)
+              //if (obsMode2 == "PSWITCHON") pksrec.srcName.append("_ps");
+              //if (obsMode2 == "PSWITCHOFF") pksrec.srcName.append("_psr");
+              if (obsMode2 == "PSWITCHON") pksrec.srcType = SrcType::PSON ;
+              if (obsMode2 == "PSWITCHOFF") pksrec.srcType = SrcType::PSOFF ;
+            } else {
+              if (obsMode2 == "FSWITCH") {
+              // for GBT frequency switch mode
+                //if (sigState) pksrec.srcName.append("_fs");
+                //if (refState) pksrec.srcName.append("_fsr");
+                if (sigState) pksrec.srcType = SrcType::FSON ;
+                if (refState) pksrec.srcType = SrcType::FSOFF ;
+              }
+            }
+          }
+        }
+        else if (cALMA) {
+          // ALMA tag
+          // split the obsType string and append a proper label
+          string substr[1] ;
+          int numSubstr = split( pksrec.obsType, substr, 1, "," );
+          String obsType = String( substr[0] );
+          int epos = obsType.find_first_of('.');
+          int nextpos = obsType.find_first_of('.',epos+1);
+          string obsMode1 = obsType.substr(0,epos);
+          string obsMode2 = obsType.substr(epos+1,nextpos-epos-1);
+          //cerr <<"obsMode2= "<<obsMode2<<endl;
+          // Current OBS_MODE format:
+          //
+          //     ON: OBSERVE_TARGET.ON_SOURCE
+          //    OFF: OBSERVE_TARGET.OFF_SOURCE
+          //
+          if (obsMode1 == "OBSERVE_TARGET") {
+            //if (obsMode2 == "ON_SOURCE") pksrec.srcName.append("_pson");
+            //if (obsMode2 == "OFF_SOURCE") pksrec.srcName.append("_psoff");
+            if (obsMode2 == "ON_SOURCE") pksrec.srcType = SrcType::PSON ;
+            if (obsMode2 == "OFF_SOURCE") pksrec.srcType = SrcType::PSOFF ;
+          }
+        }
+      }
+    }
+  }
+  // CAL state
+  // this should be apply just for GBT data?
+  Double Cal;
+  if (stateId==-1 || StateNRow==0) {
+    Cal = 0;
+  } else {
+    Cal = cCalCol(stateId);
+  }
+  if (cGBT) {
+    if (Cal > 0 && !pksrec.srcName.contains("_calon")) {
+      //pksrec.srcName.append("_calon");
+      if ( pksrec.srcType == SrcType::NOD )
+        pksrec.srcType = SrcType::NODCAL ;
+      else if ( pksrec.srcType == SrcType::PSON )
+        pksrec.srcType = SrcType::PONCAL ;
+      else if ( pksrec.srcType == SrcType::PSOFF )
+        pksrec.srcType = SrcType::POFFCAL ;
+      else if ( pksrec.srcType == SrcType::FSON )
+        pksrec.srcType = SrcType::FONCAL ;
+      else if ( pksrec.srcType == SrcType::FSOFF )
+        pksrec.srcType = SrcType::FOFFCAL ;
+      else
+        pksrec.srcName.append("_calon");
+    }
+  }
   pksrec.IFno = iIF + 1;
   Int nChan = abs(cEndChan(iIF) - cStartChan(iIF)) + 1;
   // Minimal handling on continuum data.
   Vector<Double> chanFreq = cChanFreqCol(iIF);
+  pksrec.nchan = nChan;
   if (nChan == 1) {
     cout << "The input is continuum data. "<< endl;
     pksrec.freqInc  = chanFreq(0);
+    //pksrec.freqInc  = chanFreq(0);
+    pksrec.freqInc  = cTotBWCol(iIF);
     pksrec.refFreq  = chanFreq(0);
+    pksrec.restFreq = 0.0f;
+    pksrec.restFreq.resize(1);
+    pksrec.restFreq[0] = 0.0f;
   } else {
     if (cStartChan(iIF) <= cEndChan(iIF)) {
       pksrec.freqInc = chanFreq(1) - chanFreq(0);
 …
     pksrec.refFreq  = chanFreq(cRefChan(iIF)-1);
+    pksrec.restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0));
+  }
+  pksrec.bandwidth = abs(pksrec.freqInc * nChan);
+    Bool HaveSrcRestFreq= cSrcRestFrqCol.isDefined(srcId);
+    if (HaveSrcRestFreq) {
+      //restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0));
+      //restFreq = cSrcRestFrqCol(srcId);
+      pksrec.restFreq = cSrcRestFrqCol(srcId);
+    } else {
+      pksrec.restFreq.resize(1);
+      pksrec.restFreq[0] = 0.0f;
+    }
+  }
+  //pksrec.bandwidth = abs(pksrec.freqInc * nChan);
+  pksrec.bandwidth = abs(cTotBWCol(0));
   pksrec.tcal.resize(cNPol(iIF));
   pksrec.tcal      = 0.0f;
   pksrec.tcalTime  = "";
   pksrec.azimuth   = 0.0f;
   pksrec.elevation = 0.0f;
+//  pksrec.azimuth   = 0.0f;
+//  pksrec.elevation = 0.0f;
   pksrec.parAngle  = 0.0f;
 …
   // Find the appropriate entry in the WEATHER subtable.
+  Vector<Double> wTimes = cWeatherTimeCol.getColumn();
+  Int weatherIdx;
+  for (weatherIdx = wTimes.nelements()-1; weatherIdx >= 0; weatherIdx--) {
+    if (cWeatherTimeCol(weatherIdx) <= time) {
+      break;
+    }
+  }
+  if (weatherIdx < 0) {
+  //Bool cHaveStateTab=Table::isReadable(cPKSMS.stateTableName());
+  Bool cHaveWeatherTab = Table::isReadable(cPKSMS.weatherTableName());
+  Int weatherIdx=-1;
+  if (cHaveWeatherTab) {
+    Vector<Double> wTimes = cWeatherTimeCol.getColumn();
+    for (weatherIdx = wTimes.nelements()-1; weatherIdx >= 0; weatherIdx--) {
+      if (cWeatherTimeCol(weatherIdx) <= time) {
+        break;
+      }
+    }
+  }
+  if (weatherIdx < 0 || !cHaveWeatherTab) {
     // No appropriate WEATHER entry.
     pksrec.temperature = 0.0f;
 …
   pksrec.beamNo  = ibeam + 1;
+  Matrix<Double> pointingDir = cPointingCol(fieldId);
+  pksrec.direction = pointingDir.column(0);
+  //pointing/azel
+  MVPosition mvpos(antennaCols.position()(cAntId[0]));
+  MPosition mp(mvpos);
+  Quantum<Double> qt(time,"s");
+  MVEpoch mvt(qt);
+  MEpoch me(mvt);
+  MeasFrame frame(mp, me);
+  MDirection md;
   pksrec.pCode = 0;
   pksrec.rateAge = 0.0f;
-  uInt ncols = pointingDir.ncolumn();
-  if (ncols == 1) {
-    pksrec.scanRate = 0.0f;
-  } else {
-    pksrec.scanRate(0) = pointingDir.column(1)(0);
-    pksrec.scanRate(1) = pointingDir.column(1)(1);
+  }
   pksrec.paRate = 0.0f;
+  if (cGetPointing) {
+    //cerr << "get pointing data ...." << endl;
+    ROScalarColumn<Int> pAntIdCol ;
+    ROScalarColumn<Double> psTimeCol ;
+    Table ptTable = cPKSMS.pointing() ;
+    MSPointing selPtTab( ptTable( ptTable.col("ANTENNA_ID") == cAntId[0] ) ) ;
+    pAntIdCol.attach( selPtTab, "ANTENNA_ID" ) ;
+    Vector<Int> antIds = pAntIdCol.getColumn() ;
+    psTimeCol.attach( selPtTab, "TIME" ) ;
+    Vector<Double> pTimes = psTimeCol.getColumn();
+    Bool doInterp = False ;
+    Int PtIdx=-1;
+    for (PtIdx = pTimes.nelements()-1; PtIdx >= 0; PtIdx--) {
+      if ( pTimes[PtIdx] == time ) {
+        break ;
+      }
+      else if ( pTimes[PtIdx] < time ) {
+        if ( PtIdx != pTimes.nelements()-1 ) {
+          doInterp = True ;
+        }
+        break ;
+      }
+    }
+    if ( PtIdx == -1 ) {
+      PtIdx = 0 ;
+    }
+    //cerr << "got index=" << PtIdx << endl;
+    Matrix<Double> pointingDir = cPointingCol(PtIdx);
+    ROMSPointingColumns PtCols( selPtTab ) ;
+    Vector<Double> pointingDirVec ;
+    if ( doInterp ) {
+      Double dt1 = time - pTimes[PtIdx] ;
+      Double dt2 = pTimes[PtIdx+1] - time ;
+      Vector<Double> dirVec1 = pointingDir.column(0) ;
+      Matrix<Double> pointingDir2 = cPointingCol(PtIdx+1) ;
+      Vector<Double> dirVec2 = pointingDir2.column(0) ;
+      pointingDirVec = (dt1*dirVec2+dt2*dirVec1)/(dt1+dt2) ;
+      Vector<MDirection> vmd1(1) ;
+      Vector<MDirection> vmd2(1) ;
+      PtCols.directionMeasCol().get(PtIdx,vmd1) ;
+      Vector<Double> angle1 = vmd1(0).getAngle().getValue("rad") ;
+      PtCols.directionMeasCol().get(PtIdx+1,vmd2) ;
+      Vector<Double> angle2 = vmd2(0).getAngle().getValue("rad") ;
+      Vector<Double> angle = (dt1*angle2+dt2*angle1)/(dt1+dt2) ;
+      Quantum< Vector<Double> > qangle( angle, "rad" ) ;
+      String typeStr = vmd1(0).getRefString() ;
+      //cerr << "vmd1.getRefString()=" << typeStr << endl ;
+      MDirection::Types mdType ;
+      MDirection::getType( mdType, typeStr ) ;
+      //cerr << "mdType=" << mdType << endl ;
+      md = MDirection( qangle, mdType ) ;
+      //cerr << "md=" << md.getAngle().getValue("rad") << endl ;
+    }
+    else {
+      pointingDirVec = pointingDir.column(0) ;
+      Vector<MDirection> vmd(1);
+      PtCols.directionMeasCol().get(PtIdx,vmd);
+      md = vmd[0];
+    }
+    // put J2000 coordinates in "direction"
+    if (cDirRef =="J2000") {
+      pksrec.direction = pointingDirVec ;
+    }
+    else {
+      pksrec.direction =
+        MDirection::Convert(md, MDirection::Ref(MDirection::J2000,
+                                                frame)
+                            )().getAngle("rad").getValue();
+    }
+    uInt ncols = pointingDir.ncolumn();
+    pksrec.scanRate.resize(2);
+    if (ncols == 1) {
+      pksrec.scanRate = 0.0f;
+    } else {
+      pksrec.scanRate(0) = pointingDir.column(1)(0);
+      pksrec.scanRate(1) = pointingDir.column(1)(1);
+    }
+  }
+  else {
+  // Get direction from FIELD table
+  // here, assume direction to be the field direction not pointing
+    Matrix<Double> delayDir = cFieldDelayDirCol(fieldId);
+    pksrec.direction = delayDir.column(0);
+    uInt ncols = delayDir.ncolumn();
+    pksrec.scanRate.resize(2);
+    if (ncols == 1) {
+      pksrec.scanRate = 0.0f;
+    } else {
+      pksrec.scanRate(0)  = delayDir.column(1)(0);
+      pksrec.scanRate(1)  = delayDir.column(1)(1);
+    }
+  }
+  // caluculate azimuth and elevation
+  // first, get the reference frame
+ /**
+  MVPosition mvpos(antennaCols.position()(0));
+  MPosition mp(mvpos);
+  Quantum<Double> qt(time,"s");
+  MVEpoch mvt(qt);
+  MEpoch me(mvt);
+  MeasFrame frame(mp, me);
+  **/
+  //
+  ROMSFieldColumns fldCols(cPKSMS.field());
+  Vector<MDirection> vmd(1);
+  //MDirection md;
+  fldCols.delayDirMeasCol().get(fieldId,vmd);
+  md = vmd[0];
+  //Vector<Double> dircheck = md.getAngle("rad").getValue();
+  //cerr<<"dircheck="<<dircheck<<endl;
+  Vector<Double> azel =
+        MDirection::Convert(md, MDirection::Ref(MDirection::AZEL,
+                                                frame)
+                            )().getAngle("rad").getValue();
+  //cerr<<"azel="<<azel<<endl;
+  pksrec.azimuth = azel[0];
+  pksrec.elevation = azel[1];
   // Get Tsys assuming that entries in the SYSCAL table match the main table.
 …
   cSigmaCol.get(cIdx, pksrec.sigma, True);
+  //get Tcal if available
+  if (cHaveTcal) {
+    Int nTcalColRow = cTcalCol.nrow();
+    uInt nBeam = cBeams.nelements();
+    uInt nIF = cIFs.nelements();
+    uInt nrws = nBeam * nIF;
+    if (nTcalColRow > 0) {
+    // find tcal match with the data with the data time stamp
+      Double mjds = pksrec.mjd*(24*3600);
+      Double dtcalTime;
+      if ( pksrec.mjd > lastmjd || cIdx==0 ) {
+        //Table tmptab = cSysCalTab(near(cSysCalTab.col("TIME"),mjds));
+        tmptab = cSysCalTab(near(cSysCalTab.col("TIME"),mjds), nrws);
+        //DEBUG
+        //if (cIdx == 0) {
+        //  cerr<<"inital table retrieved"<<endl;
+        //}
+      }
+      if (nBeam == 1) {
+        tmptab2 = tmptab( tmptab.col("SPECTRAL_WINDOW_ID") == iIF, 1);
+      } else {
+        tmptab2 = tmptab( tmptab.col("SPECTRAL_WINDOW_ID") == iIF &&
+                              tmptab.col("FEED_ID") == ibeam , 1);
+      }
+      //cerr<<"first subtab rows="<<tmptab.nrow()<<endl;
+      int syscalrow = tmptab2.nrow();
+      ROArrayColumn<Float> tcalCol(tmptab2, "TCAL");
+      ROScalarColumn<Double> tcalTimeCol(tmptab2, "TIME");
+      if (syscalrow==0) {
+        os << LogIO::NORMAL
+           <<"Cannot find any matching Tcal at/near the data timestamp."
+           << " Set Tcal=0.0" << LogIO::POST ;
+      } else {
+        tcalCol.get(0, pksrec.tcal);
+        tcalTimeCol.get(0,dtcalTime);
+        pksrec.tcalTime = MVTime(dtcalTime/(24*3600)).string(MVTime::YMD);
+        //DEBUG
+        //cerr<<"cIdx:"<<cIdx<<" tcal="<<tcal<<" tcalTime="<<tcalTime<<endl;
+        tmptab.markForDelete();
+        tmptab2.markForDelete();
+      }
+    }
+    lastmjd = pksrec.mjd;
+  }
   // Calibration factors (if available).
   pksrec.calFctr.resize(cNPol(iIF));
 …
     Matrix<Float> tmpData;
     Matrix<Bool>  tmpFlag;
+    cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
+    if (cHaveDataCol) {
+      Matrix<Complex> tmpCmplxData;
+      Matrix<Float> tmpReData;
+      Matrix<Float> tmpImData;
+      //cerr<<"reading spectra..."<<endl;
+      //# TODO - should have a flag to user to select DATA or CORRECTED_DATA
+      //# currently just automatically determined, --- read CORRECTED one
+      //# if the column exist.
+      if (cHaveCorrectedDataCol) {
+        cCorrectedDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True);
+      } else {
+        cDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True);
+      }
+      tmpReData = real(tmpCmplxData);
+      tmpImData = imag(tmpCmplxData);
+      tmpData = sqrt(tmpReData*tmpReData + tmpImData*tmpImData);
+    } else {
+      cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
+    }
     cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
 …
+      }
+    }
+    // Row-based flagging info. (True:1, False:0)
+    pksrec.flagrow = (cFlagRowCol(cIdx) ? 1 : 0);
+  }
   // Get cross-polarization data.
   if (cGetXPol) {
+    //cerr<<"cGetXPol="<<cGetXPol<<endl;
+    //cerr<<"cHaveXCalFctr="<<cHaveXCalFctr<<endl;
     if (cHaveXCalFctr) {
       cXCalFctrCol.get(cIdx, pksrec.xCalFctr);
 …
+    }
+    cDataCol.get(cIdx, pksrec.xPol, True);
+    if (cEndChan(iIF) < cStartChan(iIF)) {
+      Complex ctmp;
+    if(!cALMA) {
+      cDataCol.get(cIdx, pksrec.xPol, True);
+      if (cEndChan(iIF) < cStartChan(iIF)) {
+        Complex ctmp;
+        Int jchan = nChan - 1;
+        for (Int ichan = 0; ichan < nChan/2; ichan++, jchan--) {
+          ctmp = pksrec.xPol(ichan);
+          pksrec.xPol(ichan) = pksrec.xPol(jchan);
+          pksrec.xPol(jchan) = ctmp;
+        }
+      }
+    }
+  }
+  /**
+  cerr<<"scanNo="<<scanNo<<endl;
+  cerr<<"cycleNo="<<cycleNo<<endl;
+  cerr<<"mjd="<<mjd<<endl;
+  cerr<<"interval="<<interval<<endl;
+  cerr<<"fieldName="<<fieldName<<endl;
+  cerr<<"srcNmae="<<srcName<<endl;
+  cerr<<"srcDir="<<srcDir<<endl;
+  cerr<<"srcPM="<<srcPM<<endl;
+  cerr<<"srcVel="<<srcVel<<endl;
+  cerr<<"obsMode="<<obsMode<<endl;
+  cerr<<"IFno="<<IFno<<endl;
+  cerr<<"refFreq="<<refFreq<<endl;
+  cerr<<"tcal="<<tcal<<endl;
+  cerr<<"direction="<<direction<<endl;
+  cerr<<"scanRate="<<scanRate<<endl;
+  cerr<<"tsys="<<tsys<<endl;
+  cerr<<"sigma="<<sigma<<endl;
+  cerr<<"calFctr="<<calFctr<<endl;
+  cerr<<"baseLin="<<baseLin<<endl;
+  cerr<<"baseSub="<<baseSub<<endl;
+  cerr<<"spectra="<<spectra.shape()<<endl;
+  cerr<<"flagged="<<flagged.shape()<<endl;
+  cerr<<"xCalFctr="<<xCalFctr<<endl;
+  cerr<<"xPol="<<xPol<<endl;
+  **/
+  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;
+  Int iDataDesc;
+  while (True) {
+    ibeam = cBeamNoCol(cIdx);
+    //iIF   = cDataDescIdCol(cIdx);
+    iDataDesc   = cDataDescIdCol(cIdx);
+    iIF   = cSpWinIdCol(iDataDesc);
+    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(24,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;
+    if (cHaveDataCol) {
+      Matrix<Complex> tmpCmplxData;
+      cDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True);
+      tmpData = real(tmpCmplxData);
+    } else {
+      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 ichan = 0; ichan < nChan/2; ichan++, jchan--) {
+        ctmp = pksrec.xPol(ichan);
+        pksrec.xPol(ichan) = pksrec.xPol(jchan);
+        pksrec.xPol(jchan) = ctmp;
+      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);
+        }
+      }
+    }
 …
   cMSopen = False;
+}
+//-------------------------------------------------------- PKSMS2reader::splitAntenanSelectionString
+// split antenna selection string
+// delimiter is ','
+Vector<String> PKSMS2reader::splitAntennaSelectionString( const String s )
+{
+  Char delim = ',' ;
+  Int n = s.freq( delim ) + 1 ;
+  Vector<String> antlist ;
+  string sl[n] ;
+  Int numSubstr = split( s, sl, n, "," );
+  antlist.resize( numSubstr ) ;
+  for ( Int i = 0 ; i < numSubstr ; i++ ) {
+    antlist[i] = String( sl[i] ) ;
+    antlist[i].trim() ;
+  }
+  //cerr << "antlist = " << antlist << endl ;
+  return antlist ;
+}
+//-------------------------------------------------------- PKSMS2reader::setupAntennaList
+// Fill cAntenna and cAntId
+void PKSMS2reader::setupAntennaList( const String s )
+{
+  LogIO os( LogOrigin( "PKSMS2reader", "setupAntennaList()", WHERE ) ) ;
+  //cerr << "antenna specification: " << s << endl ;
+  ROMSAntennaColumns antennaCols(cPKSMS.antenna());
+  ROScalarColumn<String> antNames = antennaCols.name();
+  Int nrow = antNames.nrow() ;
+  Vector<String> antlist = splitAntennaSelectionString( s ) ;
+  Int len = antlist.size() ;
+  Vector<Int> AntId( len ) ;
+  Regex re( "[0-9]+" ) ;
+  for ( Int i = 0 ; i < len ; i++ ) {
+    if ( antlist[i].matches( re ) ) {
+      AntId[i] = atoi( antlist[i].c_str() ) ;
+      if ( AntId[i] >= nrow ) {
+        os << LogIO::SEVERE << "Antenna index out of range: " << AntId[i] << LogIO::EXCEPTION ;
+      }
+    }
+    else {
+      AntId[i] = -1 ;
+      for ( uInt j = 0 ; j < antNames.nrow() ; j++ ) {
+        if ( antlist[i] == antNames(j) ) {
+          AntId[i] = j ;
+          break ;
+        }
+      }
+      if ( AntId[i] == -1 ) {
+        os << LogIO::SEVERE << "Specified antenna name not found: " << antlist[i] << LogIO::EXCEPTION ;
+      }
+    }
+  }
+  //cerr << "AntId = " << AntId << endl ;
+  vector<Int> uniqId ;
+  uniqId.push_back( AntId(0) ) ;
+  for ( uInt i = 1 ; i < AntId.size() ; i++ ) {
+    if ( count(uniqId.begin(),uniqId.end(),AntId[i]) == 0 ) {
+      uniqId.push_back( AntId[i] ) ;
+    }
+  }
+  Vector<Int> newAntId( uniqId ) ;
+  cAntId.assign( newAntId ) ;
+  //cerr << "cAntId = " << cAntId << endl ;
+}

branches/mergetest/external/atnf/PKSIO/PKSMS2reader.h

-              r1720
+              r1779
     virtual Int open(
         const String msName,
+        const String antenna,
         Vector<Bool> &beams,
         Vector<Bool> &IFs,
 …
         String &bunit,
         Float  &equinox,
+        //String &freqRef,
         String &dopplerFrame,
         Double &mjd,
 …
         const Bool getSpectra = True,
         const Bool getXPol    = False,
+        const Bool getFeedPos = False,
+        const Bool getPointing = False,
         const Int  coordSys   = 0);
     // Find the range of the data selected in time and position.
 …
     // Read the next data record.
+/**
+    virtual Int read(
+        Int             &scanNo,
+        Int             &cycleNo,
+        Double          &mjd,
+        Double          &interval,
+        String          &fieldName,
+        String          &srcName,
+        Vector<Double>  &srcDir,
+        Vector<Double>  &srcPM,
+        Double          &srcVel,
+        String          &obsMode,
+        Int             &IFno,
+        Double          &refFreq,
+        Double          &bandwidth,
+        Double          &freqInc,
+        Vector<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);
+**/
     virtual Int read(PKSrecord &pksrec);
+    // 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);
     // Close the MS.
 …
   private:
+    Vector<String> splitAntennaSelectionString( const String s );
+    void setupAntennaList( const String s ) ;
     Bool   cHaveBaseLin, cHaveCalFctr, cHaveSrcVel, cHaveTsys, cHaveXCalFctr,
            cMSopen;
+           cMSopen, cHaveTcal, cHaveDataCol, cALMA, cHaveSysCal, cHaveCorrectedDataCol;
     Int    cCycleNo, cIdx, cNRow, cScanNo;
     Double cTime;
+    Double cTime, lastmjd;
     Vector<Int>    cEndChan, cRefChan, cStartChan;
     Vector<Bool>   cBeams, cIFs;
     Vector<Slicer> cDataSel;
+    String         cDirRef, cTelName;
     MeasurementSet cPKSMS;
+    Table          cSysCalTab, tmptab, tmptab2;
+    //Vector<String> cAntenna;
+    Vector<Int> cAntId;
     ROScalarColumn<Int>     cScanNoCol;
 …
     ROScalarColumn<Int>     cFieldIdCol;
     ROScalarColumn<String>  cFieldNameCol;
+    ROArrayColumn<Double>   cFieldDelayDirCol;
     ROScalarColumn<Int>     cSrcIdCol;
+    ROScalarColumn<Int>     cSrcId2Col;
     ROScalarColumn<String>  cSrcNameCol;
     ROArrayColumn<Double>   cSrcDirCol;
 …
     ROArrayColumn<Double>   cSrcVelCol;
     ROScalarColumn<Int>     cStateIdCol;
+    ROScalarColumn<Double>  cCalCol;
     ROScalarColumn<String>  cObsModeCol;
     ROArrayColumn<Double>   cSrcRestFrqCol;
     ROScalarColumn<Int>     cDataDescIdCol;
+    ROScalarColumn<Int>     cSpWinIdCol;
     ROArrayColumn<Double>   cChanFreqCol;
+    ROScalarColumn<Double>   cTotBWCol;
     ROScalarColumn<Double>  cWeatherTimeCol;
     ROScalarColumn<Float>   cTemperatureCol;
     ROScalarColumn<Float>   cPressureCol;
     ROScalarColumn<Float>   cHumidityCol;
+    ROArrayColumn<Float>    cTcalCol;
     ROScalarColumn<Int>     cBeamNoCol;
     ROArrayColumn<Double>   cPointingCol;
+    ROScalarColumn<Double>  cPointingTimeCol;
     ROArrayColumn<Float>    cTsysCol;
     ROArrayColumn<Float>    cSigmaCol;
 …
     ROArrayColumn<Float>    cFloatDataCol;
     ROArrayColumn<Bool>     cFlagCol;
+    ROScalarColumn<Bool>    cFlagRowCol;
     ROScalarColumn<Complex> cXCalFctrCol;
     ROArrayColumn<Complex>  cDataCol;
+    ROArrayColumn<Complex>  cCorrectedDataCol;
     ROScalarColumn<Int>     cNumReceptorCol;
+    ROScalarColumn<Bool>    cSigStateCol;
+    ROScalarColumn<Bool>    cRefStateCol;
 };

branches/mergetest/external/atnf/PKSIO/PKSMS2writer.cc

-              r1720
+              r1779
 #include <casa/BasicSL/Constants.h>
 #include <casa/Quanta/QC.h>
+#include <casa/Logging/LogIO.h>
 #include <measures/Measures/Stokes.h>
 #include <tables/Tables/ArrColDesc.h>
 …
 #include <tables/Tables/TiledShapeStMan.h>
+// Class name
+const string className = "PKSMS2writer" ;
 //------------------------------------------------- PKSMS2writer::PKSMS2writer
 …
+{
   cPKSMS = 0x0;
-  // By default, messages are written to stderr.
-  initMsg();
+}
 …
         const Bool   haveBase)
+{
+  const string methodName = "create()" ;
+  LogIO os( LogOrigin( className, methodName, WHERE ) ) ;
   if (cPKSMS) {
     logMsg("ERROR: Output MS already open, close it first.");
+    os << LogIO::SEVERE << "Output MS already open, close it first." << LogIO::POST ;
     return 1;
+  }
-  // Clear the message stack.
-  clearMsg();
   // Open a MS table.
 …
   Int maxNPol = max(cNPol);
+  cGBT = cAPEX = cSMT = cALMA = cATF = False;
+  String telName = antName;
+  // check if it is GBT data
+  if (antName.contains("GBT")) {
+    cGBT = True;
+  }
+  else if (antName.contains("APEX")) {
+    cAPEX = True;
+  }
+  else if (antName.contains("HHT") || antName.contains("SMT")) {
+    cSMT = True;
+  }
+  else if (antName.contains("ALMA")) {
+    cALMA = True;
+  }
+  else if (antName.contains("ATF")) {
+    cATF = True;
+    telName="ATF";
+  }
   // Add the non-standard CALFCTR column.
   pksDesc.addColumn(ArrayColumnDesc<Float>("CALFCTR", "Calibration factors",
 …
   // Add the optional FLOAT_DATA column.
   MS::addColumnToDesc(pksDesc, MS::FLOAT_DATA, 2);
+  //pksDesc.rwColumnDesc(MS::columnName(MS::FLOAT_DATA)).rwKeywordSet().
+  //              define("UNIT", String("Jy"));
   pksDesc.rwColumnDesc(MS::columnName(MS::FLOAT_DATA)).rwKeywordSet().
                 define("UNIT", bunit);
 …
     MS::addColumnToDesc(pksDesc, MS::DATA, 2);
+    //pksDesc.rwColumnDesc(MS::columnName(MS::DATA)).rwKeywordSet().
+    //            define("UNIT", "Jy");
     pksDesc.rwColumnDesc(MS::columnName(MS::DATA)).rwKeywordSet().
                 define("UNIT", bunit);
 …
   newtab.bindAll(incrStMan, True);
   // Use TiledShapeStMan for the FLOAT_DATA hypercube with tile size 16 kiB.
   TiledShapeStMan tiledStMan("TiledData", IPosition(3,1,128,32));
+  // Use TiledShapeStMan for the FLOAT_DATA hypercube with tile size 1 MB.
+  TiledShapeStMan tiledStMan("TiledData", IPosition(3,1,128,2048));
   newtab.bindColumn(MS::columnName(MS::FLOAT_DATA), tiledStMan);
 …
   } else if (dopplerFrame == "SOURCE") {
     MFrequency::getType(cDopplerFrame, "REST");
+  } else if (dopplerFrame == "LSRK") {
+    MFrequency::getType(cDopplerFrame, "LSRK");
+  }
 …
   addDopplerEntry();
   addFeedEntry();
+  addObservationEntry(observer, project);
+  //addObservationEntry(observer, project);
+  addObservationEntry(observer, project, telName);
   addProcessorEntry();
 …
 // Write the next data record.
+/**
+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<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 uInt            flagrow,
+        const Complex         xCalFctr,
+        const Vector<Complex> &xPol)
+**/
 Int PKSMS2writer::write(
         const PKSrecord &pksrec)
 …
     pksrec.restFreq, pksrec.srcVel);
   // Find or add the obsType to the STATE subtable.
+  // Find or add the obsMode to the STATE subtable.
   Int stateId = addStateEntry(pksrec.obsType);
   // FIELD subtable.
   Vector<Double> scanRate(2);
   scanRate(0) = pksrec.scanRate(0);
   scanRate(1) = pksrec.scanRate(1);
+  //Vector<Double> scanRate(2);
+  //scanRate(0) = pksrec.scanRate(0);
+  //scanRate(1) = pksrec.scanRate(1);
   Int fieldId = addFieldEntry(pksrec.fieldName, time, pksrec.direction,
     scanRate, srcId);
+    pksrec.scanRate, srcId);
   // POINTING subtable.
   addPointingEntry(time, pksrec.interval, pksrec.fieldName, pksrec.direction,
     scanRate);
+    pksrec.scanRate);
   // SYSCAL subtable.
   addSysCalEntry(pksrec.beamNo, iIF, time, pksrec.interval, pksrec.tcal,
     pksrec.tsys);
+    pksrec.tsys, nPol);
   // Handle weather information.
 …
   cMSCols->sigma().put(irow, pksrec.sigma);
+  Vector<Float> weight(1, 1.0f);
+  //Vector<Float> weight(1, 1.0f);
+  Vector<Float> weight(nPol, 1.0f);
   cMSCols->weight().put(irow, weight);
+  //imaging weight
+  //Vector<Float> imagingWeight(nChan);
+  //cMSCols->imagingWeight().put(irow, imagingWeight);
   // Flag information.
   Cube<Bool> flags(nPol, nChan, 1, False);
   cMSCols->flag().put(irow, flags.xyPlane(0));
+  //cMSCols->flag().put(irow, flags.xyPlane(0));
   cMSCols->flagCategory().put(irow, flags);
+  cMSCols->flagRow().put(irow, False);
+  // Row-based flagging info. (True:>0, False:0)
+  cMSCols->flagRow().put(irow, (pksrec.flagrow > 0));
   return 0;
 …
   cSysCal          = MSSysCal();
   cWeather         = MSWeather();
   // Release the main table.
   delete cPKSMS;
   cPKSMS = 0x0;
+  delete cPKSMS;
+  cPKSMS=0x0;
+}
 …
   Int n = cAntenna.nrow() - 1;
+  // do specific things for GBT
   // Data.
+  // plus some more telescopes
   cAntennaCols->name().put(n, antName);
+  cAntennaCols->station().put(n, "ATNF_PARKES");
+  //cAntennaCols->station().put(n, "ATNF_PARKES");
+  if (cGBT) {
+    cAntennaCols->station().put(n, "GREENBANK");
+    cAntennaCols->dishDiameter().put(n, 110.0);
+  }
+  else if (cAPEX) {
+    cAntennaCols->station().put(n, "CHAJNANTOR");
+    cAntennaCols->dishDiameter().put(n, 12.0);
+  }
+  else if (cALMA) {
+    // this needs to be changed in future...
+    cAntennaCols->station().put(n, "CHAJNANTOR");
+    cAntennaCols->dishDiameter().put(n, 12.0);
+  }
+  else if (cATF) {
+    //pad name for the antenna is static...
+    String stname="unknown";
+    if (antName.contains("DV")) {
+       stname="PAD001";
+    }
+    if (antName.contains("DA")) {
+       stname="PAD002";
+    }
+    cAntennaCols->station().put(n, stname);
+    cAntennaCols->dishDiameter().put(n, 12.0);
+  }
+  else if (cSMT) {
+    cAntennaCols->station().put(n, "MT_GRAHAM");
+    cAntennaCols->dishDiameter().put(n, 10.0);
+  }
+  else {
+    cAntennaCols->station().put(n, "ATNF_PARKES");
+    cAntennaCols->dishDiameter().put(n, 64.0);
+  }
   cAntennaCols->type().put(n, "GROUND-BASED");
   cAntennaCols->mount().put(n, "ALT-AZ");
 …
   Vector<Double> antOffset(3, 0.0);
   cAntennaCols->offset().put(n, antOffset);
+  cAntennaCols->dishDiameter().put(n, 64.0);
+  //cAntennaCols->dishDiameter().put(n, 64.0);
+  //if (cGBT) {
+  //  cAntennaCols->dishDiameter().put(n, 110.0);
+  //}
+  //else {
+  //  cAntennaCols->dishDiameter().put(n, 64.0);
+  //}
   // Flags.
   cAntennaCols->flagRow().put(n, False);
 …
         const Int srcId)
+{
+  ROScalarColumn<String> fldn(cField, "NAME");
+  ROScalarColumn<Int> sourceid(cField, "SOURCE_ID");
+  Int n;
+  Int nFld = cField.nrow();
+  for (n = 0; n < nFld; n++) {
+    if (fldn(n) == fieldName && sourceid(n) == srcId) {
+      break;
+    }
+  }
   // Extend the FIELD subtable.
+  cField.addRow();
+  Int n = cField.nrow() - 1;
+  // Data.
+  cFieldCols->name().put(n, fieldName);
+  cFieldCols->code().put(n, "DRIFT");
+  cFieldCols->time().put(n, time);
+  Matrix<Double> track(2, 2);
+  track.column(0) = direction;
+  track.column(1) = scanRate;
+  cFieldCols->numPoly().put(n, 1);
+  cFieldCols->delayDir().put(n, track);
+  cFieldCols->phaseDir().put(n, track);
+  cFieldCols->referenceDir().put(n, track);
+  cFieldCols->sourceId().put(n, srcId);
+  // Flags.
+  cFieldCols->flagRow().put(n, False);
+  if (n == nFld) {
+    cField.addRow();
+    //Int n = cField.nrow() - 1;
+    // Data.
+    cFieldCols->name().put(n, fieldName);
+    if (cGBT) {
+      cFieldCols->code().put(n, " ");
+    }
+    else {
+      cFieldCols->code().put(n, "DRIFT");
+    }
+    cFieldCols->time().put(n, time);
+    //Matrix<Double> track(2, 2);
+    Matrix<Double> track(2, 1);
+    track.column(0) = direction;
+    //track.column(1) = scanRate;
+    cFieldCols->numPoly().put(n, 1);
+    cFieldCols->delayDir().put(n, track);
+    cFieldCols->phaseDir().put(n, track);
+    cFieldCols->referenceDir().put(n, track);
+    cFieldCols->sourceId().put(n, srcId);
+    // Flags.
+    cFieldCols->flagRow().put(n, False);
+  }
   return n;
 …
 Int PKSMS2writer::addObservationEntry(
         const String observer,
+        const String project)
+        const String project,
+        const String antName)
+{
   // Extend the OBSERVATION subtable.
 …
   // Data.
+  cObservationCols->telescopeName().put(n, "Parkes");
+  //cObservationCols->telescopeName().put(n, "Parkes");
+  cObservationCols->telescopeName().put(n, antName);
   Vector<Double> timerange(2, 0.0);
   cObservationCols->timeRange().put(n, timerange);
 …
   Vector<String> log(1, "none");
   cObservationCols->log().put(n, log);
+  cObservationCols->scheduleType().put(n, "ATNF");
+  //cObservationCols->scheduleType().put(n, "ATNF");
+  cObservationCols->scheduleType().put(n, "");
   Vector<String> schedule(1, "Not available");
   cObservationCols->schedule().put(n, schedule);
 …
 //--------------------------------------------- PKSMS2writer::addPointingEntry
+// Modified to fill pointing data if the direction is the pointing direction.
+// So the following comment is no longer true.
 // Add an entry to the POINTING subtable.  This compulsory subtable simply
 …
         const Vector<Double> scanRate)
+{
+  // Extend the POINTING subtable.
+  cPointing.addRow();
+  Int n = cPointing.nrow() - 1;
+  // Keys.
+  cPointingCols->antennaId().put(n, 0);
+  cPointingCols->time().put(n, time);
+  cPointingCols->interval().put(n, interval);
+  // Data.
+  cPointingCols->name().put(n, fieldName);
+  cPointingCols->numPoly().put(n, 1);
+  cPointingCols->timeOrigin().put(n, time);
+  Matrix<Double> track(2, 2);
+  track.column(0) = direction;
+  track.column(1) = scanRate;
+  cPointingCols->direction().put(n, track);
+  cPointingCols->target().put(n, track);
+  cPointingCols->tracking().put(n, True);
+  ROScalarColumn<Double> tms(cPointing, "TIME");
+  Int n;
+  Int ntm = cPointing.nrow();
+  for (n = 0; n < ntm; n++) {
+    if (tms(n) == time) {
+      break;
+    }
+  }
+  if (n == ntm) {
+    // Extend the POINTING subtable.
+    cPointing.addRow();
+    //Int n = cPointing.nrow() - 1;
+    // Keys.
+    cPointingCols->antennaId().put(n, 0);
+    cPointingCols->time().put(n, time);
+    cPointingCols->interval().put(n, interval);
+    // Data.
+    cPointingCols->name().put(n, fieldName);
+    cPointingCols->numPoly().put(n, 1);
+    cPointingCols->timeOrigin().put(n, time);
+    //Matrix<Double> track(2, 2);
+    Matrix<Double> track(2, 1);
+    track.column(0) = direction;
+    //track.column(1) = scanRate;
+    cPointingCols->direction().put(n, track);
+    cPointingCols->target().put(n, track);
+    cPointingCols->tracking().put(n, True);
+  }
   return n;
+}
 …
   // Data.
   Vector<Int> corrType(2);
+  if (nPol == 1) {
+  corrType.resize(1);
+  corrType(0) = Stokes::XX;
+  }
+  else {
+  //Vector<Int> corrType(2);
   corrType(0) = Stokes::XX;
   corrType(1) = Stokes::YY;
+  }
   cPolarizationCols->corrType().put(n, corrType);
   Matrix<Int> corrProduct(2,2,1);
+  if (nPol == 1) {
+    corrProduct.resize(2,1,1);
+    corrProduct(1,0) = 0;
+  }
   if (nPol == 2) {
     corrProduct(1,0) = 0;
 …
         const Vector<Double> direction,
         const Vector<Double> properMotion,
+        const Double restFreq,
+        //const Double restFreq,
+        const Vector<Double> restFreq,
         const Double radialVelocity)
+{
 …
 //  cSourceCols->position().put(n, position);
     cSourceCols->properMotion().put(n, properMotion);
+    Vector<Double> restFrequency(1, restFreq);
+    cSourceCols->restFrequency().put(n, restFrequency);
+//  Vector<Double> restFrequency(1, restFreq);
+//  cSourceCols->restFrequency().put(n, restFrequency);
+    cSourceCols->restFrequency().put(n, restFreq);
     Vector<Double> sysvel(1, radialVelocity);
     cSourceCols->sysvel().put(n, sysvel);
 …
   // Data.
+  cSpWindowCols->name().put(n, "L-band");
+  //cSpWindowCols->name().put(n, "L-band");
+  cSpWindowCols->name().put(n, " ");
   cSpWindowCols->refFrequency().put(n, refFreq);
 …
         const Double interval,
         const Vector<Float> tcal,
+        const Vector<Float> tsys)
+{
+        const Vector<Float> tsys,
+        const Int nPol)
+{
+  LogIO os(LogOrigin("PKSMS2writer", "addSysCalEntry()", WHERE));
   // Extend the SYSCAL subtable.
   cSysCal.addRow();
   Int n = cSysCal.nrow() - 1;
+  //check fo consistency with n pol
+  //here assume size of Tcal vector = npol
+  Vector<Float> inTcal(nPol,0);
+  Int ndim = tcal.shape()(0);
+  Vector<Float> tmpTcal = tcal;
+  if (nPol != ndim) {
+    os << LogIO::WARN
+       << "Found "<< ndim <<" Tcal value(s) for the data with "<<nPol<<" polarization(s)"
+       << "(expecting one Tcal per pol)."<<endl
+       << "First "<< nPol << " Tcal value(s) will be filled." << LogIO::POST;
+    tmpTcal.resize(nPol, True);
+    inTcal = tmpTcal;
+  }
   // Keys.
   cSysCalCols->antennaId().put(n, 0);
 …
   // Data.
+  cSysCalCols->tcal().put(n, tcal);
+  //cSysCalCols->tcal().put(n, tcal);
+  cSysCalCols->tcal().put(n, inTcal);
   cSysCalCols->tsys().put(n, tsys);

branches/mergetest/external/atnf/PKSIO/PKSMS2writer.h

-              r1720
+              r1779
     // 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<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 uInt            flagrow,
+        const Complex         xCalFctr,
+        const Vector<Complex> &xPol);
+**/
     virtual Int write(
         const PKSrecord &pksrec);
 …
     ScalarColumn<Complex> *cXCalFctrCol;
+    // for handling parameters specific to GBT and other telescopes
+    Bool cGBT, cSMT, cAPEX, cALMA, cATF;
     // Add an entry to the ANTENNA subtable.
 …
     Int addObservationEntry(
         const String observer,
+        const String project);
+        const String project,
+        const String antName);
     // Add an entry to the POINTING subtable.
 …
         const Vector<Double> direction,
         const Vector<Double> properMotion,
+        const Double restFreq,
+        //const Double restFreq,
+        const Vector<Double> restFreq,
         const Double radialVelocity);
 …
         const Double interval,
         const Vector<Float> Tcal,
+        const Vector<Float> Tsys);
+        const Vector<Float> Tsys,
+        const Int nPol);
     // Add an entry to the WEATHER subtable.

branches/mergetest/external/atnf/PKSIO/PKSSDwriter.cc

-              r1735
+              r1779
 #include <atnf/PKSIO/PKSSDwriter.h>
+#include <casa/Logging/LogIO.h>
 #include <casa/stdio.h>
 #include <casa/Quanta/MVTime.h>
 …
 #include <cstring>
+// Class name
+const string className = "PKSSDwriter" ;
 //--------------------------------------------------- PKSSDwriter::PKSSDwriter
 …
 PKSSDwriter::PKSSDwriter()
+{
-  // By default, messages are written to stderr.
-  initMsg();
+}
 …
+{
   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;
+}
 …
         const Bool   haveBase)
+{
   // Clear the message stack.
   clearMsg();
+  const string methodName = "create()" ;
+  LogIO os( LogOrigin( className, methodName, WHERE ) ) ;
   double antPos[3];
 …
   cNIF = nChan.nelements();
   if (nPol.nelements() != cNIF || haveXPol.nelements() != cNIF) {
+    cerr << "PKSSDwriter::create: "
+         << "Inconsistent number of IFs for nChan, nPol, and/or haveXPol."
+         << endl;
+    os << LogIO::SEVERE << "Inconsistent number of IFs for nChan, nPol, and/or haveXPol." << LogIO::POST ;
     return 1;
+  }
 …
         (int *)cNPol.getStorage(deleteIt),
         (int *)cHaveXPol.getStorage(deleteIt), (int)cHaveBase, 1);
   logMsg(cSDwriter.getMsg());
   cSDwriter.clearMsg();
+  //logMsg(cSDwriter.getMsg());
+  //cSDwriter.clearMsg();
   if (status) {
     cSDwriter.deleteFile();
 …
         const PKSrecord &pksrec)
+{
+  const string methodName = "write()" ;
+  LogIO os( LogOrigin( className, methodName, WHERE ) ) ;
   // Do basic checks.
   Int IFno = pksrec.IFno;
   uInt iIF = IFno - 1;
   if (IFno < 1 || Int(cNIF) < IFno) {
     cerr << "PKSDwriter::write: "
          << "Invalid IF number " << IFno
          << " (maximum " << cNIF << ")." << endl;
+    os << LogIO::SEVERE
+       << "Invalid IF number " << IFno
+       << " (maximum " << cNIF << ")." << LogIO::POST ;
     return 1;
+  }
 …
   uInt nChan = pksrec.spectra.nrow();
   if (nChan != cNChan(iIF)) {
+    cerr << "PKSDwriter::write: "
+         << "Wrong number of channels for IF " << IFno << "," << endl
+         << "                   "
+    os << LogIO::SEVERE << "Wrong number of channels for IF " << IFno << "," << endl
          << "got " << nChan << " should be " << cNChan(iIF) << "." << endl;
+    os << LogIO::POST ;
     return 1;
+  }
 …
   uInt nPol = pksrec.spectra.ncolumn();
   if (nPol != cNPol(iIF)) {
+    cerr << "PKSDwriter::write: "
+         << "Wrong number of polarizations for IF " << IFno << "," << endl
+         << "                   "
+         << "got " << nPol << " should be " << cNPol(iIF) << "." << endl;
+    os << LogIO::SEVERE << "Wrong number of polarizations for IF " << IFno << "," << endl
+       << "got " << nPol << " should be " << cNPol(iIF) << "." << endl;
+    os << LogIO::POST ;
     return 1;
+  }
   // Extract calendar information from mjd.
+  // Extract calendar information frrom mjd.
   MVTime time(pksrec.mjd);
   Int year  = time.year();
 …
   mbrec.srcRA    = pksrec.srcDir(0);
   mbrec.srcDec   = pksrec.srcDir(1);
+  mbrec.restFreq = pksrec.restFreq;
+  if (pksrec.restFreq.shape()==0) {
+     mbrec.restFreq = 0;
+  }
+  else {
+     mbrec.restFreq = pksrec.restFreq(0);
+  }
   strncpy(mbrec.obsType, (char *)pksrec.obsType.chars(), 16);
 …
   Int status = cSDwriter.write(mbrec);
   logMsg(cSDwriter.getMsg());
   cSDwriter.clearMsg();
+  //logMsg(cSDwriter.getMsg());
+  //cSDwriter.clearMsg();
   if (status) {
     status = 1;
 …
+{
   cSDwriter.close();
   logMsg(cSDwriter.getMsg());
   cSDwriter.clearMsg();
+}
+  //logMsg(cSDwriter.getMsg());
+  //cSDwriter.clearMsg();
+}

branches/mergetest/external/atnf/PKSIO/PKSSDwriter.h

-              r1720
+              r1779
     virtual ~PKSSDwriter();
-    // Set message disposition.
-    virtual Int setMsg(
-        FILE *fd = 0x0);
     // Create the SDFITS file and write static data.
     virtual Int create(

branches/mergetest/external/atnf/PKSIO/PKSreader.cc

-              r1720
+              r1779
         reader = new PKSFITSreader("SDFITS");
       } else {
         // Assume it's MBFITS.
         format = "MBFITS";
         reader = new PKSFITSreader("MBFITS", retry, interpolate);
+      }
+       } else {
+         // Assume it's MBFITS.
+         format = "MBFITS";
+         reader = new PKSFITSreader("MBFITS", retry, interpolate);
+       }
+    }
 …
     format = "UNRECOGNIZED INPUT FORMAT";
+  }
   return reader;
+}
 …
 // Search a list of directories for a Parkes Multibeam dataset and return an
-// appropriate PKSreader for it.
 PKSreader* getPKSreader(
 …
 PKSreader* getPKSreader(
         const String name,
+        const String antenna,
         const Int retry,
         const Int interpolate,
 …
   // Try to open it.
   if (reader) {
     if (reader->open(name, beams, IFs, nChan, nPol, haveXPol, haveBase,
                      haveSpectra)) {
+    if (reader->open(name, antenna, beams, IFs, nChan, nPol, haveXPol,
+                     haveBase, haveSpectra)) {
       format += " OPEN ERROR";
       delete reader;
 …
 //--------------------------------------------------------------- getPKSreader
 // Search a list of directories for a Parkes Multibeam dataset and open an
+// Search a list of directories for a Parkes Multibeam dataset and return an
 // appropriate PKSreader for it.
+PKSreader* getPKSreader(
         const String name,
+PKSreader* getPKSreader(
+        const String name,
+        const String antenna,
         const Vector<String> directories,
         const Int retry,
 …
   // Try to open it.
   if (reader) {
     if (reader->open(name, beams, IFs, nChan, nPol, haveXPol, haveBase,
                      haveSpectra)) {
+    if (reader->open(name, antenna, beams, IFs, nChan, nPol, haveXPol,
+                     haveBase, haveSpectra)) {
       format += " OPEN ERROR";
       delete reader;

branches/mergetest/external/atnf/PKSIO/PKSreader.h

-              r1720
+              r1779
 #define ATNF_PKSREADER_H
-#include <atnf/PKSIO/PKSmsg.h>
 #include <atnf/PKSIO/PKSrecord.h>
+#include <atnf/PKSIO/SrcType.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>
 …
 class PKSreader* getPKSreader(
         const String name,
+        const String antenna,
         const Int retry,
         const Int interpolate,
 …
 class PKSreader* getPKSreader(
         const String name,
+        const String antenna,
         const Vector<String> directories,
         const Int retry,
 …
         Bool   &haveSpectra);
+class PKSreader : public PKSmsg
+class PKSreader
+{
   public:
 …
     virtual Int open(
         const String inName,
+        const String antenna,
         Vector<Bool> &beams,
         Vector<Bool> &IFs,
 …
         const Bool getSpectra = True,
         const Bool getXPol    = False,
+        const Bool getFeedPos = False,
+        const Bool getPointing = False,
         const Int  coordSys   = 0) = 0;
     // Find the range of the data selected in time and position.
 …
         Matrix<Double> &positions) = 0;
+    // Read the next data record.
+    // Read the next data record.
+/**
+    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,
+        Vector<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) = 0;
+**/
+/**
+    // 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;
+**/
     virtual Int read(PKSrecord &pksrec) = 0;
 …
   protected:
     Bool  cGetSpectra, cGetXPol;
+    Bool   cGetFeedPos, cGetSpectra, cGetXPol, cGetPointing;
     Int   cCoordSys;

branches/mergetest/external/atnf/PKSIO/PKSrecord.h

-              r1720
+              r1779
     Double          bandwidth;
     Double          freqInc;
+    Double          restFreq;
+    Int             nchan;
+    Vector<Double>  restFreq;
     Vector<Float>   tcal;
     String          tcalTime;
 …
     Int             pCode;
     Float           rateAge;
     Vector<Float>   scanRate;
+    Vector<Double>  scanRate;
     Float           paRate;
     Vector<Float>   tsys;
 …
     Matrix<Float>   spectra;
     Matrix<uChar>   flagged;
+    uInt            flagrow;
     Complex         xCalFctr;
     Vector<Complex> xPol;
+    Int             polNo ;
+    Int             srcType ;
 };

branches/mergetest/external/atnf/PKSIO/PKSwriter.h

-              r1720
+              r1779
 #define ATNF_PKSWRITER_H
-#include <atnf/PKSIO/PKSmsg.h>
 #include <atnf/PKSIO/PKSrecord.h>
 …
 // </summary>
 class PKSwriter : public PKSmsg
+class PKSwriter
+{
   public:

branches/mergetest/external/atnf/PKSIO/SDFITSreader.cc

-              r1737
+              r1779
 #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>
 …
 const double D2R = PI / 180.0;
+// Class name
+const string className = "SDFITSreader" ;
 //---------------------------------------------------- SDFITSreader::(statics)
 …
   cEndChan   = 0x0;
   cRefChan   = 0x0;
+  // By default, messages are written to stderr.
+  initMsg();
+  cPols      = 0x0;
+}
 …
         int    &extraSysCal)
+{
+  // Clear the message stack.
+  clearMsg();
+  const string methodName = "open()" ;
   if (cSDptr) {
 …
   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;
+  }
 …
       } 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;
 …
     cNAxes = 5;
     if (readDim(DATA, 1, &cNAxes, cNAxis)) {
       logMsg();
+      log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
       close();
       return 1;
 …
     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;
 …
     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;
 …
   if (cStatus) {
     logMsg();
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
     close();
     return 1;
 …
   char *timeCRPIX = 0;
   char *beamCRVAL = 0;
+  char *polCRVAL = 0;
   cFreqAxis   = -1;
 …
     } else if (strncmp(ctype[iaxis], "STOKES", 6) == 0) {
       cStokesAxis = iaxis;
+      polCRVAL = CRVAL[iaxis];
     } else if (strncmp(ctype[iaxis], "RA", 2) == 0) {
 …
         sprintf(cMsg, "DATA array contains a TIME axis of length %ld.",
           cNAxisTime);
+        logMsg(cMsg);
+        //logMsg(cMsg);
+        log(LogOrigin( className, methodName, WHERE ), LogIO::NORMAL, cMsg);
+      }
 …
+  }
   // 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;
 …
   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;
 …
   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.
 …
   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;
 …
                       &beamNul, beamCol, &anynul, &cStatus)) {
       delete [] beamCol;
       logMsg();
+      log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
       close();
       return 1;
 …
       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;
 …
                       &IFNul, IFCol, &anynul, &cStatus)) {
       delete [] IFCol;
       logMsg();
+      log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
       close();
       return 1;
 …
       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;
 …
             // Variable dimension array.
             if (readDim(DATA, irow+1, &cNAxes, cNAxis)) {
               logMsg();
+              log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
               close();
               return 1;
 …
           if (readDim(XPOLDATA, irow+1, &nAxis, nAxes)) {
             logMsg();
+            log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE );
             close();
             return 1;
 …
     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 ;
+  }
 …
         double &bandwidth)
+{
+  const string methodName = "getHeader()" ;
   // Has the file been opened?
   if (!cSDptr) {
 …
     // 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).
 …
           // 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 {
 …
+      }
+    }
     // Translate to FITS standard names.
     if (strncmp(dopplerFrame, "TOP", 3) == 0) {
 …
     } 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;
+  }
 …
   if (cStatus) {
     logMsg();
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
     return 1;
+  }
 …
         double* &endFreq)
+{
+  const string methodName = "getFreqInfo()" ;
   float  fqRefPix;
   double fqDelt, fqRefVal;
 …
                       &IFNul, IFCol, &anynul, &cStatus)) {
       delete [] IFCol;
       logMsg();
+      log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
       close();
       return 1;
 …
         double* &positions)
+{
+  const string methodName = "findRange()" ;
   // Has the file been opened?
   if (!cSDptr) {
 …
       delete [] beamCol;
       delete [] sel;
       logMsg();
+      log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
       return 1;
+    }
 …
       delete [] IFCol;
       delete [] sel;
       logMsg();
+      log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
       return 1;
+    }
 …
     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;
 …
         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;
 …
           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;
 …
     nSel = 0;
     delete [] positions;
     logMsg();
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
     cStatus = 0;
     return 1;
 …
         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) {
 …
               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;
+      }
 …
+  }
+  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);
 …
   int nPol = cNPol[iIF];
+  if ( cData[STOKES].colnum > 0 )
+    nPol = 1 ;
   if (cGetSpectra || cGetXPol) {
     int nxpol = cGetXPol ? 2*nChan : 0;
 …
   mbrec.nChan[0] = nChan;
   mbrec.nPol[0]  = nPol;
+  mbrec.polNo = iPol ;
   readData(FqRefPix, cRow, mbrec.fqRefPix);
 …
   if (cStatus) {
     logMsg();
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
     return 1;
+  }
 …
   if (cStatus) {
     logMsg();
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
     return 1;
+  }
 …
         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.");
+        }
 …
           blc, trc, inc, 0, mbrec.spectra[0] + iPol*nChan, &anynul,
           &cStatus)) {
         logMsg();
+        log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
         delete [] blc;
         delete [] trc;
 …
             cNAxis, blc, trc, inc, 0, mbrec.flagged[0] + iPol*nChan, &anynul,
             &cStatus)) {
           logMsg();
+          log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
           delete [] blc;
           delete [] trc;
 …
     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;
 …
   if (cStatus) {
     logMsg();
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
     return 1;
+  }
 …
   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;
+  }
 …
+}
 //------------------------------------------------------- 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 ;
+}
 …
   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;
 …
         short iPol)
+{
+  const string methodName = "alfaCal()" ;
   int  calOn;
   char chars[32];
 …
   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;

branches/mergetest/external/atnf/PKSIO/SDFITSreader.h

-              r1720
+              r1779
 #include <atnf/PKSIO/MBrecord.h>
+#include <casa/Logging/LogIO.h>
 #include <fitsio.h>
 using namespace std;
+using namespace casa;
 // <summary>
 …
     int  cBeam_1rel, cIF_1rel;
+    // for GBT
+    int *cPols ;
     enum {SCAN, CYCLE, DATE_OBS, TIME, EXPOSURE, OBJECT, OBJ_RA, OBJ_DEC,
           RESTFRQ, OBSMODE, BEAM, IF, FqRefVal, FqDelt, FqRefPix, RA, DEC,
 …
           BASELIN, BASESUB, DATA, FLAGGED, DATAXED, XPOLDATA, REFBEAM, TCAL,
           TCALTIME, AZIMUTH, ELEVATIO, PARANGLE, FOCUSAXI, FOCUSTAN, FOCUSROT,
+          TAMBIENT, PRESSURE, HUMIDITY, WINDSPEE, WINDDIRE, NDATA};
+          TAMBIENT, PRESSURE, HUMIDITY, WINDSPEE, WINDDIRE, STOKES, SIG, CAL,
+          VFRAME, RVSYS, VELDEF, NDATA};
     // Message handling.
     virtual void logMsg(const char *msg = 0x0);
+    void log(LogOrigin origin, LogIO::Command cmd, const char *msg = 0x0);
     void findData(int iData, char *name, int type);
 …
     // These are for GBT data.
     int   cGBT, cFirstScanNo;
+    double cGLastUTC[4] ;
+    int cGLastScan[4] ;
+    int cGCycleNo[4] ;
 };

branches/mergetest/external/atnf/PKSIO/SDFITSwriter.cc

-              r1736
+              r1779
 #include <atnf/PKSIO/MBrecord.h>
-#include <atnf/PKSIO/PKSmsg.h>
 #include <atnf/PKSIO/SDFITSwriter.h>
+#include <casa/Logging/LogIO.h>
 #include <casa/iostream.h>
 …
 // Factor to convert radians to degrees.
 const double R2D = 180.0 / PI;
+// Class name
+const string className = "SDFITSwriter" ;
 //------------------------------------------------- SDFITSwriter::SDFITSwriter
 …
   // Default constructor.
   cSDptr = 0x0;
-  // By default, messages are written to stderr.
-  initMsg();
+}
 …
         int    extraSysCal)
+{
+  const string methodName = "create()" ;
   if (cSDptr) {
     logMsg("ERROR: Output file already open, close it first.");
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "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.
 …
   cStatus = 0;
   if (fits_create_file(&cSDptr, sdname, &cStatus)) {
     sprintf(cMsg, "ERROR: Failed to create SDFITS file\n       %s", sdName);
     logMsg(cMsg);
+    sprintf(cMsg, "Failed to create SDFITS file\n       %s", sdName);
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, cMsg);
     return cStatus;
+  }
 …
   // Write required primary header keywords.
   if (fits_write_imghdr(cSDptr, 8, 0, 0, &cStatus)) {
     logMsg("ERROR: Failed to write required primary header keywords.");
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Failed to write required primary header keywords.");
     return cStatus;
+  }
 …
   if (cStatus) {
     logMsg("ERROR: Failed in writing primary header.");
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Failed in writing primary header.");
     return cStatus;
+  }
 …
   if (fits_create_tbl(cSDptr, BINARY_TBL, nrow, ncol, NULL, NULL, NULL,
       "SINGLE DISH", &cStatus)) {
     logMsg("ERROR: Failed to create a binary table extension.");
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Failed to create a binary table extension.");
     return 1;
+  }
 …
   if (cStatus) {
     logMsg("ERROR: Failed in writing binary table header.");
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Failed in writing binary table header.");
+  }
 …
 int SDFITSwriter::write(MBrecord &mbrec)
+{
+  const string methodName = "write()" ;
+  LogIO os( LogOrigin( className, methodName, WHERE ) ) ;
   char *cptr;
 …
   int IFno = mbrec.IFno[0];
   if (IFno < 1 || cNIF < IFno) {
+    cerr << "SDFITSwriter::write: "
+         << "Invalid IF number " << IFno
+         << " (maximum " << cNIF << ")." << endl;
+    os << LogIO::WARN
+       << "SDFITSwriter::write: "
+       << "Invalid IF number " << IFno
+       << " (maximum " << cNIF << ")." << LogIO::POST ;
     return 1;
+  }
 …
   int nChan = cNChan[iIF];
   if (mbrec.nChan[0] != nChan) {
+    cerr << "SDFITSriter::write: "
+         << "Wrong number of channels for IF " << IFno << "," << endl
+         << "                    "
+         << "got " << nChan << " should be " << mbrec.nChan[0] << "." << endl;
+    os << LogIO::WARN
+       << "SDFITSriter::write: "
+       << "Wrong number of channels for IF " << IFno << "," << endl
+       << "                    "
+       << "got " << nChan << " should be " << mbrec.nChan[0] << "." << endl;
+    os << LogIO::POST ;
     return 1;
+  }
 …
   int nPol = cNPol[iIF];
   if (mbrec.nPol[0] != nPol) {
+    cerr << "SDFITSriter::write: "
+         << "Wrong number of polarizations for IF " << IFno << "," << endl
+         << "                    "
+         << "got " << nPol << " should be " << mbrec.nPol[0] << "." << endl;
+    os << LogIO::WARN
+       << "SDFITSriter::write: "
+       << "Wrong number of polarizations for IF " << IFno << "," << endl
+       << "                    "
+       << "got " << nPol << " should be " << mbrec.nPol[0] << "." << endl;
+    os << LogIO::POST ;
     return 1;
+  }
 …
   if (cStatus) {
     logMsg("ERROR: Failed in writing binary table entry.");
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Failed in writing binary table entry.");
+  }
 …
+{
+  const string methodName = "history()" ;
   if (!cSDptr) {
     return 1;
 …
   if (fits_write_history(cSDptr, text, &cStatus)) {
     logMsg("ERROR: Failed in writing HISTORY records.");
+    log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Failed in writing HISTORY records.");
+  }
 …
 void SDFITSwriter::close()
+{
+  const string methodName = "close()" ;
   if (cSDptr) {
     cStatus = 0;
     if (fits_close_file(cSDptr, &cStatus)) {
       logMsg("ERROR: Failed to close file.");
+      log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Failed to close file.");
+    }
 …
 void SDFITSwriter::deleteFile()
+{
+  const string methodName = "deleteFile()" ;
   if (cSDptr) {
     cStatus = 0;
     if (fits_delete_file(cSDptr, &cStatus)) {
       logMsg("ERROR: Failed to close and delete file.");
+      log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Failed to close and delete file.");
+    }
 …
+}
 //------------------------------------------------------- SDFITSwriter::logMsg
+//------------------------------------------------------- SDFITSwriter::log
 // 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)
+void SDFITSwriter::log(LogOrigin origin, LogIO::Command cmd, const char *msg)
+{
+  PKSmsg::logMsg(msg);
+  LogIO os( origin ) ;
+  os << cmd << msg << endl ;
   if (cStatus) {
     fits_get_errstatus(cStatus, cMsg);
     PKSmsg::logMsg(cMsg);
+    os << cMsg << endl ;
     while (fits_read_errmsg(cMsg)) {
+      PKSmsg::logMsg(cMsg);
+    }
+  }
+      os << cMsg << endl ;
+    }
+  }
+  os << LogIO::POST ;
+}

branches/mergetest/external/atnf/PKSIO/SDFITSwriter.h

-              r1720
+              r1779
 #include <atnf/PKSIO/MBrecord.h>
 #include <atnf/PKSIO/PKSmsg.h>
+#include <casa/Logging/LogIO.h>
 #include <fitsio.h>
 using namespace std;
+using namespace casa;
 // <summary>
 …
 // </summary>
 class SDFITSwriter : public PKSmsg
+class SDFITSwriter
+{
   public:
 …
     // Message handling.
     char cMsg[256];
     virtual void logMsg(const char *msg = 0x0);
+    void log(LogOrigin origin, LogIO::Command cmd, const char *msg = 0x0);
 };

branches/mergetest/external/atnf/PKSIO/makefile

-              r1452
+              r1779
 # $Id: makefile,v 19.0 2003-07-16 03:34:05 aips2adm Exp $
+# $Id$
 XLIBLIST := CFITSIO RPFITS
 …
 # Use the generic AIPS++ class implementation makefile.
 #------------------------------------------------------
 include $(word 1, $(AIPSPATH))/code/install/makefile.imp
+include $(word 1, $(CASAPATH))/code/install/makefile.imp

Context Navigation

Legend:

branches/mergetest

branches/mergetest/external/atnf/PKSIO/FITSreader.cc

branches/mergetest/external/atnf/PKSIO/FITSreader.h

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

branches/mergetest/external/atnf/PKSIO/MBrecord.cc

branches/mergetest/external/atnf/PKSIO/MBrecord.h

branches/mergetest/external/atnf/PKSIO/PKSFITSreader.cc

branches/mergetest/external/atnf/PKSIO/PKSFITSreader.h

branches/mergetest/external/atnf/PKSIO/PKSMS2reader.cc

branches/mergetest/external/atnf/PKSIO/PKSMS2reader.h

branches/mergetest/external/atnf/PKSIO/PKSMS2writer.cc

branches/mergetest/external/atnf/PKSIO/PKSMS2writer.h

branches/mergetest/external/atnf/PKSIO/PKSSDwriter.cc

branches/mergetest/external/atnf/PKSIO/PKSSDwriter.h

branches/mergetest/external/atnf/PKSIO/PKSreader.cc

branches/mergetest/external/atnf/PKSIO/PKSreader.h

branches/mergetest/external/atnf/PKSIO/PKSrecord.h

branches/mergetest/external/atnf/PKSIO/PKSwriter.h

branches/mergetest/external/atnf/PKSIO/SDFITSreader.cc

branches/mergetest/external/atnf/PKSIO/SDFITSreader.h

branches/mergetest/external/atnf/PKSIO/SDFITSwriter.cc

branches/mergetest/external/atnf/PKSIO/SDFITSwriter.h

branches/mergetest/external/atnf/PKSIO/makefile

Download in other formats: