Changeset 1779 for branches/mergetest/external/atnf/PKSIO/PKSMS2writer.cc

Timestamp:

07/29/10 19:13:46 (14 years ago)

Author:

Kana Sugimoto

Message:

New Development: Yes

JIRA Issue: No (test merging alma branch)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s):

Description:

Location:

branches/mergetest

Files:

• . (modified) (1 prop)
• external/atnf/PKSIO/PKSMS2writer.cc (modified) (28 diffs)

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

@@ -42 +42 @@
 #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>
@@ -52 +53 @@
 #include <tables/Tables/TiledShapeStMan.h>
 
+// Class name
+const string className = "PKSMS2writer" ;
+
 //------------------------------------------------- PKSMS2writer::PKSMS2writer
 
@@ -59 +63 @@
 {
   cPKSMS = 0x0;
-
-  // By default, messages are written to stderr.
-  initMsg();
 }
 
@@ -92 +93 @@
         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.
@@ -108 +109 @@
 
   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",
@@ -116 +136 @@
   // 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);
@@ -136 +158 @@
 
     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);
@@ -152 +176 @@
   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);
 
@@ -367 +391 @@
   } else if (dopplerFrame == "SOURCE") {
     MFrequency::getType(cDopplerFrame, "REST");
+  } else if (dopplerFrame == "LSRK") {
+    MFrequency::getType(cDopplerFrame, "LSRK");
   }
 
@@ -374 +400 @@
   addDopplerEntry();
   addFeedEntry();
-  addObservationEntry(observer, project);
+  //addObservationEntry(observer, project);
+  addObservationEntry(observer, project, telName);
   addProcessorEntry();
 
@@ -384 +411 @@
 // 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)
@@ -415 +488 @@
     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.
@@ -508 +581 @@
   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;
@@ -571 +650 @@
   cSysCal          = MSSysCal();
   cWeather         = MSWeather();
-
   // Release the main table.
-  delete cPKSMS;
-  cPKSMS = 0x0;
+  delete cPKSMS; 
+  cPKSMS=0x0;
 }
 
@@ -589 +667 @@
   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");
@@ -597 +710 @@
   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);
@@ -711 +829 @@
         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;
@@ -741 +878 @@
 Int PKSMS2writer::addObservationEntry(
         const String observer,
-        const String project)
+        const String project,
+        const String antName)
 {
   // Extend the OBSERVATION subtable.
@@ -748 +886 @@
 
   // 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);
@@ -754 +893 @@
   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);
@@ -767 +907 @@
 
 //--------------------------------------------- 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
@@ -778 +921 @@
         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;
 }
@@ -821 +976 @@
   // 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;
@@ -869 +1035 @@
         const Vector<Double> direction,
         const Vector<Double> properMotion,
-        const Double restFreq,
+        //const Double restFreq,
+        const Vector<Double> restFreq,
         const Double radialVelocity)
 {
@@ -903 +1070 @@
 //  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);
@@ -933 +1101 @@
 
   // Data.
-  cSpWindowCols->name().put(n, "L-band");
+  //cSpWindowCols->name().put(n, "L-band");
+  cSpWindowCols->name().put(n, " ");
   cSpWindowCols->refFrequency().put(n, refFreq);
 
@@ -1015 +1184 @@
         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);
@@ -1029 +1214 @@
 
   // Data.
-  cSysCalCols->tcal().put(n, tcal);
+  //cSysCalCols->tcal().put(n, tcal);
+  cSysCalCols->tcal().put(n, inTcal);
   cSysCalCols->tsys().put(n, tsys);