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

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) :
@@ -86 +83 @@
   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
@@ -148 +142 @@
   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));
-    }
-    ++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 ) ;
-//       }
+//     // 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));
 //     }
-//     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 ;
+//     ++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 ;
     ++iter;
   }
@@ -227 +218 @@
       }
 
-      // 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 ) ;
-//       }
+      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");
@@ -333 +321 @@
 {
   (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
@@ -377 +362 @@
     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));
-    }
-    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;
+//     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;
+    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;
-
-    // 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;
-//     }
-//     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;
@@ -1930 +1910 @@
   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
@@ -2777 +2750 @@
                              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());
@@ -2981 +2954 @@
                     "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 
@@ -3045 +3015 @@
       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 ) {
@@ -3052 +3020 @@
         }
         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 ) ;
-            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) ) ;
+            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] ) ) ;
           }
         }
@@ -3146 +3113 @@
     // 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, ...],
@@ -3151 +3119 @@
     //           ...
     //           [ifgrpn0, ifgrpn1, ifgrpn2, ...]]
+    // freqrange[numgrp*2]
+    // freqrange: [min_grp0, max_grp0, min_grp1, max_grp1, ...]
     vector< vector<uInt> > freqgrp ;
     double freqrange = 0.0 ;
@@ -3166 +3136 @@
       freqrange = ifgfreq[2*i+1] ;
     }
-
+        
 
     // print IF groups
@@ -3355 +3325 @@
 
     // 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 ) ) ;
@@ -3374 +3350 @@
       vector<uInt> freqIdUpdate ; 
       for ( uInt irow = 0 ; irow < rows ; irow++ ) {
-        uInt ifno = ifNoVec[itable][irow] ;  // IFNO is reset by group index (IF group index)
+        uInt ifno = ifNoVec[itable][irow] ;  // IFNO is reset by group index
         double minfreq = ifgfreq[2*ifno] ;
         double maxfreq = ifgfreq[2*ifno+1] ;
@@ -3382 +3358 @@
         double resol = abcissa[1] - abcissa[0] ;
         //os << "abcissa range  : [" << abcissa[0] << "," << abcissa[nchan-1] << "]" << LogIO::POST ;
-        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 ;
+        if ( minfreq <= abcissa[0] ) 
+          nminchan = 0 ;
         else {
           //double cfreq = ( minfreq - abcissa[0] ) / resol ;
-          double cfreq = ( minfreq - abcissa[imin] + 0.5 * resol ) / resol ;
-          nminchan = imin + sigres * (int(cfreq) + ( ( cfreq - int(cfreq) <= 0.5 ) ? 0 : 1 )) ;
+          double cfreq = ( minfreq - abcissa[0] + 0.5 * resol ) / resol ;
+          nminchan = int(cfreq) + ( ( cfreq - int(cfreq) <= 0.5 ) ? 0 : 1 ) ;
         }
-        if ( maxfreq >= abcissa[imax] )
-          nmaxchan = imax;
+        if ( maxfreq >= abcissa[abcissa.size()-1] )
+          nmaxchan = abcissa.size() - 1 ;
         else {
           //double cfreq = ( abcissa[abcissa.size()-1] - maxfreq ) / resol ;
-          double cfreq = ( abcissa[imax] - maxfreq + 0.5 * resol ) / resol ;
-          nmaxchan = imax - sigres * (int(cfreq) +( ( cfreq - int(cfreq) >= 0.5 ) ? 1 : 0 )) ;
+          double cfreq = ( abcissa[abcissa.size()-1] - maxfreq + 0.5 * resol ) / resol ;
+          nmaxchan = abcissa.size() - 1 - 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 {
@@ -3436 +3383 @@
         }
       }
-//       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 ) ;
+      }    
     }
 
@@ -3484 +3431 @@
             if ( nchan < minchan ) {
               minchan = nchan ;
-              maxdnu = abs(dnu) ;
+              maxdnu = 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 ;
             }
           }
@@ -3505 +3449 @@
           //os << "   regridChannel applied to " ;
           //if ( tableid != reftable ) 
-          refreqid = newin[tableid]->frequencies().addEntry( refpixref, refvalref, maxdnu ) ;
+          refreqid = newin[tableid]->frequencies().addEntry( refpixref, refvalref, refinc ) ;
           for ( uInt irow = 0 ; irow < newin[tableid]->table().nrow() ; irow++ ) {
             uInt tfreqid = freqIdVec[tableid][irow] ;
@@ -3526 +3470 @@
           vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
           minchan = abcissa.size() ;
-          maxdnu = abs( abcissa[1] - abcissa[0]) ;
+          maxdnu = abcissa[1] - abcissa[0] ;
         }
       }
@@ -3596 +3540 @@
           vector<double> abcissa = tmpout->getAbcissa( irow ) ;
           double bw = ( abcissa[1] - abcissa[0] ) * abcissa.size() ;
-          int nchan = (int) abs( bw / gmaxdnu[igrp] ) ;
-          // All spectra will have positive frequency increments
+          int nchan = (int)( bw / gmaxdnu[igrp] ) ;
           tmpout->regridChannel( nchan, gmaxdnu[igrp], irow ) ;
           break ;
@@ -3614 +3557 @@
     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() ) {
@@ -3631 +3575 @@
       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 ;
@@ -3670 +3614 @@
 //       }
       // get a list of number of channels for each frequency group member
-      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 ;
-            }
-          }
-        }
+      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 ;
       }
       // combine spectra
@@ -3688 +3635 @@
         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++ ) {
@@ -3699 +3646 @@
             for ( uInt jrow = 0 ; jrow < iter.table().nrow() ; jrow++ ) {
               uInt ifno = ifnoCol( jrow ) ;
-              // 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 ) ;
+              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] ) {
+              if ( ifno == freqgrp[igrp][imem] && dd <= tol ) {
                 Vector<Float> spec = specCols( jrow ) ;
                 Vector<uChar> flag = flagCols( jrow ) ;
@@ -3728 +3671 @@
           specColOut.put( irow, newspec ) ;
           flagColOut.put( irow, newflag ) ;
-          // IFNO renumbering (renumbered as frequency group ID)
+          // IFNO renumbering
           ifnoColOut.put( irow, igrp ) ; 
         }
@@ -3739 +3682 @@
       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 ) == igrp && !updated[igrp] ) {
+        if ( ifnoColOut( irow ) == gmemid[igrp] && !updated[igrp] ) {
           uInt freqidOut = freqidColOut( irow ) ;
           //os << "freqgrp " << igrp << " freqidOut = " << freqidOut << LogIO::POST ;
@@ -3752 +3693 @@
           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, gmaxdnu[igrp], freqidOut ) ;
+          out->frequencies().setEntry( refpix, refval, increm, freqidOut ) ;
           updated[igrp] = true ;
         }