Changeset 1609 for branches/alma

Timestamp:

07/29/09 19:00:30 (15 years ago)

Author:

Takeshi Nakazato

Message:

New Development: No

JIRA Issue: Yes CAS-1448

Ready to Release: No

Interface Changes: Yes/No?

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

I have defined a function for folding frequency-switch data.
It is still experimental. Test data is needed.

Location:

branches/alma/src

Files:

• STMath.cpp (modified) (5 diffs)
• STMath.h (modified) (1 diff)

branches/alma/src/STMath.cpp

@@ -133 +133 @@
   uInt outrowCount = 0;
   TableIterator iter(baset, cols);
-  int count = 0 ;
+//   int count = 0 ;
   while (!iter.pastEnd()) {
     Table subt = iter.table();
@@ -1192 +1192 @@
   Table& tabout2=out2->table();
   ROScalarColumn<uInt> freqidCol1(tabout1, "FREQ_ID");
-  ROScalarColumn<uInt> freqidCol2(tabout2, "FREQ_ID");
+  ScalarColumn<uInt> freqidCol2(tabout2, "FREQ_ID");
   ROArrayColumn<Float> specCol(tabout2, "SPECTRA");
   Vector<Float> spec; specCol.get(0, spec);
@@ -1201 +1201 @@
   out1->frequencies().getEntry(rp1, rv1, inc1, freqid1);
   out2->frequencies().getEntry(rp2, rv2, inc2, freqid2);
+  //cerr << out1->frequencies().table().nrow() << " " << out2->frequencies().table().nrow() << endl ;
   if (rp1==rp2) {
     Double foffset = rv1 - rv2;
@@ -1209 +1210 @@
     }
   }
- 
+
   if (nofold) {
     std::vector< CountedPtr< Scantable > > tabs; 
@@ -1217 +1218 @@
   }
   else { //folding is not implemented yet
-    out = out1;
+    //out = out1;
+    Int choffset = static_cast<Int>((rv1-rv2)/inc2) ;
+    out = dofold( out1, out2, choffset ) ;
   }
     
   return out;
 }
+
+CountedPtr<Scantable> STMath::dofold( const CountedPtr<Scantable> &sig,
+                                      const CountedPtr<Scantable> &ref,
+                                      Int choffset )
+{
+  // output scantable
+  CountedPtr<Scantable> out = getScantable( sig, false ) ;
+
+  // get column
+  ROArrayColumn<Float> specCol1( sig->table(), "SPECTRA" ) ;
+  ROArrayColumn<Float> specCol2( ref->table(), "SPECTRA" ) ;
+  ROArrayColumn<Float> tsysCol1( sig->table(), "TSYS" ) ;
+  ROArrayColumn<Float> tsysCol2( ref->table(), "TSYS" ) ;
+  ROArrayColumn<uChar> flagCol1( sig->table(), "FLAGTRA" ) ;
+  ROArrayColumn<uChar> flagCol2( ref->table(), "FLAGTRA" ) ;
+  ROScalarColumn<Double> mjdCol1( sig->table(), "TIME" ) ;
+  ROScalarColumn<Double> mjdCol2( ref->table(), "TIME" ) ;
+  ROScalarColumn<Double> intervalCol1( sig->table(), "INTERVAL" ) ;
+  ROScalarColumn<Double> intervalCol2( ref->table(), "INTERVAL" ) ;
+
+  // check
+  if ( choffset == 0 ) {
+    cerr << "channel offset is zero, no folding" << endl ;
+    return out ;
+  }
+  int nchan = ref->nchan() ;
+  if ( abs(choffset) >= nchan ) {
+    cerr << "out-band frequency switching, no folding" << endl ;
+    return out ;
+  }
+
+  // attach column for output scantable
+  ArrayColumn<Float> specColOut( out->table(), "SPECTRA" ) ;
+  ArrayColumn<uChar> flagColOut( out->table(), "FLAGTRA" ) ;
+  ArrayColumn<Float> tsysColOut( out->table(), "TSYS" ) ;
+  ScalarColumn<Double> mjdColOut( out->table(), "TIME" ) ;
+  ScalarColumn<Double> intervalColOut( out->table(), "INTERVAL" ) ;
+
+  // for each row
+  // assume that the data order are same between sig and ref
+  RowAccumulator acc( asap::TINTSYS ) ; 
+  for ( int i = 0 ; i < sig->nrow() ; i++ ) {
+    // get values
+    Vector<Float> spsig ;
+    specCol1.get( i, spsig ) ;
+    Vector<Float> spref ;
+    specCol2.get( i, spref ) ;
+    Vector<Float> tsyssig ;
+    tsysCol1.get( i, tsyssig ) ;
+    Vector<Float> tsysref ;
+    tsysCol2.get( i, tsysref ) ;
+    Vector<uChar> flagsig ;
+    flagCol1.get( i, flagsig ) ;
+    Vector<uChar> flagref ;
+    flagCol2.get( i, flagref ) ;
+    Double timesig ;
+    mjdCol1.get( i, timesig ) ;
+    Double timeref ;
+    mjdCol2.get( i, timeref ) ;
+    Double intsig ;
+    intervalCol1.get( i, intsig ) ;
+    Double intref ;
+    intervalCol2.get( i, intref ) ;
+
+    // shift reference spectra
+    int refchan = spref.nelements() ;
+    if ( choffset > 0 ) {
+      for ( int j = 0 ; j < refchan-choffset ; j++ ) {
+        spref[j] = spref[j+choffset] ; 
+        tsysref[j] = tsysref[j+choffset] ;
+        flagref[j] = flagref[j+choffset] ;
+      }
+      for ( int j = refchan-choffset ; j < refchan ; j++ ) {
+        spref[j] = spref[j-refchan+choffset] ;
+        tsysref[j] = tsysref[j-refchan+choffset] ;
+        flagref[j] = flagref[j-refchan+choffset] ;
+      }
+    }
+    else {
+      for ( int j = 0 ; j < abs(choffset) ; j++ ) {
+        spref[j] = spref[refchan+choffset+j] ;
+        tsysref[j] = tsysref[refchan+choffset+j] ;
+        flagref[j] = flagref[refchan+choffset+j] ;
+      }
+      for ( int j = abs(choffset) ; j < refchan ; j++ ) {
+        spref[j] = spref[j+choffset] ;
+        tsysref[j] = tsysref[j+choffset] ;
+        flagref[j] = flagref[j+choffset] ;
+      }
+    }
+
+    // folding
+    acc.add( spsig, !flagsig, tsyssig, intsig, timesig ) ;
+    acc.add( spref, !flagref, tsysref, intref, timeref ) ;
+    
+    // put result
+    specColOut.put( i, acc.getSpectrum() ) ;
+    const Vector<Bool> &msk = acc.getMask() ;
+    Vector<uChar> flg( msk.shape() ) ;
+    convertArray( flg, !msk ) ;
+    flagColOut.put( i, flg ) ;
+    tsysColOut.put( i, acc.getTsys() ) ;
+    intervalColOut.put( i, acc.getInterval() ) ;
+    mjdColOut.put( i, acc.getTime() ) ;
+
+    acc.reset() ;
+  }
+
+  return out ;
+}
+
 
 CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )

branches/alma/src/STMath.h

@@ -199 +199 @@
                                     casa::Float tcal=0.0 );
 
+  /**
+   * Folding frequency-switch data
+   * @param sig
+   * @param ref
+   * @param choffset
+   **/
+  casa::CountedPtr<Scantable> dofold( const casa::CountedPtr<Scantable> &sig,
+                                      const casa::CountedPtr<Scantable> &ref,
+                                      casa::Int choffset );
 
   casa::CountedPtr<Scantable>