Changeset 1446 for branches/alma/src/Scantable.cpp

Timestamp:

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

Author:

TakTsutsumi

Message:

Merged recent updates (since 2007) from nrao-asap

File:

• branches/alma/src/Scantable.cpp (modified) (9 diffs)

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