Changes in trunk/src/STMath.cpp [2345:2467]

File:

• trunk/src/STMath.cpp (modified) (30 diffs)

trunk/src/STMath.cpp

@@ -57 +57 @@
 using namespace asap;
 
+// 2012/02/17 TN
+// Since STGrid is implemented, average doesn't consider direction 
+// when accumulating
 // tolerance for direction comparison (rad)
-#define TOL_OTF    1.0e-15
-#define TOL_POINT  2.9088821e-4  // 1 arcmin
+// #define TOL_OTF    1.0e-15
+// #define TOL_POINT  2.9088821e-4  // 1 arcmin
 
 STMath::STMath(bool insitu) :
@@ -83 +86 @@
   WeightType wtype = stringToWeight(weight);
 
+  // 2012/02/17 TN
+  // Since STGrid is implemented, average doesn't consider direction 
+  // when accumulating
   // check if OTF observation
-  String obstype = in[0]->getHeader().obstype ;
-  Double tol = 0.0 ;
-  if ( (obstype.find( "OTF" ) != String::npos) || (obstype.find( "OBSERVE_TARGET" ) != String::npos) ) {
-    tol = TOL_OTF ;
-  }
-  else {
-    tol = TOL_POINT ;
-  }
+//   String obstype = in[0]->getHeader().obstype ;
+//   Double tol = 0.0 ;
+//   if ( (obstype.find( "OTF" ) != String::npos) || (obstype.find( "OBSERVE_TARGET" ) != String::npos) ) {
+//     tol = TOL_OTF ;
+//   }
+//   else {
+//     tol = TOL_POINT ;
+//   }
 
   // output
@@ -142 +148 @@
   while (!iter.pastEnd()) {
     Table subt = iter.table();
-//     // copy the first row of this selection into the new table
-//     tout.addRow();
-//     TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
-//     // re-index to 0
-//     if ( avmode != "SCAN" && avmode != "SOURCE" ) {
-//       scanColOut.put(outrowCount, uInt(0));
+    // copy the first row of this selection into the new table
+    tout.addRow();
+    TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
+    // re-index to 0
+    if ( avmode != "SCAN" && avmode != "SOURCE" ) {
+      scanColOut.put(outrowCount, uInt(0));
+    }
+    ++outrowCount;
+    // 2012/02/17 TN
+    // Since STGrid is implemented, average doesn't consider direction 
+    // when accumulating
+//     MDirection::ScalarColumn dircol ;
+//     dircol.attach( subt, "DIRECTION" ) ;
+//     Int length = subt.nrow() ;
+//     vector< Vector<Double> > dirs ;
+//     vector<int> indexes ;
+//     for ( Int i = 0 ; i < length ; i++ ) {
+//       Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
+//       //os << << count++ << ": " ;
+//       //os << "[" << t[0] << "," << t[1] << "]" << LogIO::POST ;
+//       bool adddir = true ;
+//       for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
+//         //if ( allTrue( t == dirs[j] ) ) {
+//         Double dx = t[0] - dirs[j][0] ;
+//         Double dy = t[1] - dirs[j][1] ;
+//         Double dd = sqrt( dx * dx + dy * dy ) ;
+//         //if ( allNearAbs( t, dirs[j], tol ) ) {
+//         if ( dd <= tol ) {
+//           adddir = false ;
+//           break ;
+//         }
+//       }
+//       if ( adddir ) {
+//         dirs.push_back( t ) ;
+//         indexes.push_back( i ) ;
+//       }
 //     }
-//     ++outrowCount;
-    MDirection::ScalarColumn dircol ;
-    dircol.attach( subt, "DIRECTION" ) ;
-    Int length = subt.nrow() ;
-    vector< Vector<Double> > dirs ;
-    vector<int> indexes ;
-    for ( Int i = 0 ; i < length ; i++ ) {
-      Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
-      //os << << count++ << ": " ;
-      //os << "[" << t[0] << "," << t[1] << "]" << LogIO::POST ;
-      bool adddir = true ;
-      for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
-        //if ( allTrue( t == dirs[j] ) ) {
-        Double dx = t[0] - dirs[j][0] ;
-        Double dy = t[1] - dirs[j][1] ;
-        Double dd = sqrt( dx * dx + dy * dy ) ;
-        //if ( allNearAbs( t, dirs[j], tol ) ) {
-        if ( dd <= tol ) {
-          adddir = false ;
-          break ;
-        }
-      }
-      if ( adddir ) {
-        dirs.push_back( t ) ;
-        indexes.push_back( i ) ;
-      }
-    }
-    uInt rowNum = dirs.size() ;
-    tout.addRow( rowNum ) ;
-    for ( uInt i = 0 ; i < rowNum ; i++ ) {
-      TableCopy::copyRows( tout, subt, outrowCount+i, indexes[i], 1 ) ;
-      // re-index to 0
-      if ( avmode != "SCAN" && avmode != "SOURCE" ) {
-        scanColOut.put(outrowCount+i, uInt(0));
-      }        
-    }
-    outrowCount += rowNum ;
+//     uInt rowNum = dirs.size() ;
+//     tout.addRow( rowNum ) ;
+//     for ( uInt i = 0 ; i < rowNum ; i++ ) {
+//       TableCopy::copyRows( tout, subt, outrowCount+i, indexes[i], 1 ) ;
+//       // re-index to 0
+//       if ( avmode != "SCAN" && avmode != "SOURCE" ) {
+//         scanColOut.put(outrowCount+i, uInt(0));
+//       }        
+//     }
+//     outrowCount += rowNum ;
     ++iter;
   }
@@ -218 +227 @@
       }
 
-      vector<uInt> removeRows ;
-      uInt nrsubt = subt.nrow() ;
-      for ( uInt irow = 0 ; irow < nrsubt ; irow++ ) {
-        //if ( !allTrue((subt.col("DIRECTION").getArrayDouble(TableExprId(irow)))==rec.asArrayDouble("DIRECTION")) ) {
-        Vector<Double> x0 = (subt.col("DIRECTION").getArrayDouble(TableExprId(irow))) ;
-        Vector<Double> x1 = rec.asArrayDouble("DIRECTION") ;
-        double dx = x0[0] - x1[0];
-        double dy = x0[1] - x1[1];
-        Double dd = sqrt( dx * dx + dy * dy ) ;
-        //if ( !allNearAbs((subt.col("DIRECTION").getArrayDouble(TableExprId(irow))), rec.asArrayDouble("DIRECTION"), tol ) ) {
-        if ( dd > tol ) {
-          removeRows.push_back( irow ) ;
-        }
-      }
-      if ( removeRows.size() != 0 ) {
-        subt.removeRow( removeRows ) ;
-      }
+      // 2012/02/17 TN
+      // Since STGrid is implemented, average doesn't consider direction 
+      // when accumulating
+//       vector<uInt> removeRows ;
+//       uInt nrsubt = subt.nrow() ;
+//       for ( uInt irow = 0 ; irow < nrsubt ; irow++ ) {
+//         //if ( !allTrue((subt.col("DIRECTION").getArrayDouble(TableExprId(irow)))==rec.asArrayDouble("DIRECTION")) ) {
+//         Vector<Double> x0 = (subt.col("DIRECTION").getArrayDouble(TableExprId(irow))) ;
+//         Vector<Double> x1 = rec.asArrayDouble("DIRECTION") ;
+//         double dx = x0[0] - x1[0];
+//         double dy = x0[1] - x1[1];
+//         Double dd = sqrt( dx * dx + dy * dy ) ;
+//         //if ( !allNearAbs((subt.col("DIRECTION").getArrayDouble(TableExprId(irow))), rec.asArrayDouble("DIRECTION"), tol ) ) {
+//         if ( dd > tol ) {
+//           removeRows.push_back( irow ) ;
+//         }
+//       }
+//       if ( removeRows.size() != 0 ) {
+//         subt.removeRow( removeRows ) ;
+//       }
       
-      if ( nrsubt == removeRows.size() )
-        throw(AipsError("Averaging data is empty.")) ;
+//       if ( nrsubt == removeRows.size() )
+//         throw(AipsError("Averaging data is empty.")) ;
 
       specCol.attach(subt,"SPECTRA");
@@ -321 +333 @@
 {
   (void) mode; // currently unused
+  // 2012/02/17 TN
+  // Since STGrid is implemented, average doesn't consider direction 
+  // when accumulating
   // check if OTF observation
-  String obstype = in->getHeader().obstype ;
-  Double tol = 0.0 ;
-  if ( obstype.find( "OTF" ) != String::npos ) {
-    tol = TOL_OTF ;
-  }
-  else {
-    tol = TOL_POINT ;
-  }
+//   String obstype = in->getHeader().obstype ;
+//   Double tol = 0.0 ;
+//   if ( obstype.find( "OTF" ) != String::npos ) {
+//     tol = TOL_OTF ;
+//   }
+//   else {
+//     tol = TOL_POINT ;
+//   }
 
   // clone as this is non insitu
@@ -362 +377 @@
     flagCol.attach(subt,"FLAGTRA");
     tsysCol.attach(subt,"TSYS");
-//     tout.addRow();
-//     TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
-//     if ( avmode != "SCAN") {
-//       scanColOut.put(outrowCount, uInt(0));
+
+    tout.addRow();
+    TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
+    if ( avmode != "SCAN") {
+      scanColOut.put(outrowCount, uInt(0));
+    }
+    Vector<Float> tmp;
+    specCol.get(0, tmp);
+    uInt nchan = tmp.nelements();
+    // have to do channel by channel here as MaskedArrMath
+    // doesn't have partialMedians
+    Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
+    Vector<Float> outspec(nchan);
+    Vector<uChar> outflag(nchan,0);
+    Vector<Float> outtsys(1);/// @fixme when tsys is channel based
+    for (uInt i=0; i<nchan; ++i) {
+      Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
+      MaskedArray<Float> ma = maskedArray(specs,flags);
+      outspec[i] = median(ma);
+      if ( allEQ(ma.getMask(), False) )
+        outflag[i] = userflag;// flag data
+    }
+    outtsys[0] = median(tsysCol.getColumn());
+    specColOut.put(outrowCount, outspec);
+    flagColOut.put(outrowCount, outflag);
+    tsysColOut.put(outrowCount, outtsys);
+    Double intsum = sum(intCol.getColumn());
+    intColOut.put(outrowCount, intsum);
+    ++outrowCount;
+    ++iter;
+
+    // 2012/02/17 TN
+    // Since STGrid is implemented, average doesn't consider direction 
+    // when accumulating
+//     MDirection::ScalarColumn dircol ;
+//     dircol.attach( subt, "DIRECTION" ) ;
+//     Int length = subt.nrow() ;
+//     vector< Vector<Double> > dirs ;
+//     vector<int> indexes ;
+//     // Handle MX mode averaging
+//     if (in->nbeam() > 1 ) {      
+//       length = 1;
 //     }
-//     Vector<Float> tmp;
-//     specCol.get(0, tmp);
-//     uInt nchan = tmp.nelements();
-//     // have to do channel by channel here as MaskedArrMath
-//     // doesn't have partialMedians
-//     Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
-//     Vector<Float> outspec(nchan);
-//     Vector<uChar> outflag(nchan,0);
-//     Vector<Float> outtsys(1);/// @fixme when tsys is channel based
-//     for (uInt i=0; i<nchan; ++i) {
-//       Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
-//       MaskedArray<Float> ma = maskedArray(specs,flags);
-//       outspec[i] = median(ma);
-//       if ( allEQ(ma.getMask(), False) )
-//         outflag[i] = userflag;// flag data
+//     for ( Int i = 0 ; i < length ; i++ ) {
+//       Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
+//       bool adddir = true ;
+//       for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
+//         //if ( allTrue( t == dirs[j] ) ) {
+//         Double dx = t[0] - dirs[j][0] ;
+//         Double dy = t[1] - dirs[j][1] ;
+//         Double dd = sqrt( dx * dx + dy * dy ) ;
+//         //if ( allNearAbs( t, dirs[j], tol ) ) {
+//         if ( dd <= tol ) {
+//           adddir = false ;
+//           break ;
+//         }
+//       }
+//       if ( adddir ) {
+//         dirs.push_back( t ) ;
+//         indexes.push_back( i ) ;
+//       }
 //     }
-//     outtsys[0] = median(tsysCol.getColumn());
-//     specColOut.put(outrowCount, outspec);
-//     flagColOut.put(outrowCount, outflag);
-//     tsysColOut.put(outrowCount, outtsys);
-//     Double intsum = sum(intCol.getColumn());
-//     intColOut.put(outrowCount, intsum);
-//     ++outrowCount;
+//     uInt rowNum = dirs.size() ;
+//     tout.addRow( rowNum );
+//     for ( uInt i = 0 ; i < rowNum ; i++ ) {
+//       TableCopy::copyRows(tout, subt, outrowCount+i, indexes[i], 1) ;
+//       // Handle MX mode averaging
+//       if ( avmode != "SCAN") {
+//         scanColOut.put(outrowCount+i, uInt(0));
+//       }
+//     }
+//     MDirection::ScalarColumn dircolOut ;
+//     dircolOut.attach( tout, "DIRECTION" ) ;
+//     for ( uInt irow = 0 ; irow < rowNum ; irow++ ) {
+//       Vector<Double> t = \
+//      dircolOut(outrowCount+irow).getAngle(Unit(String("rad"))).getValue() ;
+//       Vector<Float> tmp;
+//       specCol.get(0, tmp);
+//       uInt nchan = tmp.nelements();
+//       // have to do channel by channel here as MaskedArrMath
+//       // doesn't have partialMedians
+//       Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
+//       // mask spectra for different DIRECTION
+//       for ( uInt jrow = 0 ; jrow < subt.nrow() ; jrow++ ) {
+//         Vector<Double> direction = \
+//        dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
+//         //if ( t[0] != direction[0] || t[1] != direction[1] ) {
+//         Double dx = t[0] - direction[0];
+//         Double dy = t[1] - direction[1];
+//         Double dd = sqrt(dx*dx + dy*dy);
+//         //if ( !allNearAbs( t, direction, tol ) ) {
+//         if ( dd > tol &&  in->nbeam() < 2 ) {
+//           flags[jrow] = userflag ;
+//         }
+//       }
+//       Vector<Float> outspec(nchan);
+//       Vector<uChar> outflag(nchan,0);
+//       Vector<Float> outtsys(1);/// @fixme when tsys is channel based
+//       for (uInt i=0; i<nchan; ++i) {
+//         Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
+//         MaskedArray<Float> ma = maskedArray(specs,flags);
+//         outspec[i] = median(ma);
+//         if ( allEQ(ma.getMask(), False) )
+//           outflag[i] = userflag;// flag data
+//       }
+//       outtsys[0] = median(tsysCol.getColumn());
+//       specColOut.put(outrowCount+irow, outspec);
+//       flagColOut.put(outrowCount+irow, outflag);
+//       tsysColOut.put(outrowCount+irow, outtsys);
+//       Vector<Double> integ = intCol.getColumn() ;
+//       MaskedArray<Double> mi = maskedArray( integ, flags ) ;
+//       Double intsum = sum(mi);
+//       intColOut.put(outrowCount+irow, intsum);
+//     }
+//     outrowCount += rowNum ;
 //     ++iter;
-    MDirection::ScalarColumn dircol ;
-    dircol.attach( subt, "DIRECTION" ) ;
-    Int length = subt.nrow() ;
-    vector< Vector<Double> > dirs ;
-    vector<int> indexes ;
-    // Handle MX mode averaging
-    if (in->nbeam() > 1 ) {      
-      length = 1;
-    }
-    for ( Int i = 0 ; i < length ; i++ ) {
-      Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
-      bool adddir = true ;
-      for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
-        //if ( allTrue( t == dirs[j] ) ) {
-        Double dx = t[0] - dirs[j][0] ;
-        Double dy = t[1] - dirs[j][1] ;
-        Double dd = sqrt( dx * dx + dy * dy ) ;
-        //if ( allNearAbs( t, dirs[j], tol ) ) {
-        if ( dd <= tol ) {
-          adddir = false ;
-          break ;
-        }
-      }
-      if ( adddir ) {
-        dirs.push_back( t ) ;
-        indexes.push_back( i ) ;
-      }
-    }
-    uInt rowNum = dirs.size() ;
-    tout.addRow( rowNum );
-    for ( uInt i = 0 ; i < rowNum ; i++ ) {
-      TableCopy::copyRows(tout, subt, outrowCount+i, indexes[i], 1) ;
-      // Handle MX mode averaging
-      if ( avmode != "SCAN") {
-        scanColOut.put(outrowCount+i, uInt(0));
-      }
-    }
-    MDirection::ScalarColumn dircolOut ;
-    dircolOut.attach( tout, "DIRECTION" ) ;
-    for ( uInt irow = 0 ; irow < rowNum ; irow++ ) {
-      Vector<Double> t = \
-        dircolOut(outrowCount+irow).getAngle(Unit(String("rad"))).getValue() ;
-      Vector<Float> tmp;
-      specCol.get(0, tmp);
-      uInt nchan = tmp.nelements();
-      // have to do channel by channel here as MaskedArrMath
-      // doesn't have partialMedians
-      Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
-      // mask spectra for different DIRECTION
-      for ( uInt jrow = 0 ; jrow < subt.nrow() ; jrow++ ) {
-        Vector<Double> direction = \
-          dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
-        //if ( t[0] != direction[0] || t[1] != direction[1] ) {
-        Double dx = t[0] - direction[0];
-        Double dy = t[1] - direction[1];
-        Double dd = sqrt(dx*dx + dy*dy);
-        //if ( !allNearAbs( t, direction, tol ) ) {
-        if ( dd > tol &&  in->nbeam() < 2 ) {
-          flags[jrow] = userflag ;
-        }
-      }
-      Vector<Float> outspec(nchan);
-      Vector<uChar> outflag(nchan,0);
-      Vector<Float> outtsys(1);/// @fixme when tsys is channel based
-      for (uInt i=0; i<nchan; ++i) {
-        Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
-        MaskedArray<Float> ma = maskedArray(specs,flags);
-        outspec[i] = median(ma);
-        if ( allEQ(ma.getMask(), False) )
-          outflag[i] = userflag;// flag data
-      }
-      outtsys[0] = median(tsysCol.getColumn());
-      specColOut.put(outrowCount+irow, outspec);
-      flagColOut.put(outrowCount+irow, outflag);
-      tsysColOut.put(outrowCount+irow, outtsys);
-      Vector<Double> integ = intCol.getColumn() ;
-      MaskedArray<Double> mi = maskedArray( integ, flags ) ;
-      Double intsum = sum(mi);
-      intColOut.put(outrowCount+irow, intsum);
-    }
-    outrowCount += rowNum ;
-    ++iter;
   }
   return out;
@@ -1910 +1930 @@
   ArrayColumn<Float> specCol(tout, "SPECTRA");
   ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
+  ArrayColumn<Float> tsysCol(tout, "TSYS");
+
   for (uInt i=0; i < tout.nrow(); ++i ) {
     MaskedArray<Float> main  = maskedArray(specCol(i), flagCol(i));
     MaskedArray<Float> maout;
     LatticeUtilities::bin(maout, main, 0, Int(width));
-    /// @todo implement channel based tsys binning
     specCol.put(i, maout.getArray());
     flagCol.put(i, flagsFromMA(maout));
+    if (tsysCol(i).nelements() == specCol(i).nelements()) {
+      MaskedArray<Float> matsysin = maskedArray(tsysCol(i), flagCol(i));
+      MaskedArray<Float> matsysout;
+      LatticeUtilities::bin(matsysout, matsysin, 0, Int(width));
+      tsysCol.put(i, matsysout.getArray());
+    }
     // take only the first binned spectrum's length for the deprecated
     // global header item nChan
@@ -2750 +2777 @@
                              mask, timeCol(i), !first,
                              interp, False);
-          (void) ok; // unused stop compiler nagging      
+          (void) ok; // unused stop compiler nagging
           // back into scantable
           flagOut.resize(maskOut.nelements());
@@ -2954 +2981 @@
                     "Use merge first."));
   
+  // 2012/02/17 TN
+  // Since STGrid is implemented, average doesn't consider direction 
+  // when accumulating
   // check if OTF observation
-  String obstype = in[0]->getHeader().obstype ;
-  Double tol = 0.0 ;
-  if ( obstype.find( "OTF" ) != String::npos ) {
-    tol = TOL_OTF ;
-  }
-  else {
-    tol = TOL_POINT ;
-  }
+//   String obstype = in[0]->getHeader().obstype ;
+//   Double tol = 0.0 ;
+//   if ( obstype.find( "OTF" ) != String::npos ) {
+//     tol = TOL_OTF ;
+//   }
+//   else {
+//     tol = TOL_POINT ;
+//   }
 
   CountedPtr<Scantable> out ;     // processed result 
@@ -3015 +3045 @@
       uInt rows = tmpin[itable]->nrow() ;
       uInt freqnrows = tmpin[itable]->frequencies().table().nrow() ;
+      freqIDCol.attach( tmpin[itable]->table(), "FREQ_ID" ) ;
+      //ifnoCol.attach( tmpin[itable]->table(), "IFNO" ) ;
       for ( uInt irow = 0 ; irow < rows ; irow++ ) {
         if ( freqid[itable].size() == freqnrows ) {
@@ -3020 +3052 @@
         }
         else {
-          freqIDCol.attach( tmpin[itable]->table(), "FREQ_ID" ) ;
-          ifnoCol.attach( tmpin[itable]->table(), "IFNO" ) ;
           uInt id = freqIDCol( irow ) ;
           if ( freqid[itable].size() == 0 || count( freqid[itable].begin(), freqid[itable].end(), id ) == 0 ) {
             //os << "itable = " << itable << ": IF " << id << " is included in the list" << LogIO::POST ;
             vector<double> abcissa = tmpin[itable]->getAbcissa( irow ) ;
+            double lbedge, rbedge;
             freqid[itable].push_back( id ) ;
-            iffreq[itable].push_back( abcissa[0] - 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
-            iffreq[itable].push_back( abcissa[abcissa.size()-1] + 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
+            lbedge = abcissa[0] - 0.5 * ( abcissa[1] - abcissa[0] );
+            rbedge = abcissa[abcissa.size()-1] + 0.5 * ( abcissa[1] - abcissa[0] );
+            iffreq[itable].push_back( min(lbedge, rbedge) ) ;
+            iffreq[itable].push_back( max(lbedge, rbedge) ) ;
           }
         }
@@ -3113 +3146 @@
     // Grouping continuous IF groups (without frequency gap)
     // freqgrp: list of IF group indexes in each frequency group
-    // freqrange: list of minimum and maximum frequency in each frequency group
     // freqgrp[numgrp][nummember]
     // freqgrp: [[ifgrp00, ifgrp01, ifgrp02, ...],
@@ -3119 +3151 @@
     //           ...
     //           [ifgrpn0, ifgrpn1, ifgrpn2, ...]]
-    // freqrange[numgrp*2]
-    // freqrange: [min_grp0, max_grp0, min_grp1, max_grp1, ...]
     vector< vector<uInt> > freqgrp ;
     double freqrange = 0.0 ;
@@ -3136 +3166 @@
       freqrange = ifgfreq[2*i+1] ;
     }
-        
+
 
     // print IF groups
@@ -3325 +3355 @@
 
     // save column values in the vector
-    vector< vector<uInt> > freqTableIdVec( insize ) ;
     vector< vector<uInt> > freqIdVec( insize ) ;
     vector< vector<uInt> > ifNoVec( insize ) ;
     for ( uInt itable = 0 ; itable < insize ; itable++ ) {
-      ScalarColumn<uInt> freqIDs ;
-      freqIDs.attach( newin[itable]->frequencies().table(), "ID" ) ;
       ifnoCol.attach( newin[itable]->table(), "IFNO" ) ;
       freqIDCol.attach( newin[itable]->table(), "FREQ_ID" ) ;
-      for ( uInt irow = 0 ; irow < newin[itable]->frequencies().table().nrow() ; irow++ ) {
-        freqTableIdVec[itable].push_back( freqIDs( irow ) ) ;
-      }
       for ( uInt irow = 0 ; irow < newin[itable]->table().nrow() ; irow++ ) {
         freqIdVec[itable].push_back( freqIDCol( irow ) ) ;
@@ -3350 +3374 @@
       vector<uInt> freqIdUpdate ; 
       for ( uInt irow = 0 ; irow < rows ; irow++ ) {
-        uInt ifno = ifNoVec[itable][irow] ;  // IFNO is reset by group index
+        uInt ifno = ifNoVec[itable][irow] ;  // IFNO is reset by group index (IF group index)
         double minfreq = ifgfreq[2*ifno] ;
         double maxfreq = ifgfreq[2*ifno+1] ;
@@ -3358 +3382 @@
         double resol = abcissa[1] - abcissa[0] ;
         //os << "abcissa range  : [" << abcissa[0] << "," << abcissa[nchan-1] << "]" << LogIO::POST ;
-        if ( minfreq <= abcissa[0] ) 
-          nminchan = 0 ;
+        uInt imin, imax;
+        int sigres;
+        if ( resol >= 0. ) {
+          imin = 0;
+          imax = nchan - 1 ;
+          sigres = 1 ;
+        } else {
+          imin = nchan - 1 ;
+          imax = 0 ;
+          sigres = -1 ;
+          resol = abs(resol) ;
+        }
+        if ( minfreq <= abcissa[imin] ) 
+          nminchan = imin ;
         else {
           //double cfreq = ( minfreq - abcissa[0] ) / resol ;
-          double cfreq = ( minfreq - abcissa[0] + 0.5 * resol ) / resol ;
-          nminchan = int(cfreq) + ( ( cfreq - int(cfreq) <= 0.5 ) ? 0 : 1 ) ;
+          double cfreq = ( minfreq - abcissa[imin] + 0.5 * resol ) / resol ;
+          nminchan = imin + sigres * (int(cfreq) + ( ( cfreq - int(cfreq) <= 0.5 ) ? 0 : 1 )) ;
         }
-        if ( maxfreq >= abcissa[abcissa.size()-1] )
-          nmaxchan = abcissa.size() - 1 ;
+        if ( maxfreq >= abcissa[imax] )
+          nmaxchan = imax;
         else {
           //double cfreq = ( abcissa[abcissa.size()-1] - maxfreq ) / resol ;
-          double cfreq = ( abcissa[abcissa.size()-1] - maxfreq + 0.5 * resol ) / resol ;
-          nmaxchan = abcissa.size() - 1 - int(cfreq) - ( ( cfreq - int(cfreq) >= 0.5 ) ? 1 : 0 ) ;
+          double cfreq = ( abcissa[imax] - maxfreq + 0.5 * resol ) / resol ;
+          nmaxchan = imax - sigres * (int(cfreq) +( ( cfreq - int(cfreq) >= 0.5 ) ? 1 : 0 )) ;
         }
         //os << "channel range (" << irow << "): [" << nminchan << "," << nmaxchan << "]" << LogIO::POST ;
+        if ( nmaxchan < nminchan ) {
+          int tmp = nmaxchan ;
+          nmaxchan = nminchan ;
+          nminchan = tmp ;
+        }
         if ( nmaxchan > nminchan ) {
           newin[itable]->reshapeSpectrum( nminchan, nmaxchan, irow ) ;
           int newchan = nmaxchan - nminchan + 1 ;
-          if ( count( freqIdUpdate.begin(), freqIdUpdate.end(), freqIdVec[itable][irow] ) == 0 && newchan < nchan )
+          if ( count( freqIdUpdate.begin(), freqIdUpdate.end(), freqIdVec[itable][irow] ) == 0 && newchan < nchan ) {
             freqIdUpdate.push_back( freqIdVec[itable][irow] ) ;
+
+            // Update before nminchan is lost
+            uInt freqId = freqIdVec[itable][irow] ;
+            Double refpix ;
+            Double refval ;
+            Double increment ;
+            newin[itable]->frequencies().getEntry( refpix, refval, increment, freqId ) ; 
+            //refval = refval - ( refpix - nminchan ) * increment ;
+            refval = abcissa[nminchan] ;
+            refpix = 0 ;
+            newin[itable]->frequencies().setEntry( refpix, refval, increment, freqId ) ;
+          }
         }
         else {
@@ -3383 +3436 @@
         }
       }
-      for ( uInt i = 0 ; i < freqIdUpdate.size() ; i++ ) {
-        uInt freqId = freqIdUpdate[i] ;
-        Double refpix ;
-        Double refval ;
-        Double increment ;
+//       for ( uInt i = 0 ; i < freqIdUpdate.size() ; i++ ) {
+//      uInt freqId = freqIdUpdate[i] ;
+//      Double refpix ;
+//      Double refval ;
+//      Double increment ;
         
-        // update row
-        newin[itable]->frequencies().getEntry( refpix, refval, increment, freqId ) ; 
-        refval = refval - ( refpix - nminchan ) * increment ;
-        refpix = 0 ;
-        newin[itable]->frequencies().setEntry( refpix, refval, increment, freqId ) ;
-      }    
+//      // update row
+//      newin[itable]->frequencies().getEntry( refpix, refval, increment, freqId ) ; 
+//      refval = refval - ( refpix - nminchan ) * increment ;
+//      refpix = 0 ;
+//      newin[itable]->frequencies().setEntry( refpix, refval, increment, freqId ) ;
+//       }    
     }
 
@@ -3431 +3484 @@
             if ( nchan < minchan ) {
               minchan = nchan ;
-              maxdnu = dnu ;
+              maxdnu = abs(dnu) ;
               newin[tableid]->frequencies().getEntry( refpixref, refvalref, refinc, rowid ) ;
               refreqid = rowid ;
               reftable = tableid ;
+              // Spectra are reversed when dnu < 0
+              if (dnu < 0)
+                refpixref = minchan - 1 -refpixref ;
             }
           }
@@ -3449 +3505 @@
           //os << "   regridChannel applied to " ;
           //if ( tableid != reftable ) 
-          refreqid = newin[tableid]->frequencies().addEntry( refpixref, refvalref, refinc ) ;
+          refreqid = newin[tableid]->frequencies().addEntry( refpixref, refvalref, maxdnu ) ;
           for ( uInt irow = 0 ; irow < newin[tableid]->table().nrow() ; irow++ ) {
             uInt tfreqid = freqIdVec[tableid][irow] ;
@@ -3470 +3526 @@
           vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
           minchan = abcissa.size() ;
-          maxdnu = abcissa[1] - abcissa[0] ;
+          maxdnu = abs( abcissa[1] - abcissa[0]) ;
         }
       }
@@ -3540 +3596 @@
           vector<double> abcissa = tmpout->getAbcissa( irow ) ;
           double bw = ( abcissa[1] - abcissa[0] ) * abcissa.size() ;
-          int nchan = (int)( bw / gmaxdnu[igrp] ) ;
+          int nchan = (int) abs( bw / gmaxdnu[igrp] ) ;
+          // All spectra will have positive frequency increments
           tmpout->regridChannel( nchan, gmaxdnu[igrp], irow ) ;
           break ;
@@ -3557 +3614 @@
     ScalarColumn<uInt> freqidColOut ;
     freqidColOut.attach( out->table(), "FREQ_ID" ) ;
-    MDirection::ScalarColumn dirColOut ;
-    dirColOut.attach( out->table(), "DIRECTION" ) ;
+//     MDirection::ScalarColumn dirColOut ;
+//     dirColOut.attach( out->table(), "DIRECTION" ) ;
     Table &tab = tmpout->table() ;
     Block<String> cols(1);
     cols[0] = String("POLNO") ;
     TableIterator iter( tab, cols ) ;
-    bool done = false ;
     vector< vector<uInt> > sizes( freqgrp.size() ) ;
     while( !iter.pastEnd() ) {
@@ -3575 +3631 @@
       ScalarColumn<uInt> polnos ;
       polnos.attach( iter.table(), "POLNO" ) ;
-      MDirection::ScalarColumn dircol ;
-      dircol.attach( iter.table(), "DIRECTION" ) ;
+//       MDirection::ScalarColumn dircol ;
+//       dircol.attach( iter.table(), "DIRECTION" ) ;
       uInt polno = polnos( 0 ) ;
       //os << "POLNO iteration: " << polno << LogIO::POST ;
@@ -3614 +3670 @@
 //       }
       // get a list of number of channels for each frequency group member
-      if ( !done ) {
-        for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
-          sizes[igrp].resize( freqgrp[igrp].size() ) ;
-          for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
-            for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
-              uInt ifno = ifnoCol( irow ) ;
-              if ( ifno == freqgrp[igrp][imem] ) {
-                Vector<Float> spec = specCols( irow ) ;
-                sizes[igrp][imem] = spec.nelements() ;
-                break ;
-              }                
-            }
-          }
-        }
-        done = true ;
+      for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
+        sizes[igrp].resize( freqgrp[igrp].size() ) ;
+        for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
+          for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
+            uInt ifno = ifnoCol( irow ) ;
+            if ( ifno == freqgrp[igrp][imem] ) {
+              Vector<Float> spec = specCols( irow ) ;
+              sizes[igrp][imem] = spec.nelements() ;
+              break ;
+            }
+          }
+        }
       }
       // combine spectra
@@ -3635 +3688 @@
         if ( polout == polno ) {
           uInt ifout = ifnoColOut( irow ) ;
-          Vector<Double> direction = dirColOut(irow).getAngle(Unit(String("rad"))).getValue() ;
+//           Vector<Double> direction = dirColOut(irow).getAngle(Unit(String("rad"))).getValue() ;
           uInt igrp ;
           for ( uInt jgrp = 0 ; jgrp < freqgrp.size() ; jgrp++ ) {
@@ -3646 +3699 @@
             for ( uInt jrow = 0 ; jrow < iter.table().nrow() ; jrow++ ) {
               uInt ifno = ifnoCol( jrow ) ;
-              Vector<Double> tdir = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
-              //if ( ifno == freqgrp[igrp][imem] && allTrue( tdir == direction  ) ) {
-              Double dx = tdir[0] - direction[0] ;
-              Double dy = tdir[1] - direction[1] ;
-              Double dd = sqrt( dx * dx + dy * dy ) ;
+              // 2012/02/17 TN
+              // Since STGrid is implemented, average doesn't consider direction 
+              // when accumulating
+//               Vector<Double> tdir = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
+//               //if ( ifno == freqgrp[igrp][imem] && allTrue( tdir == direction  ) ) {
+//               Double dx = tdir[0] - direction[0] ;
+//               Double dy = tdir[1] - direction[1] ;
+//               Double dd = sqrt( dx * dx + dy * dy ) ;
               //if ( ifno == freqgrp[igrp][imem] && allNearAbs( tdir, direction, tol ) ) {
-              if ( ifno == freqgrp[igrp][imem] && dd <= tol ) {
+//               if ( ifno == freqgrp[igrp][imem] && dd <= tol ) {
+              if ( ifno == freqgrp[igrp][imem] ) {
                 Vector<Float> spec = specCols( jrow ) ;
                 Vector<uChar> flag = flagCols( jrow ) ;
@@ -3671 +3728 @@
           specColOut.put( irow, newspec ) ;
           flagColOut.put( irow, newflag ) ;
-          // IFNO renumbering
+          // IFNO renumbering (renumbered as frequency group ID)
           ifnoColOut.put( irow, igrp ) ; 
         }
@@ -3682 +3739 @@
       uInt index = 0 ;
       uInt pixShift = 0 ;
+
       while ( freqgrp[igrp][index] != gmemid[igrp] ) {
         pixShift += sizes[igrp][index++] ;
       }
       for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
-        if ( ifnoColOut( irow ) == gmemid[igrp] && !updated[igrp] ) {
+        //if ( ifnoColOut( irow ) == gmemid[igrp] && !updated[igrp] ) {
+        if ( ifnoColOut( irow ) == igrp && !updated[igrp] ) {
           uInt freqidOut = freqidColOut( irow ) ;
           //os << "freqgrp " << igrp << " freqidOut = " << freqidOut << LogIO::POST ;
@@ -3693 +3752 @@
           double increm ;
           out->frequencies().getEntry( refpix, refval, increm, freqidOut ) ;
+          if (increm < 0)
+            refpix = sizes[igrp][index] - 1 - refpix ; // reversed
           refpix += pixShift ;
-          out->frequencies().setEntry( refpix, refval, increm, freqidOut ) ;
+          out->frequencies().setEntry( refpix, refval, gmaxdnu[igrp], freqidOut ) ;
           updated[igrp] = true ;
         }