Changeset 2443 for branches/hpc33

Timestamp:

03/23/12 18:20:50 (13 years ago)

Author:

Takeshi Nakazato

Message:

New Development: No

JIRA Issue: No

Ready for Test: No

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs:

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Test tuning ALMA calibration code.

Location:

branches/hpc33/src

Files:

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

branches/hpc33/src/STMath.cpp

@@ -4450 +4450 @@
 }
 
+vector<float> STMath::getSpectrumFromTime( double reftime, 
+                                           CountedPtr<Scantable>& s,
+                                           string mode ) 
+{
+  LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ;
+  vector<float> sp ;
+
+  if ( s->nrow() == 0 ) {
+    os << LogIO::SEVERE << "No spectra in the input scantable. Return empty spectrum." << LogIO::POST ;
+    return sp ;
+  }
+  else if ( s->nrow() == 1 ) {
+    //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ;
+    return s->getSpectrum( 0 ) ; 
+  }
+  else {
+    ROScalarColumn<Double> timeCol( s->table(), "TIME" ) ;
+    Vector<Double> timeVec = timeCol.getColumn() ;
+    vector<int> idx = getRowIdFromTime( reftime, timeVec ) ;
+    if ( mode == "before" ) {
+      int id = -1 ;
+      if ( idx[0] != -1 ) {
+        id = idx[0] ;
+      }
+      else if ( idx[1] != -1 ) {
+        os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
+        id = idx[1] ;
+      }
+      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
+      sp = s->getSpectrum( id ) ;
+    }
+    else if ( mode == "after" ) {
+      int id = -1 ;
+      if ( idx[1] != -1 ) {
+        id = idx[1] ;
+      }
+      else if ( idx[0] != -1 ) {
+        os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
+        id = idx[1] ;
+      }
+      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
+      sp = s->getSpectrum( id ) ;
+    }
+    else if ( mode == "nearest" ) {
+      int id = -1 ;
+      if ( idx[0] == -1 ) {
+        id = idx[1] ;
+      }
+      else if ( idx[1] == -1 ) {
+        id = idx[0] ;
+      }
+      else if ( idx[0] == idx[1] ) {
+        id = idx[0] ;
+      }
+      else {
+        //double t0 = getMJD( s->getTime( idx[0] ) ) ;
+        //double t1 = getMJD( s->getTime( idx[1] ) ) ;
+//         double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
+//         double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
+        double t0 = timeVec[idx[0]] ;
+        double t1 = timeVec[idx[1]] ;
+//         cout << "t0-t0c=" << t0-t0c << endl ;
+//         cout << "t1-t1c=" << t1-t1c << endl ;
+//         double tref = getMJD( reftime ) ;
+        double tref = reftime ;
+        if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
+          id = idx[1] ;
+        }
+        else {
+          id = idx[0] ;
+        }
+      }
+      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
+      sp = s->getSpectrum( id ) ;      
+    }
+    else if ( mode == "linear" ) {
+      if ( idx[0] == -1 ) {
+        // use after
+        os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
+        int id = idx[1] ;
+        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
+        sp = s->getSpectrum( id ) ;
+      }
+      else if ( idx[1] == -1 ) {
+        // use before
+        os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
+        int id = idx[0] ;
+        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
+        sp = s->getSpectrum( id ) ;
+      }
+      else if ( idx[0] == idx[1] ) {
+        // use before
+        //os << "No need to interporate." << LogIO::POST ;
+        int id = idx[0] ;
+        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
+        sp = s->getSpectrum( id ) ;
+      }
+      else {
+        // do interpolation
+        //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
+        //double t0 = getMJD( s->getTime( idx[0] ) ) ;
+        //double t1 = getMJD( s->getTime( idx[1] ) ) ;
+//         double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
+//         double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
+        double t0 = timeVec[idx[0]] ;
+        double t1 = timeVec[idx[1]] ;
+//         cout << "t0-t0c=" << t0-t0c << endl ;
+//         cout << "t1-t1c=" << t1-t1c << endl ;
+//         double tref = getMJD( reftime ) ;
+        double tref = reftime ;
+        vector<float> sp0 = s->getSpectrum( idx[0] ) ;
+        vector<float> sp1 = s->getSpectrum( idx[1] ) ;
+        for ( unsigned int i = 0 ; i < sp0.size() ; i++ ) {
+          float v = ( sp1[i] - sp0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + sp0[i] ;
+          sp.push_back( v ) ;
+        }
+      }
+    }
+    else {
+      os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
+    }
+    return sp ;
+  }
+}
+
 double STMath::getMJD( string strtime ) 
 {
@@ -4508 +4633 @@
   v.push_back( just_before ) ;
   v.push_back( just_after ) ;
+
+  return v ;
+}
+
+vector<int> STMath::getRowIdFromTime( double reftime, Vector<Double> &t )
+{
+//   double reft = reftime ;
+  double dtmin = 1.0e100 ;
+  double dtmax = -1.0e100 ;
+//   vector<double> dt ;
+  int just_before = -1 ;
+  int just_after = -1 ;
+  //cout << setprecision(24) << reft << endl ;
+//   ROScalarColumn<Double> timeCol( s->table(), "TIME" ) ;
+//   for ( int i = 0 ; i < s->nrow() ; i++ ) {
+//     cout << setprecision(24) << timeCol(i) << endl ;
+//     //dt.push_back( getMJD( s->getTime( i ) ) - reft ) ; 
+//     dt.push_back( timeCol(i) - reft ) ; 
+//   }
+  Vector<Double> dt = t - reftime ;
+  for ( unsigned int i = 0 ; i < dt.size() ; i++ ) {
+    if ( dt[i] > 0.0 ) {
+      // after reftime
+      if ( dt[i] < dtmin ) {
+        just_after = i ;
+        dtmin = dt[i] ;
+      }
+    }
+    else if ( dt[i] < 0.0 ) {
+      // before reftime
+      if ( dt[i] > dtmax ) {
+        just_before = i ;
+        dtmax = dt[i] ;
+      }
+    }
+    else {
+      // just a reftime
+      just_before = i ;
+      just_after = i ;
+      dtmax = 0 ;
+      dtmin = 0 ;
+      break ;
+    }
+  }
+
+  vector<int> v(2) ;
+  v[0] = just_before ;
+  v[1] = just_after ;
 
   return v ;
@@ -4833 +5006 @@
                                             int index ) 
 {
-  string reftime = on->getTime( index ) ;
+//   string reftime = on->getTime( index ) ;
+  ROTableColumn timeCol( on->table(), "TIME" ) ;
+  double reftime = timeCol.asdouble(index) ;
   vector<int> ii( 1, on->getIF( index ) ) ;
   vector<int> ib( 1, on->getBeam( index ) ) ;

branches/hpc33/src/STMath.h

@@ -403 +403 @@
 
   vector<float> getSpectrumFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode = "before" ) ;
+  vector<float> getSpectrumFromTime( double reftime, casa::CountedPtr<Scantable>& s, string mode = "before" ) ;
   vector<float> getTcalFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ;
   vector<float> getTsysFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ;
   vector<int> getRowIdFromTime( string reftime, casa::CountedPtr<Scantable>& s ) ;
+  vector<int> getRowIdFromTime( double reftime, casa::Vector<casa::Double>& t ) ;
 
   // Chopper-Wheel type calibration