Changeset 2462


Ignore:
Timestamp:
04/12/12 15:17:35 (12 years ago)
Author:
Kana Sugimoto
Message:

New Development: No (a bug fix)

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs: sdsmooth unit tests (test07-09)

Put in Release Notes: No

Module(s):

Description: fixed a bug in Scantable::regridChannel.

The function now works properly for spectra with negative frequency increments.


File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/Scantable.cpp

    r2437 r2462  
    7070#include "STUpgrade.h"
    7171#include "Scantable.h"
     72
     73#define debug 1
    7274
    7375using namespace casa;
    … …  
    19181920       ***/
    19191921      //refval = refval - ( refpix + 0.5 * (1 - factor) ) * increment ;
    1920       refpix = (refpix + 0.5)/factor - 0.5;
     1922      if (factor >= 0 ) {
     1923        refpix = (refpix + 0.5)/factor - 0.5;
     1924      } else {
     1925        refpix = (abcissa.size() + 0.5 - refpix)/abs(factor) - 0.5;
     1926      }
    19211927      freqTable_.setEntry( refpix, refval, increment*factor, currId ) ;
    19221928      os << "ID" << currId << ": channel width (Orig) = " << oldincr << " [" << freqTable_.getUnitString() << "], scale factor = " << factor << LogIO::POST ;
    … …  
    19871993  int oldsize = abcissa.size() ;
    19881994  double olddnu = abcissa[1] - abcissa[0] ;
    1989   int ichan = 0 ;
     1995  //int ichan = 0 ;
    19901996  double wsum = 0.0 ;
    19911997  Vector<double> zi( nChan+1 ) ;
    19921998  Vector<double> yi( oldsize + 1 ) ;
    1993   zi[0] = abcissa[0] - 0.5 * olddnu ;
    1994   for ( int ii = 1 ; ii < nChan ; ii++ )
    1995     zi[ii] = zi[0] + dnu * ii ;
    1996   zi[nChan] = zi[nChan-1] + dnu ;
    19971999  yi[0] = abcissa[0] - 0.5 * olddnu ;
    19982000  for ( int ii = 1 ; ii < oldsize ; ii++ )
    … …  
    20002002  yi[oldsize] = abcissa[oldsize-1] \
    20012003    + 0.5 * (abcissa[oldsize-1] - abcissa[oldsize-2]) ;
    2002   if ( dnu > 0.0 ) {
    2003     for ( int ii = 0 ; ii < nChan ; ii++ ) {
    2004       double zl = zi[ii] ;
    2005       double zr = zi[ii+1] ;
    2006       for ( int j = ichan ; j < oldsize ; j++ ) {
    2007         double yl = yi[j] ;
    2008         double yr = yi[j+1] ;
    2009         if ( yl <= zl ) {
    2010           if ( yr <= zl ) {
    2011             continue ;
    2012           }
    2013           else if ( yr <= zr ) {
    2014             if (!oldflag[j]) {
    2015               newspec[ii] += oldspec[j] * ( yr - zl ) ;
    2016               if (regridTsys) newtsys[ii] += oldtsys[j] * ( yr - zl ) ;
    2017               wsum += ( yr - zl ) ;
    2018             }
    2019             newflag[ii] = newflag[ii] && oldflag[j] ;
    2020           }
    2021           else {
    2022             if (!oldflag[j]) {
    2023               newspec[ii] += oldspec[j] * dnu ;
    2024               if (regridTsys) newtsys[ii] += oldtsys[j] * dnu ;
    2025               wsum += dnu ;
    2026             }
    2027             newflag[ii] = newflag[ii] && oldflag[j] ;
    2028             ichan = j ;
    2029             break ;
    2030           }
    2031         }
    2032         else if ( yl < zr ) {
    2033           if ( yr <= zr ) {
    2034             if (!oldflag[j]) {
    2035               newspec[ii] += oldspec[j] * ( yr - yl ) ;
    2036               if (regridTsys) newtsys[ii] += oldtsys[j] * ( yr - yl ) ;
    2037               wsum += ( yr - yl ) ;
    2038             }
    2039             newflag[ii] = newflag[ii] && oldflag[j] ;
    2040           }
    2041           else {
    2042             if (!oldflag[j]) {
    2043               newspec[ii] += oldspec[j] * ( zr - yl ) ;
    2044               if (regridTsys) newtsys[ii] += oldtsys[j] * ( zr - yl ) ;
    2045               wsum += ( zr - yl ) ;
    2046             }
    2047             newflag[ii] = newflag[ii] && oldflag[j] ;
    2048             ichan = j ;
    2049             break ;
    2050           }
    2051         }
    2052         else {
    2053           ichan = j - 1 ;
    2054           break ;
    2055         }
    2056       }
    2057       if ( wsum != 0.0 ) {
    2058         newspec[ii] /= wsum ;
    2059         if (regridTsys) newtsys[ii] /= wsum ;
    2060       }
    2061       wsum = 0.0 ;
    2062     }
    2063   }
    2064   else if ( dnu < 0.0 ) {
    2065     for ( int ii = 0 ; ii < nChan ; ii++ ) {
    2066       double zl = zi[ii] ;
    2067       double zr = zi[ii+1] ;
    2068       for ( int j = ichan ; j < oldsize ; j++ ) {
    2069         double yl = yi[j] ;
    2070         double yr = yi[j+1] ;
    2071         if ( yl >= zl ) {
    2072           if ( yr >= zl ) {
    2073             continue ;
    2074           }
    2075           else if ( yr >= zr ) {
    2076             if (!oldflag[j]) {
    2077               newspec[ii] += oldspec[j] * abs( yr - zl ) ;
    2078               if (regridTsys) newtsys[ii] += oldtsys[j] * abs( yr - zl ) ;
    2079               wsum += abs( yr - zl ) ;
    2080             }
    2081             newflag[ii] = newflag[ii] && oldflag[j] ;
    2082           }
    2083           else {
    2084             if (!oldflag[j]) {
    2085               newspec[ii] += oldspec[j] * abs( dnu ) ;
    2086               if (regridTsys) newtsys[ii] += oldtsys[j] * abs( dnu ) ;
    2087               wsum += abs( dnu ) ;
    2088             }
    2089             newflag[ii] = newflag[ii] && oldflag[j] ;
    2090             ichan = j ;
    2091             break ;
    2092           }
    2093         }
    2094         else if ( yl > zr ) {
    2095           if ( yr >= zr ) {
    2096             if (!oldflag[j]) {
    2097               newspec[ii] += oldspec[j] * abs( yr - yl ) ;
    2098               if (regridTsys) newtsys[ii] += oldtsys[j] * abs( yr - yl ) ;
    2099               wsum += abs( yr - yl ) ;
    2100             }
    2101             newflag[ii] = newflag[ii] && oldflag[j] ;
    2102           }
    2103           else {
    2104             if (!oldflag[j]) {
    2105               newspec[ii] += oldspec[j] * abs( zr - yl ) ;
    2106               if (regridTsys) newtsys[ii] += oldtsys[j] * abs( zr - yl ) ;
    2107               wsum += abs( zr - yl ) ;
    2108             }
    2109             newflag[ii] = newflag[ii] && oldflag[j] ;
    2110             ichan = j ;
    2111             break ;
    2112           }
    2113         }
    2114         else {
    2115           ichan = j - 1 ;
    2116           break ;
    2117         }
    2118       }
    2119       if ( wsum != 0.0 ) {
    2120         newspec[ii] /= wsum ;
    2121         if (regridTsys) newtsys[ii] /= wsum ;
    2122       }
    2123       wsum = 0.0 ;
    2124     }
    2125   }
    2126 //   //ofs << "olddnu = " << olddnu << ", dnu = " << dnu << endl ;
    2127 //   pile += dnu ;
    2128 //   wedge = olddnu * ( refChan + 1 ) ;
    2129 //   while ( wedge < pile ) {
    2130 //     newspec[0] += olddnu * oldspec[refChan] ;
    2131 //     newflag[0] = newflag[0] || oldflag[refChan] ;
    2132 //     //ofs << "channel " << refChan << " is included in new channel 0" << endl ;
    2133 //     refChan++ ;
    2134 //     wedge += olddnu ;
    2135 //     wsum += olddnu ;
    2136 //     //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
     2004  //zi[0] = abcissa[0] - 0.5 * olddnu ;
     2005  zi[0] = ((olddnu*dnu > 0) ? yi[0] : yi[oldsize]) ;
     2006  for ( int ii = 1 ; ii < nChan ; ii++ )
     2007    zi[ii] = zi[0] + dnu * ii ;
     2008  zi[nChan] = zi[nChan-1] + dnu ;
     2009  // Access zi and yi in ascending order
     2010  int izs = ((dnu > 0) ? 0 : nChan ) ;
     2011  int ize = ((dnu > 0) ? nChan : 0 ) ;
     2012  int izincr = ((dnu > 0) ? 1 : -1 ) ;
     2013  int ichan =  ((olddnu > 0) ? 0 : oldsize ) ;
     2014  int iye = ((olddnu > 0) ? oldsize : 0 ) ;
     2015  int iyincr = ((olddnu > 0) ? 1 : -1 ) ;
     2016  //for ( int ii = izs ; ii != ize ; ii+=izincr ){
     2017  int ii = izs ;
     2018  while (ii != ize) {
     2019    // always zl < zr
     2020    double zl = zi[ii] ;
     2021    double zr = zi[ii+izincr] ;
     2022    // Need to access smaller index for the new spec, flag, and tsys.
     2023    // Values between zi[k] and zi[k+1] should be stored in newspec[k], etc.
     2024    int i = min(ii, ii+izincr) ;
     2025    //for ( int jj = ichan ; jj != iye ; jj+=iyincr ) {
     2026    int jj = ichan ;
     2027    while (jj != iye) {
     2028      // always yl < yr
     2029      double yl = yi[jj] ;
     2030      double yr = yi[jj+iyincr] ;
     2031      // Need to access smaller index for the original spec, flag, and tsys.
     2032      // Values between yi[k] and yi[k+1] are stored in oldspec[k], etc.
     2033      int j = min(jj, jj+iyincr) ;
     2034      if ( yr <= zl ) {
     2035        jj += iyincr ;
     2036        continue ;
     2037      }
     2038      else if ( yl <= zl ) {
     2039        if ( yr < zr ) {
     2040          if (!oldflag[j]) {
     2041            newspec[i] += oldspec[j] * ( yr - zl ) ;
     2042            if (regridTsys) newtsys[i] += oldtsys[j] * ( yr - zl ) ;
     2043            wsum += ( yr - zl ) ;
     2044          }
     2045          newflag[i] = newflag[i] && oldflag[j] ;
     2046        }
     2047        else {
     2048          if (!oldflag[j]) {
     2049            newspec[i] += oldspec[j] * abs(dnu) ;
     2050            if (regridTsys) newtsys[i] += oldtsys[j] * abs(dnu) ;
     2051            wsum += abs(dnu) ;
     2052          }
     2053          newflag[i] = newflag[i] && oldflag[j] ;
     2054          ichan = jj ;
     2055          break ;
     2056        }
     2057      }
     2058      else if ( yl < zr ) {
     2059        if ( yr <= zr ) {
     2060          if (!oldflag[j]) {
     2061            newspec[i] += oldspec[j] * ( yr - yl ) ;
     2062            if (regridTsys) newtsys[i] += oldtsys[j] * ( yr - yl ) ;
     2063            wsum += ( yr - yl ) ;
     2064          }
     2065          newflag[i] = newflag[i] && oldflag[j] ;
     2066        }
     2067        else {
     2068          if (!oldflag[j]) {
     2069            newspec[i] += oldspec[j] * ( zr - yl ) ;
     2070            if (regridTsys) newtsys[i] += oldtsys[j] * ( zr - yl ) ;
     2071            wsum += ( zr - yl ) ;
     2072          }
     2073          newflag[i] = newflag[i] && oldflag[j] ;
     2074          ichan = jj ;
     2075          break ;
     2076        }
     2077      }
     2078      else {
     2079        ichan = jj - iyincr ;
     2080        break ;
     2081      }
     2082      jj += iyincr ;
     2083    }
     2084    if ( wsum != 0.0 ) {
     2085      newspec[i] /= wsum ;
     2086      if (regridTsys) newtsys[i] /= wsum ;
     2087    }
     2088    wsum = 0.0 ;
     2089    ii += izincr ;
     2090  }
     2091//   if ( dnu > 0.0 ) {
     2092//     for ( int ii = 0 ; ii < nChan ; ii++ ) {
     2093//       double zl = zi[ii] ;
     2094//       double zr = zi[ii+1] ;
     2095//       for ( int j = ichan ; j < oldsize ; j++ ) {
     2096//         double yl = yi[j] ;
     2097//         double yr = yi[j+1] ;
     2098//         if ( yl <= zl ) {
     2099//           if ( yr <= zl ) {
     2100//             continue ;
     2101//           }
     2102//           else if ( yr <= zr ) {
     2103//          if (!oldflag[j]) {
     2104//            newspec[ii] += oldspec[j] * ( yr - zl ) ;
     2105//            if (regridTsys) newtsys[ii] += oldtsys[j] * ( yr - zl ) ;
     2106//            wsum += ( yr - zl ) ;
     2107//          }
     2108//          newflag[ii] = newflag[ii] && oldflag[j] ;
     2109//           }
     2110//           else {
     2111//          if (!oldflag[j]) {
     2112//            newspec[ii] += oldspec[j] * dnu ;
     2113//            if (regridTsys) newtsys[ii] += oldtsys[j] * dnu ;
     2114//            wsum += dnu ;
     2115//          }
     2116//          newflag[ii] = newflag[ii] && oldflag[j] ;
     2117//             ichan = j ;
     2118//             break ;
     2119//           }
     2120//         }
     2121//         else if ( yl < zr ) {
     2122//           if ( yr <= zr ) {
     2123//          if (!oldflag[j]) {
     2124//            newspec[ii] += oldspec[j] * ( yr - yl ) ;
     2125//            if (regridTsys) newtsys[ii] += oldtsys[j] * ( yr - yl ) ;
     2126//               wsum += ( yr - yl ) ;
     2127//          }
     2128//          newflag[ii] = newflag[ii] && oldflag[j] ;
     2129//           }
     2130//           else {
     2131//          if (!oldflag[j]) {
     2132//            newspec[ii] += oldspec[j] * ( zr - yl ) ;
     2133//            if (regridTsys) newtsys[ii] += oldtsys[j] * ( zr - yl ) ;
     2134//            wsum += ( zr - yl ) ;
     2135//          }
     2136//          newflag[ii] = newflag[ii] && oldflag[j] ;
     2137//             ichan = j ;
     2138//             break ;
     2139//           }
     2140//         }
     2141//         else {
     2142//           ichan = j - 1 ;
     2143//           break ;
     2144//         }
     2145//       }
     2146//       if ( wsum != 0.0 ) {
     2147//         newspec[ii] /= wsum ;
     2148//      if (regridTsys) newtsys[ii] /= wsum ;
     2149//       }
     2150//       wsum = 0.0 ;
     2151//     }
    21372152//   }
    2138 //   frac = ( wedge - pile ) / olddnu ;
    2139 //   wsum += ( 1.0 - frac ) * olddnu ;
    2140 //   newspec[0] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
    2141 //   newflag[0] = newflag[0] || oldflag[refChan] ;
    2142 //   //ofs << "channel " << refChan << " is partly included in new channel 0" << " with fraction of " << ( 1.0 - frac ) << endl ;
    2143 //   //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
    2144 //   newspec[0] /= wsum ;
    2145 //   //ofs << "newspec[0] = " << newspec[0] << endl ;
    2146 //   //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
    2147 
    2148 //   /***
    2149 //    * ichan = 1 - nChan-2
    2150 //    ***/
    2151 //   for ( int ichan = 1 ; ichan < nChan - 1 ; ichan++ ) {
    2152 //     pile += dnu ;
    2153 //     newspec[ichan] += frac * olddnu * oldspec[refChan] ;
    2154 //     newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
    2155 //     //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << frac << endl ;
    2156 //     refChan++ ;
    2157 //     wedge += olddnu ;
    2158 //     wsum = frac * olddnu ;
    2159 //     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
    2160 //     while ( wedge < pile ) {
    2161 //       newspec[ichan] += olddnu * oldspec[refChan] ;
    2162 //       newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
    2163 //       //ofs << "channel " << refChan << " is included in new channel " << ichan << endl ;
    2164 //       refChan++ ;
    2165 //       wedge += olddnu ;
    2166 //       wsum += olddnu ;
    2167 //       //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
     2153//   else if ( dnu < 0.0 ) {
     2154//     for ( int ii = 0 ; ii < nChan ; ii++ ) {
     2155//       double zl = zi[ii] ;
     2156//       double zr = zi[ii+1] ;
     2157//       for ( int j = ichan ; j < oldsize ; j++ ) {
     2158//         double yl = yi[j] ;
     2159//         double yr = yi[j+1] ;
     2160//         if ( yl >= zl ) {
     2161//           if ( yr >= zl ) {
     2162//             continue ;
     2163//           }
     2164//           else if ( yr >= zr ) {
     2165//          if (!oldflag[j]) {
     2166//            newspec[ii] += oldspec[j] * abs( yr - zl ) ;
     2167//            if (regridTsys) newtsys[ii] += oldtsys[j] * abs( yr - zl ) ;
     2168//            wsum += abs( yr - zl ) ;
     2169//          }
     2170//          newflag[ii] = newflag[ii] && oldflag[j] ;
     2171//           }
     2172//           else {
     2173//          if (!oldflag[j]) {
     2174//            newspec[ii] += oldspec[j] * abs( dnu ) ;
     2175//            if (regridTsys) newtsys[ii] += oldtsys[j] * abs( dnu ) ;
     2176//            wsum += abs( dnu ) ;
     2177//          }
     2178//          newflag[ii] = newflag[ii] && oldflag[j] ;
     2179//             ichan = j ;
     2180//             break ;
     2181//           }
     2182//         }
     2183//         else if ( yl > zr ) {
     2184//           if ( yr >= zr ) {
     2185//          if (!oldflag[j]) {
     2186//            newspec[ii] += oldspec[j] * abs( yr - yl ) ;
     2187//            if (regridTsys) newtsys[ii] += oldtsys[j] * abs( yr - yl ) ;
     2188//            wsum += abs( yr - yl ) ;
     2189//          }
     2190//          newflag[ii] = newflag[ii] && oldflag[j] ;
     2191//           }
     2192//           else {
     2193//          if (!oldflag[j]) {
     2194//            newspec[ii] += oldspec[j] * abs( zr - yl ) ;
     2195//            if (regridTsys) newtsys[ii] += oldtsys[j] * abs( zr - yl ) ;
     2196//            wsum += abs( zr - yl ) ;
     2197//          }
     2198//          newflag[ii] = newflag[ii] && oldflag[j] ;
     2199//             ichan = j ;
     2200//             break ;
     2201//           }
     2202//         }
     2203//         else {
     2204//           ichan = j - 1 ;
     2205//           break ;
     2206//         }
     2207//       }
     2208//       if ( wsum != 0.0 ) {
     2209//         newspec[ii] /= wsum ;
     2210//      if (regridTsys) newtsys[ii] /= wsum ;
     2211//       }
     2212//       wsum = 0.0 ;
    21682213//     }
    2169 //     frac = ( wedge - pile ) / olddnu ;
    2170 //     wsum += ( 1.0 - frac ) * olddnu ;
    2171 //     newspec[ichan] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
    2172 //     newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
    2173 //     //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << ( 1.0 - frac ) << endl ;
    2174 //     //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
    2175 //     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
    2176 //     newspec[ichan] /= wsum ;
    2177 //     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << endl ;
    21782214//   }
    2179 
    2180 //   /***
    2181 //    * ichan = nChan-1
    2182 //    ***/
    2183 //   // NOTE: Assumed that all spectra have the same bandwidth
    2184 //   pile += dnu ;
    2185 //   newspec[nChan-1] += frac * olddnu * oldspec[refChan] ;
    2186 //   newflag[nChan-1] = newflag[nChan-1] || oldflag[refChan] ;
    2187 //   //ofs << "channel " << refChan << " is partly included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
    2188 //   refChan++ ;
    2189 //   wedge += olddnu ;
    2190 //   wsum = frac * olddnu ;
    2191 //   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
    2192 //   for ( int jchan = refChan ; jchan < oldsize ; jchan++ ) {
    2193 //     newspec[nChan-1] += olddnu * oldspec[jchan] ;
    2194 //     newflag[nChan-1] = newflag[nChan-1] || oldflag[jchan] ;
    2195 //     wsum += olddnu ;
    2196 //     //ofs << "channel " << jchan << " is included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
    2197 //     //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
    2198 //   }
    2199 //   //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
    2200 //   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
    2201 //   newspec[nChan-1] /= wsum ;
    2202 //   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << endl ;
    2203 
    2204 //   // ofs.close() ;
     2215// //   //ofs << "olddnu = " << olddnu << ", dnu = " << dnu << endl ;
     2216// //   pile += dnu ;
     2217// //   wedge = olddnu * ( refChan + 1 ) ;
     2218// //   while ( wedge < pile ) {
     2219// //     newspec[0] += olddnu * oldspec[refChan] ;
     2220// //     newflag[0] = newflag[0] || oldflag[refChan] ;
     2221// //     //ofs << "channel " << refChan << " is included in new channel 0" << endl ;
     2222// //     refChan++ ;
     2223// //     wedge += olddnu ;
     2224// //     wsum += olddnu ;
     2225// //     //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
     2226// //   }
     2227// //   frac = ( wedge - pile ) / olddnu ;
     2228// //   wsum += ( 1.0 - frac ) * olddnu ;
     2229// //   newspec[0] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
     2230// //   newflag[0] = newflag[0] || oldflag[refChan] ;
     2231// //   //ofs << "channel " << refChan << " is partly included in new channel 0" << " with fraction of " << ( 1.0 - frac ) << endl ;
     2232// //   //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
     2233// //   newspec[0] /= wsum ;
     2234// //   //ofs << "newspec[0] = " << newspec[0] << endl ;
     2235// //   //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
     2236
     2237// //   /***
     2238// //    * ichan = 1 - nChan-2
     2239// //    ***/
     2240// //   for ( int ichan = 1 ; ichan < nChan - 1 ; ichan++ ) {
     2241// //     pile += dnu ;
     2242// //     newspec[ichan] += frac * olddnu * oldspec[refChan] ;
     2243// //     newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
     2244// //     //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << frac << endl ;
     2245// //     refChan++ ;
     2246// //     wedge += olddnu ;
     2247// //     wsum = frac * olddnu ;
     2248// //     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
     2249// //     while ( wedge < pile ) {
     2250// //       newspec[ichan] += olddnu * oldspec[refChan] ;
     2251// //       newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
     2252// //       //ofs << "channel " << refChan << " is included in new channel " << ichan << endl ;
     2253// //       refChan++ ;
     2254// //       wedge += olddnu ;
     2255// //       wsum += olddnu ;
     2256// //       //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
     2257// //     }
     2258// //     frac = ( wedge - pile ) / olddnu ;
     2259// //     wsum += ( 1.0 - frac ) * olddnu ;
     2260// //     newspec[ichan] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
     2261// //     newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
     2262// //     //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << ( 1.0 - frac ) << endl ;
     2263// //     //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
     2264// //     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
     2265// //     newspec[ichan] /= wsum ;
     2266// //     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << endl ;
     2267// //   }
     2268
     2269// //   /***
     2270// //    * ichan = nChan-1
     2271// //    ***/
     2272// //   // NOTE: Assumed that all spectra have the same bandwidth
     2273// //   pile += dnu ;
     2274// //   newspec[nChan-1] += frac * olddnu * oldspec[refChan] ;
     2275// //   newflag[nChan-1] = newflag[nChan-1] || oldflag[refChan] ;
     2276// //   //ofs << "channel " << refChan << " is partly included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
     2277// //   refChan++ ;
     2278// //   wedge += olddnu ;
     2279// //   wsum = frac * olddnu ;
     2280// //   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
     2281// //   for ( int jchan = refChan ; jchan < oldsize ; jchan++ ) {
     2282// //     newspec[nChan-1] += olddnu * oldspec[jchan] ;
     2283// //     newflag[nChan-1] = newflag[nChan-1] || oldflag[jchan] ;
     2284// //     wsum += olddnu ;
     2285// //     //ofs << "channel " << jchan << " is included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
     2286// //     //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
     2287// //   }
     2288// //   //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
     2289// //   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
     2290// //   newspec[nChan-1] /= wsum ;
     2291// //   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << endl ;
     2292
     2293// //   // ofs.close() ;
    22052294
    22062295  specCol_.put( irow, newspec ) ;
Note: See TracChangeset for help on using the changeset viewer.