Changeset 1446 for branches/alma/src

Timestamp:

11/12/08 17:04:01 (16 years ago)

Author:

TakTsutsumi

Message:

Merged recent updates (since 2007) from nrao-asap

Location:

branches/alma/src

Files:

• Makefile (added)
• MathUtils.cpp (modified) (2 diffs)
• MathUtils.h (modified) (1 diff)
• RowAccumulator.cpp (modified) (1 diff)
• RowAccumulator.h (modified) (1 diff)
• STAsciiWriter.cpp (modified) (1 diff)
• STFiller.cpp (modified) (2 diffs)
• STFrequencies.cpp (modified) (3 diffs)
• STFrequencies.h (modified) (3 diffs)
• STMath.cpp (modified) (5 diffs)
• STMath.h (modified) (6 diffs)
• STMathWrapper.h (modified) (1 diff)
• STMolecules.cpp (modified) (7 diffs)
• STMolecules.h (modified) (3 diffs)
• STWriter.cpp (modified) (8 diffs)
• Scantable.cpp (modified) (9 diffs)
• Scantable.h (modified) (6 diffs)
• ScantableWrapper.h (modified) (4 diffs)
• Templates.cpp (modified) (1 diff)
• python_STMath.cpp (modified) (1 diff)
• python_Scantable.cpp (modified) (3 diffs)

branches/alma/src/MathUtils.cpp

@@ -37 +37 @@
 #include <casa/BasicSL/String.h>
 #include <scimath/Mathematics/MedianSlider.h>
+#include <casa/Exceptions/Error.h>
 
 #include "MathUtils.h"
@@ -47 +48 @@
    String str(which);
    str.upcase();
-   if (str.contains(String("MIN"))) {
+   if (str.matches(String("MIN"))) {
       return min(data);
-   } else if (str.contains(String("MAX"))) {
+   } else if (str.matches(String("MAX"))) {
       return max(data);
-   } else if (str.contains(String("SUMSQ"))) {
+   } else if (str.matches(String("SUMSQ"))) {
       return sumsquares(data);
-   } else if (str.contains(String("SUM"))) {
+   } else if (str.matches(String("SUM"))) {
       return sum(data);
-   } else if (str.contains(String("MEAN"))) {
+   } else if (str.matches(String("MEAN"))) {
       return mean(data);
-   } else if (str.contains(String("VAR"))) {
+   } else if (str.matches(String("VAR"))) {
       return variance(data);
-   } else if (str.contains(String("STDDEV"))) {
+      } else if (str.matches(String("STDDEV"))) {
       return stddev(data);
-   } else if (str.contains(String("AVDEV"))) {
+   } else if (str.matches(String("AVDEV"))) {
       return avdev(data);
-   } else if (str.contains(String("RMS"))) {
+   } else if (str.matches(String("RMS"))) {
       uInt n = data.nelementsValid();
       return sqrt(sumsquares(data)/n);
-   } else if (str.contains(String("MED"))) {
+   } else if (str.matches(String("MEDIAN"))) {
       return median(data);
-   }
+   } else {
+      String msg = str + " is not a valid type of statistics";
+      throw(AipsError(msg));
+   } 
    return 0.0;
 }

branches/alma/src/MathUtils.h

@@ -82 +82 @@
 
 /**
- * Convert casa implementations to stl
- * @param in casa string
- * @return a std vector of std strings
+ * Convert a std::vector of std::string
+ * to a casa::Vector casa::String
+ * @param in
+ * @return
  */
 std::vector<std::string> tovectorstring(const casa::Vector<casa::String>& in);
 
 /**
- * convert stl implementations to casa versions
+ * Convert a casa::Vector of casa::String
+ * to a stl std::vector of stl std::string
  * @param in
  * @return

branches/alma/src/RowAccumulator.cpp

@@ -152 +152 @@
   userMask_ = m;
 }
+
+// Added by TT  check the state of RowAccumulator
+casa::Bool RowAccumulator::state() const
+{
+  return initialized_;
+}

branches/alma/src/RowAccumulator.h

@@ -85 +85 @@
     */
   void reset();
+  /**
+    * check the initialization state 
+    */ 
+  casa::Bool state() const;
 
 private:

branches/alma/src/STAsciiWriter.cpp

@@ -122 +122 @@
     String wcs = stable.frequencies().print(rec.asuInt("FREQ_ID"), True);
     addLine(of, "WCS", wcs);
-    addLine(of, "Rest Freq.", 
-            stable.molecules().getRestFrequency(rec.asuInt("MOLECULE_ID") ));
+    std::vector<double> restfreqs= stable.molecules().getRestFrequency(rec.asuInt("MOLECULE_ID"));
+    int nf = restfreqs.size();
+    //addLine(of, "Rest Freq.", 
+    //        stable.molecules().getRestFrequency(rec.asuInt("MOLECULE_ID") ));
+    addLine(of, "Rest Freq.", restfreqs[0]);
+    for ( unsigned int i=1; i<nf; ++i) {
+        addLine(of, " ", restfreqs[i]);
+    } 
     of << setfill('#') << setw(70) << "" << setfill(' ') << endl;
 

branches/alma/src/STFiller.cpp

@@ -205 +205 @@
     }
   }
+  String freqFrame = header_->freqref;
+  //translate frequency reference frame back to 
+  //MS style (as PKSMS2reader converts the original frame 
+  //in FITS standard style)
+  if (freqFrame == "TOPOCENT") {
+    freqFrame = "TOPO";
+  } else if (freqFrame == "GEOCENER") {
+    freqFrame = "GEO";
+  } else if (freqFrame == "BARYCENT") {
+    freqFrame = "BARY";
+  } else if (freqFrame == "GALACTOC") {
+    freqFrame = "GALACTO";
+  } else if (freqFrame == "LOCALGRP") {
+    freqFrame = "LGROUP";
+  } else if (freqFrame == "CMBDIPOL") {
+    freqFrame = "CMB";
+  } else if (freqFrame == "SOURCE") {
+    freqFrame = "REST";
+  }
+  table_->frequencies().setFrame(freqFrame);
      
 }
@@ -222 +242 @@
   Float  azimuth, elevation, focusAxi, focusRot, focusTan,
     humidity, parAngle, pressure, temperature, windAz, windSpeed;
-  Double bandwidth, freqInc, interval, mjd, refFreq, restFreq, srcVel;
+  Double bandwidth, freqInc, interval, mjd, refFreq, srcVel;
   String          fieldName, srcName, tcalTime, obsType;
   Vector<Float>   calFctr, sigma, tcal, tsys;
   Matrix<Float>   baseLin, baseSub;
-  Vector<Double>  direction(2), scanRate(2), srcDir(2), srcPM(2);
+  Vector<Double>  direction(2), scanRate(2), srcDir(2), srcPM(2), restFreq;
   Matrix<Float>   spectra;
   Matrix<uChar>   flagtra;

branches/alma/src/STFrequencies.cpp

@@ -128 +128 @@
 }
 
+void STFrequencies::setEntry( Double refpix, Double refval, Double inc, uInt id )
+{
+  Table t = table_(table_.col("ID") == Int(id) );
+  if (t.nrow() == 0 ) {
+    throw(AipsError("STFrequencies::getEntry - freqID out of range"));
+  }
+  for ( uInt i = 0 ; i < table_.nrow() ; i++ ) {
+    uInt fid ;
+    idCol_.get( i, fid ) ;
+    if ( fid == id ) {
+      refpixCol_.put( i, refpix ) ;
+      refvalCol_.put( i, refval ) ;
+      incrCol_.put( i, inc ) ;
+    }
+  }
+}
+
 SpectralCoordinate STFrequencies::getSpectralCoordinate( uInt id ) const
 {
@@ -140 +157 @@
   // get first row - there should only be one matching id
   const TableRecord& rec = row.get(0);
+
   return SpectralCoordinate( getFrame(true), rec.asDouble("REFVAL"),
                              rec.asDouble("INCREMENT"),
@@ -145 +163 @@
 }
 
+/**
 SpectralCoordinate
-  STFrequencies::getSpectralCoordinate( const MDirection& md,
-                                        const MPosition& mp,
-                                        const MEpoch& me,
-                                        Double restfreq, uInt id ) const
+  asap::STFrequencies::getSpectralCoordinate( const MDirection& md,
+                                              const MPosition& mp,
+                                              const MEpoch& me,
+                                              Double restfreq, uInt id ) const
+**/
+SpectralCoordinate
+  asap::STFrequencies::getSpectralCoordinate( const MDirection& md,
+                                              const MPosition& mp,
+                                              const MEpoch& me,
+                                              Vector<Double> restfreq, uInt id ) const
 {
   SpectralCoordinate spc = getSpectralCoordinate(id);
-  spc.setRestFrequency(restfreq, True);
+  //spc.setRestFrequency(restfreq, True);
+  // for now just use the first rest frequency 
+  if (restfreq.nelements()==0 ) {
+    restfreq.resize(1);
+    restfreq[0] = 0;
+  }
+  spc.setRestFrequency(restfreq[0], True);
   if ( !spc.setReferenceConversion(getFrame(), me, mp, md) ) {
     throw(AipsError("Couldn't convert frequency frame."));

branches/alma/src/STFrequencies.h

@@ -59 +59 @@
                  casa::Double& inc, casa::uInt id );
 
+  /***
+   * Set the frequency values for a specific id via references
+   * @param refpix the reference pixel
+   * @param refval the reference value
+   * @param inc    the increment
+   * @param id     the identifier
+   *
+   * 17/09/2008 Takeshi Nakazato
+   ***/
+  void setEntry( casa::Double refpix, casa::Double refval, 
+                 casa::Double inc, casa::uInt id ) ;
+
 
   bool conformant(const STFrequencies& other) const;
@@ -69 +81 @@
   casa::SpectralCoordinate getSpectralCoordinate( casa::uInt freqID ) const;
 
+  /**
   casa::SpectralCoordinate getSpectralCoordinate( const casa::MDirection& md,
                                                   const casa::MPosition& mp,
@@ -75 +88 @@
                                                   casa::uInt freqID
                                                   ) const;
+  **/
+  casa::SpectralCoordinate getSpectralCoordinate( const casa::MDirection& md,
+                                                  const casa::MPosition& mp,
+                                                  const casa::MEpoch& me,
+                                                  casa::Vector<casa::Double> restfreq,
+                                                  casa::uInt freqID
+                                                  ) const;
 
   /**

branches/alma/src/STMath.cpp

@@ -191 +191 @@
     }
     //write out
-    Vector<uChar> flg(msk.shape());
-    convertArray(flg, !msk);
-    flagColOut.put(i, flg);
-    specColOut.put(i, acc.getSpectrum());
-    tsysColOut.put(i, acc.getTsys());
-    intColOut.put(i, acc.getInterval());
-    mjdColOut.put(i, acc.getTime());
-    // we should only have one cycle now -> reset it to be 0
-    // frequency switched data has different CYCLENO for different IFNO
-    // which requires resetting this value
-    cycColOut.put(i, uInt(0));
+    if (acc.state()) {
+      Vector<uChar> flg(msk.shape());
+      convertArray(flg, !msk);
+      flagColOut.put(i, flg);
+      specColOut.put(i, acc.getSpectrum());
+      tsysColOut.put(i, acc.getTsys());
+      intColOut.put(i, acc.getInterval());
+      mjdColOut.put(i, acc.getTime());
+      // we should only have one cycle now -> reset it to be 0
+      // frequency switched data has different CYCLENO for different IFNO
+      // which requires resetting this value
+      cycColOut.put(i, uInt(0));
+    } else {
+      ostringstream oss;
+      oss << "For output row="<<i<<", all input rows of data are flagged. no averaging" << endl;
+      pushLog(String(oss));
+    }
     acc.reset();
   }
@@ -947 +953 @@
 {
 
-
+  
   STSelector sel;
   CountedPtr< Scantable > ws = getScantable(s, false);
   CountedPtr< Scantable > sig, sigwcal, ref, refwcal;
-  CountedPtr< Scantable > calsig, calref, out;
+  CountedPtr< Scantable > calsig, calref, out, out1, out2;
+  Bool nofold=False;
 
   //split the data
@@ -973 +980 @@
   calref = dototalpower(refwcal, ref, tcal=tcal);
 
-  out=dosigref(calsig,calref,smoothref,tsysv,tau);
-
+  out1=dosigref(calsig,calref,smoothref,tsysv,tau); 
+  out2=dosigref(calref,calsig,smoothref,tsysv,tau); 
+
+  Table& tabout1=out1->table();
+  Table& tabout2=out2->table();
+  ROScalarColumn<uInt> freqidCol1(tabout1, "FREQ_ID");
+  ROScalarColumn<uInt> freqidCol2(tabout2, "FREQ_ID");
+  ROArrayColumn<Float> specCol(tabout2, "SPECTRA");
+  Vector<Float> spec; specCol.get(0, spec);
+  uInt nchan = spec.nelements();
+  uInt freqid1; freqidCol1.get(0,freqid1);
+  uInt freqid2; freqidCol2.get(0,freqid2);
+  Double rp1, rp2, rv1, rv2, inc1, inc2;
+  out1->frequencies().getEntry(rp1, rv1, inc1, freqid1);
+  out2->frequencies().getEntry(rp2, rv2, inc2, freqid2);
+  if (rp1==rp2) {
+    Double foffset = rv1 - rv2;
+    uInt choffset = static_cast<uInt>(foffset/abs(inc2));
+    if (choffset >= nchan) {
+      cerr<<"out-band frequency switching, no folding"<<endl;
+      nofold = True;
+    }
+  }
+ 
+  if (nofold) {
+    std::vector< CountedPtr< Scantable > > tabs; 
+    tabs.push_back(out1);
+    tabs.push_back(out2);
+    out = merge(tabs);
+  }
+  else { //folding is not implemented yet
+    out = out1;
+  }
+    
   return out;
 }
-
 
 CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
@@ -1565 +1603 @@
           id = out->frequencies().addEntry(rp, rv, inc);
           freqidcol.put(k,id);
-          String name,fname;Double rf;
+          //String name,fname;Double rf;
+          Vector<String> name,fname;Vector<Double> rf;
           (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
           id = out->molecules().addEntry(rf, name, fname);
@@ -1981 +2020 @@
   return out;
 }
+
+// Averaging spectra with different channel/resolution
+CountedPtr<Scantable>
+STMath::new_average( const std::vector<CountedPtr<Scantable> >& in,
+                     const bool& compel,
+                     const std::vector<bool>& mask,
+                     const std::string& weight,
+                     const std::string& avmode )
+  throw ( casa::AipsError )
+{
+  if ( avmode == "SCAN" && in.size() != 1 )
+    throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
+                    "Use merge first."));
+  
+  CountedPtr<Scantable> out ;     // processed result 
+  if ( compel ) {
+    std::vector< CountedPtr<Scantable> > newin ; // input for average process
+    uInt insize = in.size() ;    // number of scantables 
+
+    // TEST: do normal average in each table before IF grouping
+    cout << "Do preliminary averaging" << endl ;
+    vector< CountedPtr<Scantable> > tmpin( insize ) ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      vector< CountedPtr<Scantable> > v( 1, in[itable] ) ;
+      tmpin[itable] = average( v, mask, weight, avmode ) ;
+    }
+
+    // warning
+    cout << "Average spectra with different spectral resolution" << endl ;
+    cout << endl ;
+
+    // temporarily set coordinfo
+    vector<string> oldinfo( insize ) ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      vector<string> coordinfo = in[itable]->getCoordInfo() ;
+      oldinfo[itable] = coordinfo[0] ;
+      coordinfo[0] = "Hz" ;
+      tmpin[itable]->setCoordInfo( coordinfo ) ;
+    }
+
+    // columns
+    ScalarColumn<uInt> freqIDCol ;
+    ScalarColumn<uInt> ifnoCol ;
+    ScalarColumn<uInt> scannoCol ;
+
+
+    // check IF frequency coverage
+    //cout << "Check IF settings in each table" << endl ;
+    vector< vector<uInt> > freqid( insize );
+    vector< vector<double> > iffreq( insize ) ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      uInt rows = tmpin[itable]->nrow() ;
+      uInt freqnrows = tmpin[itable]->frequencies().table().nrow() ;
+      for ( uInt irow = 0 ; irow < rows ; irow++ ) {
+        if ( freqid[itable].size() == freqnrows ) {
+          break ;
+        }
+        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 ) {
+            //cout << "itable = " << itable << ": IF " << id << " is included in the list" << endl ;
+            vector<double> abcissa = tmpin[itable]->getAbcissa( irow ) ;
+            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] ) ) ;
+          }
+        }
+      }
+    }
+
+    // debug
+    //cout << "IF settings summary:" << endl ;
+    //for ( uInt i = 0 ; i < freqid.size() ; i++ ) {
+    //cout << "   Table" << i << endl ;
+    //for ( uInt j = 0 ; j < freqid[i].size() ; j++ ) {
+    //cout << "      id = " << freqid[i][j] << " (min,max) = (" << iffreq[i][2*j] << "," << iffreq[i][2*j+1] << ")" << endl ;
+    //}
+    //}
+    //cout << endl ;
+
+    // IF grouping based on their frequency coverage
+    //cout << "IF grouping based on their frequency coverage" << endl ;
+
+    // IF group 
+    // ifgrp[numgrp][nummember*2]
+    // ifgrp = [table00, freqrow00, table01, freqrow01, ...],
+    //         [table11, freqrow11, table11, freqrow11, ...],
+    //         ...
+    // ifgfreq[numgrp*2] 
+    // ifgfreq = [min0, max0, min1, max1, ...]
+    vector< vector<uInt> > ifgrp ;
+    vector<double> ifgfreq ;
+
+    // parameter for IF grouping
+    // groupmode = OR    retrieve all region 
+    //             AND   only retrieve overlaped region
+    //string groupmode = "AND" ;
+    string groupmode = "OR" ;
+    uInt sizecr = 0 ;
+    if ( groupmode == "AND" ) 
+      sizecr = 2 ;
+    else if ( groupmode == "OR" ) 
+      sizecr = 0 ;
+
+    vector<double> sortedfreq ;
+    for ( uInt i = 0 ; i < iffreq.size() ; i++ ) {
+      for ( uInt j = 0 ; j < iffreq[i].size() ; j++ ) {
+        if ( count( sortedfreq.begin(), sortedfreq.end(), iffreq[i][j] ) == 0 )
+          sortedfreq.push_back( iffreq[i][j] ) ;
+      }
+    }
+    sort( sortedfreq.begin(), sortedfreq.end() ) ;
+    for ( vector<double>::iterator i = sortedfreq.begin() ; i != sortedfreq.end()-1 ; i++ ) {
+      ifgfreq.push_back( *i ) ;
+      ifgfreq.push_back( *(i+1) ) ;
+    }
+    ifgrp.resize( ifgfreq.size()/2 ) ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      for ( uInt iif = 0 ; iif < freqid[itable].size() ; iif++ ) {
+        double range0 = iffreq[itable][2*iif] ;
+        double range1 = iffreq[itable][2*iif+1] ;
+        for ( uInt j = 0 ; j < ifgrp.size() ; j++ ) {
+          double fmin = max( range0, ifgfreq[2*j] ) ;
+          double fmax = min( range1, ifgfreq[2*j+1] ) ;
+          if ( fmin < fmax ) {
+            ifgrp[j].push_back( itable ) ;
+            ifgrp[j].push_back( freqid[itable][iif] ) ;
+          }
+        }
+      }
+    }
+    vector< vector<uInt> >::iterator fiter = ifgrp.begin() ; 
+    vector<double>::iterator giter = ifgfreq.begin() ;
+    while( fiter != ifgrp.end() ) {
+      if ( fiter->size() <= sizecr ) {
+        fiter = ifgrp.erase( fiter ) ;
+        giter = ifgfreq.erase( giter ) ;
+        giter = ifgfreq.erase( giter ) ;
+      }
+      else {
+        fiter++ ;
+        advance( giter, 2 ) ;
+      }
+    }
+
+    // freqgrp[numgrp][nummember]
+    // freqgrp = [ifgrp00, ifgrp01, ifgrp02, ...],
+    //           [ifgrp10, ifgrp11, ifgrp12, ...],
+    //           ...
+    // freqrange[numgrp*2]
+    // freqrange = [min0, max0, min1, max1, ...]
+    vector< vector<uInt> > freqgrp ;
+    double freqrange = 0.0 ;
+    uInt grpnum = 0 ;
+    for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
+      // Assumed that ifgfreq was sorted
+      if ( grpnum != 0 && freqrange == ifgfreq[2*i] ) {
+        freqgrp[grpnum-1].push_back( i ) ;
+      }
+      else {
+        vector<uInt> grp0( 1, i ) ;
+        freqgrp.push_back( grp0 ) ;
+        grpnum++ ;
+      }
+      freqrange = ifgfreq[2*i+1] ;
+    }
+        
+
+    // print IF groups
+    cout << "IF Group summary: " << endl ;
+    cout << "   GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
+    for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
+      cout << "   GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
+      for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
+        cout << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ; 
+      }
+      cout << endl ;
+    }
+    cout << endl ;
+    
+    // print frequency group
+    cout << "Frequency Group summary: " << endl ;
+    cout << "   GROUP_ID [FREQ_MIN, FREQ_MAX]: IF_GROUP_ID" << endl ;
+    for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
+      cout << "   GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*freqgrp[i][0]] << "," << ifgfreq[2*freqgrp[i][freqgrp[i].size()-1]+1] << "]: " ;
+      for ( uInt j = 0 ; j < freqgrp[i].size() ; j++ ) {
+        cout << freqgrp[i][j] << " " ;
+      }
+      cout << endl ;
+    }
+    cout << endl ;
+
+    // membership check
+    // groups[numtable][numIF][nummembership]
+    vector< vector< vector<uInt> > > groups( insize ) ;
+    for ( uInt i = 0 ; i < insize ; i++ ) {
+      groups[i].resize( freqid[i].size() ) ;
+    }
+    for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
+      for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
+        uInt tableid = ifgrp[igrp][2*imem] ;
+        vector<uInt>::iterator iter = find( freqid[tableid].begin(), freqid[tableid].end(), ifgrp[igrp][2*imem+1] ) ;
+        if ( iter != freqid[tableid].end() ) {
+          uInt rowid = distance( freqid[tableid].begin(), iter ) ;
+          groups[tableid][rowid].push_back( igrp ) ;
+        }
+      }
+    }
+
+    // print membership
+    //for ( uInt i = 0 ; i < insize ; i++ ) {
+    //cout << "Table " << i << endl ;
+    //for ( uInt j = 0 ; j < groups[i].size() ; j++ ) {
+    //cout << "   FREQ_ID " <<  setw( 2 ) << freqid[i][j] << ": " ;
+    //for ( uInt k = 0 ; k < groups[i][j].size() ; k++ ) {
+    //cout << setw( 2 ) << groups[i][j][k] << " " ;
+    //}
+    //cout << endl ;
+    //}
+    //}
+
+    // set back coordinfo
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      vector<string> coordinfo = tmpin[itable]->getCoordInfo() ;
+      coordinfo[0] = oldinfo[itable] ;
+      tmpin[itable]->setCoordInfo( coordinfo ) ;
+    }
+
+    // Create additional table if needed
+    bool oldInsitu = insitu_ ;
+    setInsitu( false ) ;
+    vector< vector<uInt> > addrow( insize ) ;
+    vector<uInt> addtable( insize, 0 ) ;
+    vector<uInt> newtableids( insize ) ;
+    vector<uInt> newifids( insize, 0 ) ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      //cout << "Table " << setw(2) << itable << ": " ;
+      for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
+        addrow[itable].push_back( groups[itable][ifrow].size()-1 ) ;
+        //cout << addrow[itable][ifrow] << " " ;
+      }
+      addtable[itable] = *max_element( addrow[itable].begin(), addrow[itable].end() ) ;
+      //cout << "(" << addtable[itable] << ")" << endl ;
+    }
+    newin.resize( insize ) ;
+    copy( tmpin.begin(), tmpin.end(), newin.begin() ) ;
+    for ( uInt i = 0 ; i < insize ; i++ ) {
+      newtableids[i] = i ;
+    }
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
+        CountedPtr<Scantable> add = getScantable( newin[itable], false ) ;
+        vector<int> freqidlist ;
+        for ( uInt i = 0 ; i < groups[itable].size() ; i++ ) {
+          if ( groups[itable][i].size() > iadd + 1 ) {
+            freqidlist.push_back( freqid[itable][i] ) ;
+          }
+        }
+        stringstream taqlstream ;
+        taqlstream << "SELECT FROM $1 WHERE FREQ_ID IN [" ;
+        for ( uInt i = 0 ; i < freqidlist.size() ; i++ ) {
+          taqlstream << i ;
+          if ( i < freqidlist.size() - 1 )
+            taqlstream << "," ;
+          else
+            taqlstream << "]" ;
+        }
+        string taql = taqlstream.str() ;
+        //cout << "taql = " << taql << endl ;
+        STSelector selector = STSelector() ;
+        selector.setTaQL( taql ) ;
+        add->setSelection( selector ) ;
+        newin.push_back( add ) ;
+        newtableids.push_back( itable ) ;
+        newifids.push_back( iadd + 1 ) ;
+      }
+    }
+
+    // udpate ifgrp
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
+        for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
+          if ( groups[itable][ifrow].size() > iadd + 1 ) {
+            uInt igrp = groups[itable][ifrow][iadd+1] ;
+            for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
+              if ( ifgrp[igrp][2*imem] == newtableids[iadd+insize] && ifgrp[igrp][2*imem+1] == freqid[newtableids[iadd+insize]][ifrow] ) {
+                ifgrp[igrp][2*imem] = insize + iadd ;
+              }
+            }
+          }
+        }
+      }
+    }
+
+    // print IF groups again for debug
+    //cout << "IF Group summary: " << endl ;
+    //cout << "   GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
+    //for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
+    //cout << "   GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
+    //for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
+    //cout << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ; 
+    //}
+    //cout << endl ;
+    //}
+    //cout << endl ;
+
+    // reset SCANNO and IF: IF number is reset by the result of sortation 
+    cout << "All scan number is set to 0" << endl ;
+    //cout << "All IF number is set to IF group index" << endl ;
+    cout << endl ;
+    insize = newin.size() ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      uInt rows = newin[itable]->nrow() ;
+      Table &tmpt = newin[itable]->table() ;
+      freqIDCol.attach( tmpt, "FREQ_ID" ) ;
+      scannoCol.attach( tmpt, "SCANNO" ) ;
+      ifnoCol.attach( tmpt, "IFNO" ) ;
+      for ( uInt irow=0 ; irow < rows ; irow++ ) {
+        scannoCol.put( irow, 0 ) ;
+        uInt freqID = freqIDCol( irow ) ;
+        vector<uInt>::iterator iter = find( freqid[newtableids[itable]].begin(), freqid[newtableids[itable]].end(), freqID ) ;
+        if ( iter != freqid[newtableids[itable]].end() ) {
+          uInt index = distance( freqid[newtableids[itable]].begin(), iter ) ;
+          ifnoCol.put( irow, groups[newtableids[itable]][index][newifids[itable]] ) ;
+        }
+        else {
+          throw(AipsError("IF grouping was wrong in additional tables.")) ;
+        } 
+      }
+    }
+    oldinfo.resize( insize ) ;
+    setInsitu( oldInsitu ) ;
+
+    // temporarily set coordinfo
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      vector<string> coordinfo = newin[itable]->getCoordInfo() ;
+      oldinfo[itable] = coordinfo[0] ;
+      coordinfo[0] = "Hz" ;
+      newin[itable]->setCoordInfo( coordinfo ) ;
+    }
+
+    // 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 ) ) ;
+        ifNoVec[itable].push_back( ifnoCol( irow ) ) ;
+      }
+    }
+
+    // reset spectra and flagtra: pick up common part of frequency coverage
+    //cout << "Pick common frequency range and align resolution" << endl ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      uInt rows = newin[itable]->nrow() ;
+      int nminchan = -1 ;
+      int nmaxchan = -1 ;
+      vector<uInt> freqIdUpdate ; 
+      for ( uInt irow = 0 ; irow < rows ; irow++ ) {
+        uInt ifno = ifNoVec[itable][irow] ;  // IFNO is reset by group index
+        double minfreq = ifgfreq[2*ifno] ;
+        double maxfreq = ifgfreq[2*ifno+1] ;
+        //cout << "frequency range: [" << minfreq << "," << maxfreq << "]" << endl ;
+        vector<double> abcissa = newin[itable]->getAbcissa( irow ) ;
+        int nchan = abcissa.size() ;
+        double resol = abcissa[1] - abcissa[0] ;
+        //cout << "abcissa range  : [" << abcissa[0] << "," << abcissa[nchan-1] << "]" << endl ;
+        if ( minfreq <= abcissa[0] ) 
+          nminchan = 0 ;
+        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 ) ;
+        }
+        if ( maxfreq >= abcissa[abcissa.size()-1] )
+          nmaxchan = abcissa.size() - 1 ;
+        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 ) ;
+        }
+        //cout << "channel range (" << irow << "): [" << nminchan << "," << nmaxchan << "]" << endl ;
+        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 )
+            freqIdUpdate.push_back( freqIdVec[itable][irow] ) ;
+        }
+        else {
+          throw(AipsError("Failed to pick up common part of frequency range.")) ;
+        }
+      }
+      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 ) ;
+      }    
+    }
+
+    
+    // reset spectra and flagtra: align spectral resolution
+    //cout << "Align spectral resolution" << endl ;
+    vector<double> gmaxdnu( freqgrp.size(), 0.0 ) ;
+    vector<uInt> gmemid( freqgrp.size(), 0 ) ;
+    for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
+      double maxdnu = 0.0 ;       // maximum (coarsest) frequency resolution 
+      int minchan = INT_MAX ;     // minimum channel number
+      Double refpixref = -1 ;     // reference of 'reference pixel'
+      Double refvalref = -1 ;     // reference of 'reference frequency'
+      Double refinc = -1 ;        // reference frequency resolution 
+      uInt refreqid ;
+      uInt reftable = INT_MAX;
+      // process only if group member > 1
+      if ( ifgrp[igrp].size() > 2 ) {
+        // find minchan and maxdnu in each group
+        for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
+          uInt tableid = ifgrp[igrp][2*imem] ;
+          uInt rowid = ifgrp[igrp][2*imem+1] ;
+          vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
+          if ( iter != freqIdVec[tableid].end() ) {
+            uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
+            vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
+            int nchan = abcissa.size() ;
+            double dnu = abcissa[1] - abcissa[0] ;
+            //cout << "GROUP " << igrp << " (" << tableid << "," << rowid << "): nchan = " << nchan << " (minchan = " << minchan << ")" << endl ; 
+            if ( nchan < minchan ) {
+              minchan = nchan ;
+              maxdnu = dnu ;
+              newin[tableid]->frequencies().getEntry( refpixref, refvalref, refinc, rowid ) ;
+              refreqid = rowid ;
+              reftable = tableid ;
+            }
+          }
+        }
+        // regrid spectra in each group
+        cout << "GROUP " << igrp << endl ;
+        cout << "   Channel number is adjusted to " << minchan << endl ;
+        cout << "   Corresponding frequency resolution is " << maxdnu << "Hz" << endl ;
+        for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
+          uInt tableid = ifgrp[igrp][2*imem] ;
+          uInt rowid = ifgrp[igrp][2*imem+1] ;
+          freqIDCol.attach( newin[tableid]->table(), "FREQ_ID" ) ;
+          //cout << "tableid = " << tableid << " rowid = " << rowid << ": " << endl ;
+          //cout << "   regridChannel applied to " ;
+          if ( tableid != reftable ) 
+            refreqid = newin[tableid]->frequencies().addEntry( refpixref, refvalref, refinc ) ;
+          for ( uInt irow = 0 ; irow < newin[tableid]->table().nrow() ; irow++ ) {
+            uInt tfreqid = freqIdVec[tableid][irow] ;
+            if ( tfreqid == rowid ) {     
+              //cout << irow << " " ;
+              newin[tableid]->regridChannel( minchan, maxdnu, irow ) ;
+              freqIDCol.put( irow, refreqid ) ;
+              freqIdVec[tableid][irow] = refreqid ;
+            }
+          }
+          //cout << endl ;
+        } 
+      }
+      else {
+        uInt tableid = ifgrp[igrp][0] ;
+        uInt rowid = ifgrp[igrp][1] ;
+        vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
+        if ( iter != freqIdVec[tableid].end() ) {
+          uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
+          vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
+          minchan = abcissa.size() ;
+          maxdnu = abcissa[1] - abcissa[0] ;
+        }
+      }
+      for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
+        if ( count( freqgrp[i].begin(), freqgrp[i].end(), igrp ) > 0 ) {
+          if ( maxdnu > gmaxdnu[i] ) {
+            gmaxdnu[i] = maxdnu ;
+            gmemid[i] = igrp ;
+          }
+          break ;
+        }
+      }
+    }
+    cout << endl ;
+
+    // set back coordinfo
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      vector<string> coordinfo = newin[itable]->getCoordInfo() ;
+      coordinfo[0] = oldinfo[itable] ;
+      newin[itable]->setCoordInfo( coordinfo ) ;
+    }      
+
+    // accumulate all rows into the first table 
+    // NOTE: assumed in.size() = 1
+    vector< CountedPtr<Scantable> > tmp( 1 ) ;
+    if ( newin.size() == 1 )
+      tmp[0] = newin[0] ;
+    else 
+      tmp[0] = merge( newin ) ;
+
+    //return tmp[0] ;
+
+    // average
+    CountedPtr<Scantable> tmpout = average( tmp, mask, weight, avmode ) ;
+
+    //return tmpout ;
+
+    // combine frequency group
+    cout << "Combine spectra based on frequency grouping" << endl ;
+    cout << "IFNO is renumbered as frequency group ID (see above)" << endl ;
+    vector<string> coordinfo = tmpout->getCoordInfo() ;
+    oldinfo[0] = coordinfo[0] ;
+    coordinfo[0] = "Hz" ;
+    tmpout->setCoordInfo( coordinfo ) ;
+    // create proformas of output table
+    stringstream taqlstream ;
+    taqlstream << "SELECT FROM $1 WHERE IFNO IN [" ;
+    for ( uInt i = 0 ; i < gmemid.size() ; i++ ) {
+      taqlstream << gmemid[i] ;
+      if ( i < gmemid.size() - 1 )
+        taqlstream << "," ;
+      else
+        taqlstream << "]" ;
+    }
+    string taql = taqlstream.str() ;
+    //cout << "taql = " << taql << endl ;
+    STSelector selector = STSelector() ;
+    selector.setTaQL( taql ) ;
+    oldInsitu = insitu_ ;
+    setInsitu( false ) ;
+    out = getScantable( tmpout, false ) ;
+    setInsitu( oldInsitu ) ;
+    out->setSelection( selector ) ;
+    // regrid rows
+    ifnoCol.attach( tmpout->table(), "IFNO" ) ;
+    for ( uInt irow = 0 ; irow < tmpout->table().nrow() ; irow++ ) {
+      uInt ifno = ifnoCol( irow ) ;
+      for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
+        if ( count( freqgrp[igrp].begin(), freqgrp[igrp].end(), ifno ) > 0 ) {
+          vector<double> abcissa = tmpout->getAbcissa( irow ) ;
+          double bw = ( abcissa[1] - abcissa[0] ) * abcissa.size() ;
+          int nchan = (int)( bw / gmaxdnu[igrp] ) ;
+          tmpout->regridChannel( nchan, gmaxdnu[igrp], irow ) ;
+          break ;
+        }
+      }
+    }
+    // combine spectra
+    ArrayColumn<Float> specColOut ;
+    specColOut.attach( out->table(), "SPECTRA" ) ;
+    ArrayColumn<uChar> flagColOut ;
+    flagColOut.attach( out->table(), "FLAGTRA" ) ;
+    ScalarColumn<uInt> ifnoColOut ;
+    ifnoColOut.attach( out->table(), "IFNO" ) ;
+    ScalarColumn<uInt> polnoColOut ;
+    polnoColOut.attach( out->table(), "POLNO" ) ;
+    ScalarColumn<uInt> freqidColOut ;
+    freqidColOut.attach( out->table(), "FREQ_ID" ) ;
+    Table &tab = tmpout->table() ;
+    TableIterator iter( tab, "POLNO" ) ;
+    vector< vector<uInt> > sizes( freqgrp.size() ) ;
+    while( !iter.pastEnd() ) {
+      vector< vector<Float> > specout( freqgrp.size() ) ;
+      vector< vector<uChar> > flagout( freqgrp.size() ) ;
+      ArrayColumn<Float> specCols ;
+      specCols.attach( iter.table(), "SPECTRA" ) ;
+      ArrayColumn<uChar> flagCols ;
+      flagCols.attach( iter.table(), "FLAGTRA" ) ;
+      ifnoCol.attach( iter.table(), "IFNO" ) ;
+      ScalarColumn<uInt> polnos ;
+      polnos.attach( iter.table(), "POLNO" ) ;
+      uInt polno = polnos( 0 ) ;
+      //cout << "POLNO iteration: " << polno << endl ;
+      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 ) ;
+              Vector<uChar> flag = flagCols( irow ) ;
+              vector<Float> svec ;
+              spec.tovector( svec ) ;
+              vector<uChar> fvec ;
+              flag.tovector( fvec ) ;
+              //cout << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << endl ;
+              specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ;
+              flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ;
+              //cout << "specout[" << igrp << "].size() = " << specout[igrp].size() << endl ;
+              sizes[igrp][imem] = spec.nelements() ;
+            }
+          }
+        }
+        for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
+          uInt ifout = ifnoColOut( irow ) ;
+          uInt polout = polnoColOut( irow ) ;
+          if ( ifout == gmemid[igrp] && polout == polno ) {
+            // set SPECTRA and FRAGTRA
+            Vector<Float> newspec( specout[igrp] ) ;
+            Vector<uChar> newflag( flagout[igrp] ) ;
+            specColOut.put( irow, newspec ) ;
+            flagColOut.put( irow, newflag ) ;
+            // IFNO renumbering
+            ifnoColOut.put( irow, igrp ) ;
+          }
+        } 
+      }
+      iter++ ;
+    }
+    // update FREQUENCIES subtable
+    vector<bool> updated( freqgrp.size(), false ) ;
+    for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
+      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] ) {
+          uInt freqidOut = freqidColOut( irow ) ;
+          //cout << "freqgrp " << igrp << " freqidOut = " << freqidOut << endl ;
+          double refpix ;
+          double refval ;
+          double increm ;
+          out->frequencies().getEntry( refpix, refval, increm, freqidOut ) ;
+          refpix += pixShift ;
+          out->frequencies().setEntry( refpix, refval, increm, freqidOut ) ;
+          updated[igrp] = true ;
+        }
+      }
+    }
+
+    //out = tmpout ;
+
+    coordinfo = tmpout->getCoordInfo() ;
+    coordinfo[0] = oldinfo[0] ;
+    tmpout->setCoordInfo( coordinfo ) ;
+  }
+  else {
+    // simple average
+    out =  average( in, mask, weight, avmode ) ;
+  }
+  
+  return out ;
+}

branches/alma/src/STMath.h

@@ -65 +65 @@
     * average a vector of Scantables
     * @param in the vector of Scantables to average
-    * @param mask an optional mask to apply on specific weights
+    * @param an optional mask to apply on specific weights
     * @param weight weighting scheme
     * @param avmode the mode ov averaging. Per "SCAN" or "ALL".
@@ -134 +134 @@
   /**
     * Calibrate total power scans (translated from GBTIDL)
-    * @param calon uncalibrated Scantable with CAL noise signal
+    * @param calon uncalibrated Scantable with CAL noise signal 
     * @param caloff uncalibrated Scantable with no CAL signal
     * @param tcal optional scalar Tcal, CAL temperature (K)
-    * @return casa::CountedPtr<Scantable> which holds a calibrated Scantable
+    * @return casa::CountedPtr<Scantable> which holds a calibrated Scantable 
     * (spectrum - average of the two CAL on and off spectra;
     * tsys - mean Tsys = <caloff>*Tcal/<calon-caloff> + Tcal/2)
-    */
+    */              
   casa::CountedPtr<Scantable> dototalpower( const casa::CountedPtr<Scantable>& calon,
                                             const casa::CountedPtr<Scantable>& caloff,
@@ -151 +151 @@
     * @param smoothref optional Boxcar smooth width of the reference scans
     * default: no smoothing (=1)
-    * @param tsysv optional scalar Tsys value at the zenith, required to
-    * set tau, as well
+    * @param tsysv optional scalar Tsys value at the zenith, required to 
+    * set tau, as well 
     * @param tau optional scalar Tau value
     * @return casa::CountedPtr<Scantable> which holds combined scans
@@ -163 +163 @@
                                         casa::Float tau=0.0 );
 
- /**
+  /**
     * Calibrate GBT Nod scan pairs (translated from GBTIDL)
     * @param s Scantable which contains Nod scans
@@ -170 +170 @@
     * @param tsysv optional scalar Tsys value at the zenith, required to
     * set tau, as well
-    * @param tau optional scalar Tau value
+    * @param tau optional scalar Tau value 
     * @param tcal optional scalar Tcal, CAL temperature (K)
     * @return casa::CountedPtr<Scantable> which holds calibrated scans
@@ -266 +266 @@
               double width);
 
+  // test for average spectra with different channel/resolution
+  casa::CountedPtr<Scantable>
+    new_average( const std::vector<casa::CountedPtr<Scantable> >& in,
+                 const bool& compel, 
+                 const std::vector<bool>& mask = std::vector<bool>(),
+                 const std::string& weight = "NONE",
+                 const std::string& avmode = "SCAN" )
+    throw (casa::AipsError) ;
+
+
 private:
   casa::CountedPtr<Scantable>  applyToPol( const casa::CountedPtr<Scantable>& in,

branches/alma/src/STMathWrapper.h

@@ -188 +188 @@
   { return ScantableWrapper(STMath::lagFlag(in.getCP(), frequency, width)); }
 
+  // test for average spectra with different channel/resolution
+  ScantableWrapper
+    new_average( const std::vector<ScantableWrapper>& in,
+                 const bool& compel,
+                 const std::vector<bool>& mask,
+                 const std::string& weight,
+                 const std::string& avmode )
+  {
+    std::vector<casa::CountedPtr<Scantable> > sts;
+    for (unsigned int i=0; i<in.size(); ++i) sts.push_back(in[i].getCP());
+    return ScantableWrapper(STMath::new_average(sts, compel, mask, weight, avmode));
+  }
+
 };
 

branches/alma/src/STMolecules.cpp

@@ -14 +14 @@
 #include <tables/Tables/SetupNewTab.h>
 #include <tables/Tables/ScaColDesc.h>
+#include <tables/Tables/ArrColDesc.h>
 #include <tables/Tables/TableRecord.h>
 #include <tables/Tables/TableParse.h>
 #include <tables/Tables/TableRow.h>
 #include <casa/Containers/RecordField.h>
+
+#include <tables/Tables/TableProxy.h>
 
 #include "STMolecules.h"
@@ -59 +62 @@
 {
   // add to base class table
-  table_.addColumn(ScalarColumnDesc<Double>("RESTFREQUENCY"));
-  table_.addColumn(ScalarColumnDesc<String>("NAME"));
-  table_.addColumn(ScalarColumnDesc<String>("FORMATTEDNAME"));
+  //table_.addColumn(ScalarColumnDesc<Double>("RESTFREQUENCY"));
+  table_.addColumn(ArrayColumnDesc<Double>("RESTFREQUENCY"));
+  //table_.addColumn(ScalarColumnDesc<String>("NAME"));
+  table_.addColumn(ArrayColumnDesc<String>("NAME"));
+  //table_.addColumn(ScalarColumnDesc<String>("FORMATTEDNAME"));
+  table_.addColumn(ArrayColumnDesc<String>("FORMATTEDNAME"));
   table_.rwKeywordSet().define("UNIT", String("Hz"));
   // new cached columns
@@ -69 +75 @@
 }
 
+/***
 uInt STMolecules::addEntry( Double restfreq, const String& name,
                             const String& formattedname )
@@ -94 +101 @@
   return resultid;
 }
-
+***/
+uInt STMolecules::addEntry( Vector<Double> restfreq, const Vector<String>& name,
+                            const Vector<String>& formattedname )
+{
+// How to handle this...?
+  Table result =
+    table_( nelements(table_.col("RESTFREQUENCY")) == restfreq.size() && 
+            all(table_.col("RESTFREQUENCY")== restfreq) );
+  uInt resultid = 0;
+  if ( result.nrow() > 0) {
+    ROScalarColumn<uInt> c(result, "ID");
+    c.get(0, resultid);
+  } else {
+    uInt rno = table_.nrow();
+    table_.addRow();
+    // get last assigned _id and increment
+    if ( rno > 0 ) {
+      idCol_.get(rno-1, resultid);
+      resultid++;
+    }
+    restfreqCol_.put(rno, restfreq);
+    nameCol_.put(rno, name);
+    formattednameCol_.put(rno, formattedname);
+    idCol_.put(rno, resultid);
+  }
+  return resultid;
+}
+
+
+
+
+/***
 void STMolecules::getEntry( Double& restfreq, String& name,
                             String& formattedname, uInt id ) const
@@ -104 +142 @@
   ROTableRow row(t);
   // get first row - there should only be one matching id
+
   const TableRecord& rec = row.get(0);
   restfreq = rec.asDouble("RESTFREQUENCY");
@@ -109 +148 @@
   formattedname = rec.asString("FORMATTEDNAME");
 }
+***/
+void STMolecules::getEntry( Vector<Double>& restfreq, Vector<String>& name,
+                            Vector<String>& formattedname, uInt id ) const
+{
+  Table t = table_(table_.col("ID") == Int(id) );
+  if (t.nrow() == 0 ) {
+    throw(AipsError("STMolecules::getEntry - id out of range"));
+  }
+  ROTableRow row(t);
+  // get first row - there should only be one matching id
+
+  const TableRecord& rec = row.get(0);
+  //restfreq = rec.asDouble("RESTFREQUENCY");
+  restfreq = rec.asArrayDouble("RESTFREQUENCY");
+  //name = rec.asString("NAME");
+  name = rec.asArrayString("NAME");
+  //formattedname = rec.asString("FORMATTEDNAME");
+  formattedname = rec.asArrayString("FORMATTEDNAME");
+}
 
 std::vector< double > asap::STMolecules::getRestFrequencies( ) const
 {
   std::vector<double> out;
+  //TableProxy itsTable(table_);
+  //Record rec;
   Vector<Double> rfs = restfreqCol_.getColumn();
   rfs.tovector(out);
+  //rec = itsTable.getVarColumn("RESTFREQUENCY", 0, 1, 1);
   return out;
 }
 
-double asap::STMolecules::getRestFrequency( uInt id ) const
-{
+std::vector< double > asap::STMolecules::getRestFrequency( uInt id ) const
+{
+  std::vector<double> out;
   Table t = table_(table_.col("ID") == Int(id) );
   if (t.nrow() == 0 ) {
@@ -126 +188 @@
   ROTableRow row(t);
   const TableRecord& rec = row.get(0);
-  return double(rec.asDouble("RESTFREQUENCY"));
-}
-
+  //return double(rec.asDouble("RESTFREQUENCY"));
+  Vector<Double> rfs = rec.asArrayDouble("RESTFREQUENCY");
+  rfs.tovector(out);
+  return out;
+}
+
+int asap::STMolecules::nrow() const
+{
+  return int(table_.nrow());
+}
 
 }//namespace

branches/alma/src/STMolecules.h

@@ -17 +17 @@
 #include <tables/Tables/Table.h>
 #include <tables/Tables/ScalarColumn.h>
+#include <tables/Tables/ArrayColumn.h>
+#include <casa/Arrays/Array.h>
 
 #include "STSubTable.h"
@@ -37 +39 @@
   STMolecules& operator=(const STMolecules& other);
 
+/***
   casa::uInt addEntry( casa::Double restfreq, const casa::String& name="",
                        const casa::String& formattedname="");
+***/
 
+  casa::uInt addEntry( casa::Vector<casa::Double> restfreq, const casa::Vector<casa::String>& name=casa::Vector<casa::String>(0),
+                       const casa::Vector<casa::String>& formattedname=casa::Vector<casa::String>(0));
+
+/***
   void getEntry( casa::Double& restfreq, casa::String& name,
                  casa::String& formattedname, casa::uInt id) const;
+***/
+  void getEntry( casa::Vector<casa::Double>& restfreq, casa::Vector<casa::String>& name,
+                 casa::Vector<casa::String>& formattedname, casa::uInt id) const;
 
   std::vector<double> getRestFrequencies() const;
-  double getRestFrequency( casa::uInt id ) const;
+  std::vector<double> getRestFrequency( casa::uInt id ) const;
   const casa::String& name() const { return name_; }
+  int nrow() const;
 
 private:
@@ -52 +64 @@
   //casa::Table table_;
   //casa::ScalarColumn<casa::uInt> freqidCol_;
-  casa::ScalarColumn<casa::Double> restfreqCol_;
-  casa::ScalarColumn<casa::String> nameCol_;
-  casa::ScalarColumn<casa::String> formattednameCol_; // e.g. latex
+  //casa::ScalarColumn<casa::Double> restfreqCol_;
+  casa::ArrayColumn<casa::Double> restfreqCol_;
+  //casa::ScalarColumn<casa::String> nameCol_;
+  casa::ArrayColumn<casa::String> nameCol_;
+  //casa::ScalarColumn<casa::String> formattednameCol_; // e.g. latex
+  casa::ArrayColumn<casa::String> formattednameCol_; // e.g. latex
 
 };

branches/alma/src/STWriter.cpp

@@ -165 +165 @@
       Table btable = beamit.table();
       // position only varies by beam
+      // No, we want to pointing data which varies by cycle!
       MDirection::ScalarColumn dirCol(btable, "DIRECTION");
       Vector<Double> direction = dirCol(0).getAngle("rad").getValue();
@@ -181 +182 @@
       while (!cycit.pastEnd() ) {
         Table ctable = cycit.table();
+        //MDirection::ScalarColumn dirCol(ctable, "DIRECTION");
+        //Vector<Double> direction = dirCol(0).getAngle("rad").getValue();
         TableIterator ifit(ctable, "IFNO");
         Int ifno = 1;
@@ -191 +194 @@
           ifno = rec.asuInt("IFNO")+1;
           uInt nchan = specCol(0).nelements();
-          Double cdelt,crval,crpix, restfreq;
+          //Double cdelt,crval,crpix, restfreq;
+          Double cdelt,crval,crpix;
+          Vector<Double> restfreq;
           Float focusAxi, focusTan, focusRot,
                 temperature, pressure, humidity, windSpeed, windAz;
           Float tmp0,tmp1,tmp2,tmp3,tmp4;
           Vector<Float> tcalval;
-          String stmp0,stmp1, tcalt;
+          //String stmp0,stmp1, tcalt;
+          String tcalt;
+          Vector<String> stmp0, stmp1;
           in->frequencies().getEntry(crpix,crval,cdelt, rec.asuInt("FREQ_ID"));
           in->focus().getEntry(focusAxi, focusTan, focusRot,
@@ -215 +222 @@
           Vector<Float> tsys = tsysFromTable(itable);
           // dummy data
+          //uInt npol; 
+          //if ( hdr.antennaname == "GBT" ) {
+          //  npol = nPolUsed;
+          //}
+          //else {   
+          //  npol = specs.ncolumn();
+          //}
           uInt npol = specs.ncolumn();
 
@@ -231 +245 @@
                                       refFreqNew, nchan*abs(cdelt), cdelt,
                                       restfreq,
-                                      tcal,
+                                      tcalval,
                                       tcalt,
                                       rec.asFloat("AZIMUTH"),
@@ -253 +267 @@
           }
           ++count;
-          //++ifno;
+          ++ifno;
           ++ifit;
         }
@@ -259 +273 @@
         ++cycit;
       }
-      //++beamno;
+      ++beamno;
       ++beamit;
     }
@@ -272 +286 @@
   if ( format_ == "MS2" ) {
     replacePtTab(table, filename);
-  }
+  } 
   return 0;
 }

branches/alma/src/Scantable.cpp

@@ -11 +11 @@
 //
 #include <map>
+#include <fstream>
 
 #include <casa/aips.h>
@@ -24 +25 @@
 #include <casa/Arrays/Vector.h>
 #include <casa/Arrays/VectorSTLIterator.h>
+#include <casa/Arrays/Slice.h>
 #include <casa/BasicMath/Math.h>
 #include <casa/BasicSL/Constants.h>
@@ -533 +535 @@
     return int(n);
   } else {
-    // take the first SCANNO,POLNO,BEAMNO,CYCLENO as nbeam shouldn't
-    // vary with these
+    // take the first SCANNO,POLNO,BEAMNO,CYCLENO as nbeam shouldn't vary with these
     Table t = table_(table_.col("IFNO") == ifno);
     if ( t.nrow() == 0 ) return 0;
@@ -786 +787 @@
   table_.keywordSet().get("FluxUnit", tmp);
   oss << setw(15) << "Flux Unit:" << tmp << endl;
-  Vector<Double> vec(moleculeTable_.getRestFrequencies());
+  //Vector<Double> vec(moleculeTable_.getRestFrequencies());
+  int nid = moleculeTable_.nrow();
+  Bool firstline = True;
   oss << setw(15) << "Rest Freqs:";
-  if (vec.nelements() > 0) {
-      oss << setprecision(10) << vec << " [Hz]" << endl;
-  } else {
-      oss << "none" << endl;
+  for (int i=0; i<nid; i++) {
+      Table t = table_(table_.col("MOLECULE_ID") == i);
+      if (t.nrow() >  0) {
+          Vector<Double> vec(moleculeTable_.getRestFrequency(i));
+          if (vec.nelements() > 0) {
+               if (firstline) { 
+                   oss << setprecision(10) << vec << " [Hz]" << endl;
+                   firstline=False;
+               }
+               else{ 
+                   oss << setw(15)<<" " << setprecision(10) << vec << " [Hz]" << endl;
+               } 
+          } else {
+              oss << "none" << endl;
+          }
+      }
   }
 
@@ -843 +858 @@
         oss << setw(10) << "";
         oss << setw(3) << std::right << irec.asuInt("IFNO") << std::left
-            << setw(2) << "" << frequencies().print(irec.asuInt("FREQ_ID"))
-            << endl;
+            << setw(2) << "" << frequencies().print(irec.asuInt("FREQ_ID"));
 
         ++iiter;
@@ -891 +905 @@
   const MDirection& md = getDirection(whichrow);
   const MEpoch& me = timeCol_(whichrow);
-  Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
+  //Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
+  Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
   SpectralCoordinate spc =
     freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
@@ -950 +965 @@
   const MDirection& md = getDirection(whichrow);
   const MEpoch& me = timeCol_(whichrow);
-  const Double& rf = mmolidCol_(whichrow);
+  //const Double& rf = mmolidCol_(whichrow);
+  const Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
   SpectralCoordinate spc =
     freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
@@ -967 +983 @@
 }
 
-void Scantable::setRestFrequencies( double rf, const std::string& name,
+/**
+void asap::Scantable::setRestFrequencies( double rf, const std::string& name,
                                           const std::string& unit )
+**/
+void Scantable::setRestFrequencies( vector<double> rf, const vector<std::string>& name,
+                                          const std::string& unit )
+
 {
   ///@todo lookup in line table to fill in name and formattedname
   Unit u(unit);
-  Quantum<Double> urf(rf, u);
-  uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), name, "");
+  //Quantum<Double> urf(rf, u);
+  Quantum<Vector<Double> >urf(rf, u);
+  Vector<String> formattedname(0); 
+  //cerr<<"Scantable::setRestFrequnecies="<<urf<<endl;
+  
+  //uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), name, "");
+  uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), mathutil::toVectorString(name), formattedname);
   TableVector<uInt> tabvec(table_, "MOLECULE_ID");
   tabvec = id;
 }
 
-void Scantable::setRestFrequencies( const std::string& name )
+/**
+void asap::Scantable::setRestFrequencies( const std::string& name )
 {
   throw(AipsError("setRestFrequencies( const std::string& name ) NYI"));
+  ///@todo implement
+}
+**/
+void Scantable::setRestFrequencies( const vector<std::string>& name )
+{
+  throw(AipsError("setRestFrequencies( const vector<std::string>& name ) NYI"));
   ///@todo implement
 }
@@ -1117 +1150 @@
 }
 
-}
- //namespace asap
+void asap::Scantable::reshapeSpectrum( int nmin, int nmax )
+  throw( casa::AipsError )
+{
+  // assumed that all rows have same nChan
+  Vector<Float> arr = specCol_( 0 ) ;
+  int nChan = arr.nelements() ;
+  
+  // if nmin < 0 or nmax < 0, nothing to do
+  if (  nmin < 0 ) {
+    throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
+    }
+  if (  nmax < 0  ) {
+    throw( casa::indexError<int>( nmax, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
+  }
+  
+  // if nmin > nmax, exchange values
+  if ( nmin > nmax ) {
+    int tmp = nmax ;
+    nmax = nmin ;
+    nmin = tmp ;
+    cout << "Swap values. Applied range is [" 
+         << nmin << ", " << nmax << "]" << endl ;
+  }
+  
+  // if nmin exceeds nChan, nothing to do
+  if ( nmin >= nChan ) {
+    throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Specified minimum exceeds nChan." ) ) ;
+  }
+  
+  // if nmax exceeds nChan, reset nmax to nChan
+  if ( nmax >= nChan ) {
+    if ( nmin == 0 ) {
+      // nothing to do
+      cout << "Whole range is selected. Nothing to do." << endl ;
+      return ;
+    }
+    else {
+      cout << "Specified maximum exceeds nChan. Applied range is ["
+           << nmin << ", " << nChan-1 << "]." << endl ;
+      nmax = nChan - 1 ;
+    }
+  }
+  
+  // reshape specCol_ and flagCol_
+  for ( int irow = 0 ; irow < nrow() ; irow++ ) {
+    reshapeSpectrum( nmin, nmax, irow ) ;
+  }
+
+  // update FREQUENCIES subtable
+  Double refpix ;
+  Double refval ;
+  Double increment ;
+  int freqnrow = freqTable_.table().nrow() ;
+  Vector<uInt> oldId( freqnrow ) ;
+  Vector<uInt> newId( freqnrow ) ;
+  for ( int irow = 0 ; irow < freqnrow ; irow++ ) {
+    freqTable_.getEntry( refpix, refval, increment, irow ) ;
+    /***
+     * need to shift refpix to nmin
+     * note that channel nmin in old index will be channel 0 in new one
+     ***/
+    refval = refval - ( refpix - nmin ) * increment ;
+    refpix = 0 ;  
+    freqTable_.setEntry( refpix, refval, increment, irow ) ;
+  }
+  
+  // update nchan
+  int newsize = nmax - nmin + 1 ;
+  table_.rwKeywordSet().define( "nChan", newsize ) ;
+  
+  // update bandwidth
+  // assumed all spectra in the scantable have same bandwidth
+  table_.rwKeywordSet().define( "Bandwidth", increment * newsize ) ;
+  
+  return ;
+}
+
+void asap::Scantable::reshapeSpectrum( int nmin, int nmax, int irow ) 
+{
+  // reshape specCol_ and flagCol_
+  Vector<Float> oldspec = specCol_( irow ) ;
+  Vector<uChar> oldflag = flagsCol_( irow ) ;
+  uInt newsize = nmax - nmin + 1 ;
+  specCol_.put( irow, oldspec( Slice( nmin, newsize, 1 ) ) ) ;
+  flagsCol_.put( irow, oldflag( Slice( nmin, newsize, 1 ) ) ) ;
+
+  return ;
+}
+
+void asap::Scantable::regridChannel( int nChan, double dnu )
+{
+  cout << "Regrid abcissa with channel number " << nChan << " and spectral resoultion " << dnu << "Hz." << endl ;
+  // assumed that all rows have same nChan
+  Vector<Float> arr = specCol_( 0 ) ;
+  int oldsize = arr.nelements() ;
+
+  // if oldsize == nChan, nothing to do
+  if ( oldsize == nChan ) {
+    cout << "Specified channel number is same as current one. Nothing to do." << endl ;
+    return ;
+  }
+
+  // if oldChan < nChan, unphysical operation
+  if ( oldsize < nChan ) {
+    cout << "Unphysical operation. Nothing to do." << endl ;
+    return ;
+  }
+
+  // change channel number for specCol_ and flagCol_
+  Vector<Float> newspec( nChan, 0 ) ;
+  Vector<uChar> newflag( nChan, false ) ;
+  vector<string> coordinfo = getCoordInfo() ;
+  string oldinfo = coordinfo[0] ;
+  coordinfo[0] = "Hz" ;
+  setCoordInfo( coordinfo ) ;
+  for ( int irow = 0 ; irow < nrow() ; irow++ ) {
+    regridChannel( nChan, dnu, irow ) ;
+  } 
+  coordinfo[0] = oldinfo ;
+  setCoordInfo( coordinfo ) ;
+
+
+  // NOTE: this method does not update metadata such as 
+  //       FREQUENCIES subtable, nChan, Bandwidth, etc.
+
+  return ; 
+}
+
+void asap::Scantable::regridChannel( int nChan, double dnu, int irow ) 
+{
+  // logging
+  //ofstream ofs( "average.log", std::ios::out | std::ios::app ) ;
+  //ofs << "IFNO = " << getIF( irow ) << " irow = " << irow << endl ;
+
+  Vector<Float> oldspec = specCol_( irow ) ;
+  Vector<uChar> oldflag = flagsCol_( irow ) ;
+  Vector<Float> newspec( nChan, 0 ) ;
+  Vector<uChar> newflag( nChan, false ) ;
+  
+  // regrid
+  vector<double> abcissa = getAbcissa( irow ) ;
+  int oldsize = abcissa.size() ;
+  double olddnu = abcissa[1] - abcissa[0] ;
+  int refChan = 0 ;
+  double frac = 0.0 ;
+  double wedge = 0.0 ;
+  double pile = 0.0 ;
+  double wsum = 0.0 ;
+  /*** 
+   * ichan = 0
+   ***/
+  //ofs << "olddnu = " << olddnu << ", dnu = " << dnu << endl ;
+  pile += dnu ;
+  wedge = olddnu * ( refChan + 1 ) ;
+  while ( wedge < pile ) {
+    newspec[0] += olddnu * oldspec[refChan] ;
+    newflag[0] = newflag[0] || oldflag[refChan] ;
+    //ofs << "channel " << refChan << " is included in new channel 0" << endl ;
+    refChan++ ;
+    wedge += olddnu ;
+    wsum += olddnu ;
+    //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
+  }
+  frac = ( wedge - pile ) / olddnu ;
+  wsum += ( 1.0 - frac ) * olddnu ;
+  newspec[0] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
+  newflag[0] = newflag[0] || oldflag[refChan] ;
+  //ofs << "channel " << refChan << " is partly included in new channel 0" << " with fraction of " << ( 1.0 - frac ) << endl ;
+  //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
+  newspec[0] /= wsum ;
+  //ofs << "newspec[0] = " << newspec[0] << endl ;
+  //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
+
+  /***
+   * ichan = 1 - nChan-2
+   ***/
+  for ( int ichan = 1 ; ichan < nChan - 1 ; ichan++ ) {
+    pile += dnu ;
+    newspec[ichan] += frac * olddnu * oldspec[refChan] ;
+    newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
+    //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << frac << endl ;
+    refChan++ ;
+    wedge += olddnu ;
+    wsum = frac * olddnu ;
+    //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
+    while ( wedge < pile ) {
+      newspec[ichan] += olddnu * oldspec[refChan] ;
+      newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
+      //ofs << "channel " << refChan << " is included in new channel " << ichan << endl ;
+      refChan++ ;
+      wedge += olddnu ;
+      wsum += olddnu ;
+      //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
+    }
+    frac = ( wedge - pile ) / olddnu ;
+    wsum += ( 1.0 - frac ) * olddnu ;
+    newspec[ichan] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
+    newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
+    //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << ( 1.0 - frac ) << endl ;
+    //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
+    //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
+    newspec[ichan] /= wsum ;
+    //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << endl ;
+  }
+
+  /***
+   * ichan = nChan-1
+   ***/
+  // NOTE: Assumed that all spectra have the same bandwidth
+  pile += dnu ; 
+  newspec[nChan-1] += frac * olddnu * oldspec[refChan] ;
+  newflag[nChan-1] = newflag[nChan-1] || oldflag[refChan] ;
+  //ofs << "channel " << refChan << " is partly included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
+  refChan++ ;
+  wedge += olddnu ;
+  wsum = frac * olddnu ;
+  //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
+  for ( int jchan = refChan ; jchan < oldsize ; jchan++ ) {
+    newspec[nChan-1] += olddnu * oldspec[jchan] ;
+    newflag[nChan-1] = newflag[nChan-1] || oldflag[jchan] ;
+    wsum += olddnu ;
+    //ofs << "channel " << jchan << " is included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
+    //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
+  }
+  //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
+  //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
+  newspec[nChan-1] /= wsum ;
+  //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << endl ;
+  
+  specCol_.put( irow, newspec ) ;
+  flagsCol_.put( irow, newflag ) ;
+
+  // ofs.close() ;
+
+  return ;
+}
+
+}
+//namespace asap

branches/alma/src/Scantable.h

@@ -22 +22 @@
 #include <casa/Utilities/CountedPtr.h>
 
+#include <casa/Exceptions/Error.h>
+
 #include <coordinates/Coordinates/SpectralCoordinate.h>
 
@@ -29 +31 @@
 
 #include <measures/TableMeasures/ScalarMeasColumn.h>
+
+#include <casa/Arrays/Vector.h>
+#include <casa/Quanta/Quantum.h>
 
 #include "Logger.h"
@@ -167 +172 @@
 
         /**
-         * get the direction type as a string, e.g. "J2000"
+         * get the direction as a string
          * @param[in] whichrow the row number
          * return the direction string
@@ -277 +282 @@
   std::vector<uint> getScanNos() { return getNumbers(scanCol_); }
   int getScan(int whichrow) const { return scanCol_(whichrow); }
+
+  //TT addition
+  std::vector<uint> getMolNos() {return getNumbers(mmolidCol_); }
 
   /**
@@ -352 +360 @@
   std::vector<double> getRestFrequencies() const
     { return moleculeTable_.getRestFrequencies(); }
-
+  std::vector<double> getRestFrequency(int id) const
+    { return moleculeTable_.getRestFrequency(id); }
+
+  /**
   void setRestFrequencies(double rf, const std::string& name = "",
                           const std::string& = "Hz");
-  void setRestFrequencies(const std::string& name);
+  **/
+  // Modified by Takeshi Nakazato 05/09/2008
+  /***
+  void setRestFrequencies(vector<double> rf, const vector<std::string>& name = "",
+                          const std::string& = "Hz");
+  ***/
+  void setRestFrequencies(vector<double> rf,
+                          const vector<std::string>& name = vector<std::string>(1,""),
+                          const std::string& = "Hz");
+
+  //void setRestFrequencies(const std::string& name);
+  void setRestFrequencies(const vector<std::string>& name);
 
   void shift(int npix);
@@ -393 +415 @@
    * against the information found in GBT_GO table for
    * scan number orders to get correct pairs.
-   * @param[in] scan list
-   * @return status
+   * @param[in] scan list 
+   * @return status 
    */
   int checkScanInfo(const std::vector<int>& scanlist) const;
 
   /**
-   * Get the direction as a vector, for a specific row
+   * Get the direction as a vector, for a specific row  
    * @param[in] whichrow the row numbyyer
-   * @return the direction in a vector
+   * @return the direction in a vector 
    */
   std::vector<double> getDirectionVector(int whichrow) const;
 
+  /**
+   * Reshape spectrum
+   * @param[in] nmin, nmax minimum and maximum channel
+   * @param[in] irow       row number
+   * 
+   * 30/07/2008 Takeshi Nakazato  
+   **/
+  void reshapeSpectrum( int nmin, int nmax ) throw( casa::AipsError );
+  void reshapeSpectrum( int nmin, int nmax, int irow ) ;
+
+  /**
+   * Change channel number under fixed bandwidth
+   * @param[in] nchan, dnu new channel number and spectral resolution
+   * @param[in] irow       row number
+   *
+   * 27/08/2008 Takeshi Nakazato
+   **/
+  void regridChannel( int nchan, double dnu ) ;
+  void regridChannel( int nchan, double dnu, int irow ) ;
+ 
 private:
 

branches/alma/src/ScantableWrapper.h

@@ -133 +133 @@
   std::vector<uint> getScanNos() { return table_->getScanNos(); }
   int getScan(int whichrow) const {return table_->getScan(whichrow);}
+  std::vector<uint> getMolNos() { return table_->getMolNos();}
 
   STSelector getSelection() const { return table_->getSelection(); }
@@ -157 +158 @@
   { table_->shift(npix); }
 
+/**
+  commented out by TT
   void setRestFrequencies(double rf, const std::string& name,
                           const std::string& unit)
     { table_->setRestFrequencies(rf, name, unit); }
+**/
+  void setRestFrequencies(vector<double> rf, const vector<std::string>& name,
+                          const std::string& unit)
+    { table_->setRestFrequencies(rf, name, unit); }
+
 /*
   void setRestFrequencies(const std::string& name) {
@@ -166 +174 @@
 */
 
+/*
   std::vector<double> getRestFrequencies() const
     { return table_->getRestFrequencies(); }
+*/
+  std::vector<double> getRestFrequency(int id) const
+    { return table_->getRestFrequency(id); }
 
   void setCoordInfo(std::vector<string> theinfo) {
@@ -208 +220 @@
   int checkScanInfo(const vector<int>& scanlist) const
     { return table_->checkScanInfo(scanlist); }
-
+ 
   std::vector<double>  getDirectionVector(int whichrow) const
     { return table_->getDirectionVector(whichrow); }
+
+  void reshapeSpectrum( int nmin, int nmax )
+  { table_->reshapeSpectrum( nmin, nmax ); }
 
 private:

branches/alma/src/Templates.cpp

@@ -41 +41 @@
 template class casa::PtrRep<asap::Scantable>;
 template class casa::SimpleCountedConstPtr<asap::Scantable>;
+template class casa::SimpleCountedConstPtr<asap::STPol>;
 template class casa::SimpleCountedPtr<asap::Scantable>;
 template class casa::SimpleCountedPtr<asap::STPol>;

branches/alma/src/python_STMath.cpp

@@ -73 +73 @@
         .def("_mx_extract", &STMathWrapper::mxExtract)
         .def("_lag_flag", &STMathWrapper::lagFlag)
+        // testing average spectra with different channel/resolution
+        .def("_new_average", &STMathWrapper::new_average)
           ;
     };

branches/alma/src/python_Scantable.cpp

@@ -58 +58 @@
     .def("getscannos", &ScantableWrapper::getScanNos)
     .def("getcycle", &ScantableWrapper::getCycle)
+    .def("getmolnos", &ScantableWrapper::getMolNos)
     .def("nif", &ScantableWrapper::nif,
          (boost::python::arg("scanno")=-1) )
@@ -107 +108 @@
     .def("_summary",  &ScantableWrapper::summary,
          (boost::python::arg("verbose")=true) )
-    .def("_getrestfreqs",  &ScantableWrapper::getRestFrequencies)
+    //.def("_getrestfreqs",  &ScantableWrapper::getRestFrequencies)
+    .def("_getrestfreqs",  &ScantableWrapper::getRestFrequency)
     .def("_setrestfreqs",  &ScantableWrapper::setRestFrequencies)
     .def("shift_refpix", &ScantableWrapper::shift)
@@ -123 +125 @@
     .def("_setsourcetype", &ScantableWrapper::setSourceType)
     .def("_getdirectionvec", &ScantableWrapper::getDirectionVector)
+    .def("_reshape", &ScantableWrapper::reshapeSpectrum, 
+         (boost::python::arg("nmin")=-1, 
+          boost::python::arg("nmax")=-1) )
   ;
 };