Changeset 1610 for branches

Timestamp:

08/03/09 11:34:53 (15 years ago)

Author:

Takeshi Nakazato

Message:

New Development: No

JIRA Issue: No

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Update of spectral regrigging algorithm.

File:

• branches/alma/src/Scantable.cpp (modified) (1 diff)

branches/alma/src/Scantable.cpp

@@ -1308 +1308 @@
   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 ;
+  //int refChan = 0 ;
+  //double frac = 0.0 ;
+  //double wedge = 0.0 ;
+  //double pile = 0.0 ;
+  int ichan = 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 ;
+  Vector<Float> z( nChan ) ;
+  z[0] = abcissa[0] - 0.5 * olddnu + 0.5 * dnu ;
+  for ( int ii = 1 ; ii < nChan ; ii++ )
+    z[ii] = z[ii-1] + dnu ;
+  Vector<Float> zi( nChan+1 ) ;
+  Vector<Float> yi( oldsize + 1 ) ;
+  zi[0] = z[0] - 0.5 * dnu ;
+  zi[1] = z[0] + 0.5 * dnu ;
+  for ( int ii = 2 ; ii < nChan ; ii++ ) 
+    zi[ii] = zi[ii-1] + dnu ;
+  zi[nChan] = z[nChan-1] + 0.5 * dnu ;
+  yi[0] = abcissa[0] - 0.5 * olddnu ;
+  yi[1] = abcissa[1] + 0.5 * olddnu ;
+  for ( int ii = 2 ; ii < oldsize ; ii++ )
+    yi[ii] = abcissa[ii-1] + olddnu ;
+  yi[oldsize] = abcissa[oldsize-1] + 0.5 * olddnu ; 
+  if ( dnu > 0.0 ) {
+    for ( int ii = 0 ; ii < nChan ; ii++ ) {
+      double zl = zi[ii] ;
+      double zr = zi[ii+1] ;
+      for ( int j = ichan ; j < oldsize ; j++ ) {
+        double yl = yi[j] ;
+        double yr = yi[j+1] ;
+        if ( yl <= zl ) {
+          if ( yr <= zl ) {
+            continue ;
+          }
+          else if ( yr <= zr ) {
+            newspec[ii] += oldspec[j] * ( yr - zl ) ;
+            newflag[ii] = newflag[ii] || oldflag[j] ;
+            wsum += ( yr - zl ) ;
+          }
+          else {
+            newspec[ii] += oldspec[j] * dnu ;
+            newflag[ii] = newflag[ii] || oldflag[j] ;
+            wsum += dnu ;
+            ichan = j ;
+            break ;
+          }
+        }
+        else if ( yl < zr ) {
+          if ( yr <= zr ) {
+              newspec[ii] += oldspec[j] * ( yr - yl ) ;
+              newflag[ii] = newflag[ii] || oldflag[j] ;
+              wsum += ( yr - yl ) ;
+          }
+          else {
+            newspec[ii] += oldspec[j] * ( zr - yl ) ;
+            newflag[ii] = newflag[ii] || oldflag[j] ;
+            wsum += ( zr - yl ) ;
+            ichan = j ;
+            break ;
+          }
+        }
+        else {
+          ichan = j - 1 ;
+          break ;
+        }
+      }
+      newspec[ii] /= wsum ;
+      wsum = 0.0 ;
+    }
+  }
+  else if ( dnu < 0.0 ) {
+    for ( int ii = 0 ; ii < nChan ; ii++ ) {
+      double zl = zi[ii] ;
+      double zr = zi[ii+1] ;
+      for ( int j = ichan ; j < oldsize ; j++ ) {
+        double yl = yi[j] ;
+        double yr = yi[j+1] ;
+        if ( yl >= zl ) {
+          if ( yr >= zl ) {
+            continue ;
+          }
+          else if ( yr >= zr ) {
+            newspec[ii] += oldspec[j] * abs( yr - zl ) ;
+            newflag[ii] = newflag[ii] || oldflag[j] ;
+            wsum += abs( yr - zl ) ;
+          }
+          else {
+            newspec[ii] += oldspec[j] * abs( dnu ) ;
+            newflag[ii] = newflag[ii] || oldflag[j] ;
+            wsum += abs( dnu ) ;
+            ichan = j ;
+            break ;
+          }
+        }
+        else if ( yl > zr ) {
+          if ( yr >= zr ) {
+            newspec[ii] += oldspec[j] * abs( yr - yl ) ;
+            newflag[ii] = newflag[ii] || oldflag[j] ;
+            wsum += abs( yr - yl ) ;
+          }
+          else {
+            newspec[ii] += oldspec[j] * abs( zr - yl ) ;
+            newflag[ii] = newflag[ii] || oldflag[j] ;
+            wsum += abs( zr - yl ) ;
+            ichan = j ;
+            break ;
+          }
+        }
+        else {
+          ichan = j - 1 ;
+          break ;
+        }
+      }
+      newspec[ii] /= wsum ;
+      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() ;
+//   specCol_.put( irow, newspec ) ;
+//   flagsCol_.put( irow, newflag ) ;
+
+//   // ofs.close() ;
+
 
   return ;