Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/Scantable.cpp

    r2348 r2463  
    7070#include "STUpgrade.h"
    7171#include "Scantable.h"
     72
     73#define debug 1
    7274
    7375using namespace casa;
    … …  
    18851887}
    18861888
     1889void asap::Scantable::regridSpecChannel( double dnu, int nChan )
     1890{
     1891  LogIO os( LogOrigin( "Scantable", "regridChannel()", WHERE ) ) ;
     1892  os << "Regrid abcissa with spectral resoultion " << dnu << " " << freqTable_.getUnitString() << " with channel number " << ((nChan>0)? String(nChan) : "covering band width")<< LogIO::POST ;
     1893  int freqnrow = freqTable_.table().nrow() ;
     1894  Vector<bool> firstTime( freqnrow, true ) ;
     1895  double oldincr, factor;
     1896  uInt currId;
     1897  Double refpix ;
     1898  Double refval ;
     1899  Double increment ;
     1900  for ( int irow = 0 ; irow < nrow() ; irow++ ) {
     1901    currId = mfreqidCol_(irow);
     1902    vector<double> abcissa = getAbcissa( irow ) ;
     1903    if (nChan < 0) {
     1904      int oldsize = abcissa.size() ;
     1905      double bw = (abcissa[oldsize-1]-abcissa[0]) +                     \
     1906        0.5 * (abcissa[1]-abcissa[0] + abcissa[oldsize-1]-abcissa[oldsize-2]) ;
     1907      nChan = int( ceil( abs(bw/dnu) ) ) ;
     1908    }
     1909    // actual regridding
     1910    regridChannel( nChan, dnu, irow ) ;
     1911
     1912    // update FREQUENCIES subtable
     1913    if (firstTime[currId]) {
     1914      oldincr = abcissa[1]-abcissa[0] ;
     1915      factor = dnu/oldincr ;
     1916      firstTime[currId] = false ;
     1917      freqTable_.getEntry( refpix, refval, increment, currId ) ;
     1918
     1919      //refval = refval - ( refpix + 0.5 * (1 - factor) ) * increment ;
     1920      if (factor > 0 ) {
     1921        refpix = (refpix + 0.5)/factor - 0.5;
     1922      } else {
     1923        refpix = (abcissa.size() - 0.5 - refpix)/abs(factor) - 0.5;
     1924      }
     1925      freqTable_.setEntry( refpix, refval, increment*factor, currId ) ;
     1926      //os << "ID" << currId << ": channel width (Orig) = " << oldincr << " [" << freqTable_.getUnitString() << "], scale factor = " << factor << LogIO::POST ;
     1927      //os << "     frequency increment (Orig) = " << increment << "-> (New) " << increment*factor << LogIO::POST ;
     1928    }
     1929  }
     1930}
     1931
    18871932void asap::Scantable::regridChannel( int nChan, double dnu )
    18881933{
    … …  
    19051950  }
    19061951
    1907   // change channel number for specCol_ and flagCol_
    1908   Vector<Float> newspec( nChan, 0 ) ;
    1909   Vector<uChar> newflag( nChan, false ) ;
     1952  // change channel number for specCol_, flagCol_, and tsysCol_ (if necessary)
    19101953  vector<string> coordinfo = getCoordInfo() ;
    19111954  string oldinfo = coordinfo[0] ;
    … …  
    19331976  Vector<Float> oldspec = specCol_( irow ) ;
    19341977  Vector<uChar> oldflag = flagsCol_( irow ) ;
     1978  Vector<Float> oldtsys = tsysCol_( irow ) ;
    19351979  Vector<Float> newspec( nChan, 0 ) ;
    1936   Vector<uChar> newflag( nChan, false ) ;
     1980  Vector<uChar> newflag( nChan, true ) ;
     1981  Vector<Float> newtsys ;
     1982  bool regridTsys = false ;
     1983  if (oldtsys.size() == oldspec.size()) {
     1984    regridTsys = true ;
     1985    newtsys.resize(nChan,false) ;
     1986    newtsys = 0 ;
     1987  }
    19371988
    19381989  // regrid
    … …  
    19401991  int oldsize = abcissa.size() ;
    19411992  double olddnu = abcissa[1] - abcissa[0] ;
    1942   //int refChan = 0 ;
    1943   //double frac = 0.0 ;
    1944   //double wedge = 0.0 ;
    1945   //double pile = 0.0 ;
    1946   int ichan = 0 ;
     1993  //int ichan = 0 ;
    19471994  double wsum = 0.0 ;
    1948   Vector<Float> zi( nChan+1 ) ;
    1949   Vector<Float> yi( oldsize + 1 ) ;
    1950   zi[0] = abcissa[0] - 0.5 * olddnu ;
    1951   zi[1] = zi[1] + dnu ;
    1952   for ( int ii = 2 ; ii < nChan ; ii++ )
     1995  Vector<double> zi( nChan+1 ) ;
     1996  Vector<double> yi( oldsize + 1 ) ;
     1997  yi[0] = abcissa[0] - 0.5 * olddnu ;
     1998  for ( int ii = 1 ; ii < oldsize ; ii++ )
     1999    yi[ii] = 0.5* (abcissa[ii-1] + abcissa[ii]) ;
     2000  yi[oldsize] = abcissa[oldsize-1] \
     2001    + 0.5 * (abcissa[oldsize-1] - abcissa[oldsize-2]) ;
     2002  //zi[0] = abcissa[0] - 0.5 * olddnu ;
     2003  zi[0] = ((olddnu*dnu > 0) ? yi[0] : yi[oldsize]) ;
     2004  for ( int ii = 1 ; ii < nChan ; ii++ )
    19532005    zi[ii] = zi[0] + dnu * ii ;
    19542006  zi[nChan] = zi[nChan-1] + dnu ;
    1955   yi[0] = abcissa[0] - 0.5 * olddnu ;
    1956   yi[1] = abcissa[1] + 0.5 * olddnu ;
    1957   for ( int ii = 2 ; ii < oldsize ; ii++ )
    1958     yi[ii] = abcissa[ii-1] + olddnu ;
    1959   yi[oldsize] = abcissa[oldsize-1] + 0.5 * olddnu ;
    1960   if ( dnu > 0.0 ) {
    1961     for ( int ii = 0 ; ii < nChan ; ii++ ) {
    1962       double zl = zi[ii] ;
    1963       double zr = zi[ii+1] ;
    1964       for ( int j = ichan ; j < oldsize ; j++ ) {
    1965         double yl = yi[j] ;
    1966         double yr = yi[j+1] ;
    1967         if ( yl <= zl ) {
    1968           if ( yr <= zl ) {
    1969             continue ;
    1970           }
    1971           else if ( yr <= zr ) {
    1972             newspec[ii] += oldspec[j] * ( yr - zl ) ;
    1973             newflag[ii] = newflag[ii] || oldflag[j] ;
    1974             wsum += ( yr - zl ) ;
    1975           }
    1976           else {
    1977             newspec[ii] += oldspec[j] * dnu ;
    1978             newflag[ii] = newflag[ii] || oldflag[j] ;
    1979             wsum += dnu ;
    1980             ichan = j ;
    1981             break ;
    1982           }
    1983         }
    1984         else if ( yl < zr ) {
    1985           if ( yr <= zr ) {
    1986               newspec[ii] += oldspec[j] * ( yr - yl ) ;
    1987               newflag[ii] = newflag[ii] || oldflag[j] ;
    1988               wsum += ( yr - yl ) ;
    1989           }
    1990           else {
    1991             newspec[ii] += oldspec[j] * ( zr - yl ) ;
    1992             newflag[ii] = newflag[ii] || oldflag[j] ;
    1993             wsum += ( zr - yl ) ;
    1994             ichan = j ;
    1995             break ;
    1996           }
    1997         }
    1998         else {
    1999           ichan = j - 1 ;
    2000           break ;
    2001         }
    2002       }
    2003       if ( wsum != 0.0 )
    2004         newspec[ii] /= wsum ;
    2005       wsum = 0.0 ;
    2006     }
    2007   }
    2008   else if ( dnu < 0.0 ) {
    2009     for ( int ii = 0 ; ii < nChan ; ii++ ) {
    2010       double zl = zi[ii] ;
    2011       double zr = zi[ii+1] ;
    2012       for ( int j = ichan ; j < oldsize ; j++ ) {
    2013         double yl = yi[j] ;
    2014         double yr = yi[j+1] ;
    2015         if ( yl >= zl ) {
    2016           if ( yr >= zl ) {
    2017             continue ;
    2018           }
    2019           else if ( yr >= zr ) {
    2020             newspec[ii] += oldspec[j] * abs( yr - zl ) ;
    2021             newflag[ii] = newflag[ii] || oldflag[j] ;
    2022             wsum += abs( yr - zl ) ;
    2023           }
    2024           else {
    2025             newspec[ii] += oldspec[j] * abs( dnu ) ;
    2026             newflag[ii] = newflag[ii] || oldflag[j] ;
    2027             wsum += abs( dnu ) ;
    2028             ichan = j ;
    2029             break ;
    2030           }
    2031         }
    2032         else if ( yl > zr ) {
    2033           if ( yr >= zr ) {
    2034             newspec[ii] += oldspec[j] * abs( yr - yl ) ;
    2035             newflag[ii] = newflag[ii] || oldflag[j] ;
    2036             wsum += abs( yr - yl ) ;
    2037           }
    2038           else {
    2039             newspec[ii] += oldspec[j] * abs( zr - yl ) ;
    2040             newflag[ii] = newflag[ii] || oldflag[j] ;
    2041             wsum += abs( zr - yl ) ;
    2042             ichan = j ;
    2043             break ;
    2044           }
    2045         }
    2046         else {
    2047           ichan = j - 1 ;
    2048           break ;
    2049         }
    2050       }
    2051       if ( wsum != 0.0 )
    2052         newspec[ii] /= wsum ;
    2053       wsum = 0.0 ;
    2054     }
    2055   }
    2056 //    * ichan = 0
    2057 //    ***/
    2058 //   //ofs << "olddnu = " << olddnu << ", dnu = " << dnu << endl ;
    2059 //   pile += dnu ;
    2060 //   wedge = olddnu * ( refChan + 1 ) ;
    2061 //   while ( wedge < pile ) {
    2062 //     newspec[0] += olddnu * oldspec[refChan] ;
    2063 //     newflag[0] = newflag[0] || oldflag[refChan] ;
    2064 //     //ofs << "channel " << refChan << " is included in new channel 0" << endl ;
    2065 //     refChan++ ;
    2066 //     wedge += olddnu ;
    2067 //     wsum += olddnu ;
    2068 //     //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
     2007  // Access zi and yi in ascending order
     2008  int izs = ((dnu > 0) ? 0 : nChan ) ;
     2009  int ize = ((dnu > 0) ? nChan : 0 ) ;
     2010  int izincr = ((dnu > 0) ? 1 : -1 ) ;
     2011  int ichan =  ((olddnu > 0) ? 0 : oldsize ) ;
     2012  int iye = ((olddnu > 0) ? oldsize : 0 ) ;
     2013  int iyincr = ((olddnu > 0) ? 1 : -1 ) ;
     2014  //for ( int ii = izs ; ii != ize ; ii+=izincr ){
     2015  int ii = izs ;
     2016  while (ii != ize) {
     2017    // always zl < zr
     2018    double zl = zi[ii] ;
     2019    double zr = zi[ii+izincr] ;
     2020    // Need to access smaller index for the new spec, flag, and tsys.
     2021    // Values between zi[k] and zi[k+1] should be stored in newspec[k], etc.
     2022    int i = min(ii, ii+izincr) ;
     2023    //for ( int jj = ichan ; jj != iye ; jj+=iyincr ) {
     2024    int jj = ichan ;
     2025    while (jj != iye) {
     2026      // always yl < yr
     2027      double yl = yi[jj] ;
     2028      double yr = yi[jj+iyincr] ;
     2029      // Need to access smaller index for the original spec, flag, and tsys.
     2030      // Values between yi[k] and yi[k+1] are stored in oldspec[k], etc.
     2031      int j = min(jj, jj+iyincr) ;
     2032      if ( yr <= zl ) {
     2033        jj += iyincr ;
     2034        continue ;
     2035      }
     2036      else if ( yl <= zl ) {
     2037        if ( yr < zr ) {
     2038          if (!oldflag[j]) {
     2039            newspec[i] += oldspec[j] * ( yr - zl ) ;
     2040            if (regridTsys) newtsys[i] += oldtsys[j] * ( yr - zl ) ;
     2041            wsum += ( yr - zl ) ;
     2042          }
     2043          newflag[i] = newflag[i] && oldflag[j] ;
     2044        }
     2045        else {
     2046          if (!oldflag[j]) {
     2047            newspec[i] += oldspec[j] * abs(dnu) ;
     2048            if (regridTsys) newtsys[i] += oldtsys[j] * abs(dnu) ;
     2049            wsum += abs(dnu) ;
     2050          }
     2051          newflag[i] = newflag[i] && oldflag[j] ;
     2052          ichan = jj ;
     2053          break ;
     2054        }
     2055      }
     2056      else if ( yl < zr ) {
     2057        if ( yr <= zr ) {
     2058          if (!oldflag[j]) {
     2059            newspec[i] += oldspec[j] * ( yr - yl ) ;
     2060            if (regridTsys) newtsys[i] += oldtsys[j] * ( yr - yl ) ;
     2061            wsum += ( yr - yl ) ;
     2062          }
     2063          newflag[i] = newflag[i] && oldflag[j] ;
     2064        }
     2065        else {
     2066          if (!oldflag[j]) {
     2067            newspec[i] += oldspec[j] * ( zr - yl ) ;
     2068            if (regridTsys) newtsys[i] += oldtsys[j] * ( zr - yl ) ;
     2069            wsum += ( zr - yl ) ;
     2070          }
     2071          newflag[i] = newflag[i] && oldflag[j] ;
     2072          ichan = jj ;
     2073          break ;
     2074        }
     2075      }
     2076      else {
     2077        ichan = jj - iyincr ;
     2078        break ;
     2079      }
     2080      jj += iyincr ;
     2081    }
     2082    if ( wsum != 0.0 ) {
     2083      newspec[i] /= wsum ;
     2084      if (regridTsys) newtsys[i] /= wsum ;
     2085    }
     2086    wsum = 0.0 ;
     2087    ii += izincr ;
     2088  }
     2089//   if ( dnu > 0.0 ) {
     2090//     for ( int ii = 0 ; ii < nChan ; ii++ ) {
     2091//       double zl = zi[ii] ;
     2092//       double zr = zi[ii+1] ;
     2093//       for ( int j = ichan ; j < oldsize ; j++ ) {
     2094//         double yl = yi[j] ;
     2095//         double yr = yi[j+1] ;
     2096//         if ( yl <= zl ) {
     2097//           if ( yr <= zl ) {
     2098//             continue ;
     2099//           }
     2100//           else if ( yr <= zr ) {
     2101//          if (!oldflag[j]) {
     2102//            newspec[ii] += oldspec[j] * ( yr - zl ) ;
     2103//            if (regridTsys) newtsys[ii] += oldtsys[j] * ( yr - zl ) ;
     2104//            wsum += ( yr - zl ) ;
     2105//          }
     2106//          newflag[ii] = newflag[ii] && oldflag[j] ;
     2107//           }
     2108//           else {
     2109//          if (!oldflag[j]) {
     2110//            newspec[ii] += oldspec[j] * dnu ;
     2111//            if (regridTsys) newtsys[ii] += oldtsys[j] * dnu ;
     2112//            wsum += dnu ;
     2113//          }
     2114//          newflag[ii] = newflag[ii] && oldflag[j] ;
     2115//             ichan = j ;
     2116//             break ;
     2117//           }
     2118//         }
     2119//         else if ( yl < zr ) {
     2120//           if ( yr <= zr ) {
     2121//          if (!oldflag[j]) {
     2122//            newspec[ii] += oldspec[j] * ( yr - yl ) ;
     2123//            if (regridTsys) newtsys[ii] += oldtsys[j] * ( yr - yl ) ;
     2124//               wsum += ( yr - yl ) ;
     2125//          }
     2126//          newflag[ii] = newflag[ii] && oldflag[j] ;
     2127//           }
     2128//           else {
     2129//          if (!oldflag[j]) {
     2130//            newspec[ii] += oldspec[j] * ( zr - yl ) ;
     2131//            if (regridTsys) newtsys[ii] += oldtsys[j] * ( zr - yl ) ;
     2132//            wsum += ( zr - yl ) ;
     2133//          }
     2134//          newflag[ii] = newflag[ii] && oldflag[j] ;
     2135//             ichan = j ;
     2136//             break ;
     2137//           }
     2138//         }
     2139//         else {
     2140//           ichan = j - 1 ;
     2141//           break ;
     2142//         }
     2143//       }
     2144//       if ( wsum != 0.0 ) {
     2145//         newspec[ii] /= wsum ;
     2146//      if (regridTsys) newtsys[ii] /= wsum ;
     2147//       }
     2148//       wsum = 0.0 ;
     2149//     }
    20692150//   }
    2070 //   frac = ( wedge - pile ) / olddnu ;
    2071 //   wsum += ( 1.0 - frac ) * olddnu ;
    2072 //   newspec[0] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
    2073 //   newflag[0] = newflag[0] || oldflag[refChan] ;
    2074 //   //ofs << "channel " << refChan << " is partly included in new channel 0" << " with fraction of " << ( 1.0 - frac ) << endl ;
    2075 //   //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
    2076 //   newspec[0] /= wsum ;
    2077 //   //ofs << "newspec[0] = " << newspec[0] << endl ;
    2078 //   //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
    2079 
    2080 //   /***
    2081 //    * ichan = 1 - nChan-2
    2082 //    ***/
    2083 //   for ( int ichan = 1 ; ichan < nChan - 1 ; ichan++ ) {
    2084 //     pile += dnu ;
    2085 //     newspec[ichan] += frac * olddnu * oldspec[refChan] ;
    2086 //     newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
    2087 //     //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << frac << endl ;
    2088 //     refChan++ ;
    2089 //     wedge += olddnu ;
    2090 //     wsum = frac * olddnu ;
    2091 //     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
    2092 //     while ( wedge < pile ) {
    2093 //       newspec[ichan] += olddnu * oldspec[refChan] ;
    2094 //       newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
    2095 //       //ofs << "channel " << refChan << " is included in new channel " << ichan << endl ;
    2096 //       refChan++ ;
    2097 //       wedge += olddnu ;
    2098 //       wsum += olddnu ;
    2099 //       //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
     2151//   else if ( dnu < 0.0 ) {
     2152//     for ( int ii = 0 ; ii < nChan ; ii++ ) {
     2153//       double zl = zi[ii] ;
     2154//       double zr = zi[ii+1] ;
     2155//       for ( int j = ichan ; j < oldsize ; j++ ) {
     2156//         double yl = yi[j] ;
     2157//         double yr = yi[j+1] ;
     2158//         if ( yl >= zl ) {
     2159//           if ( yr >= zl ) {
     2160//             continue ;
     2161//           }
     2162//           else if ( yr >= zr ) {
     2163//          if (!oldflag[j]) {
     2164//            newspec[ii] += oldspec[j] * abs( yr - zl ) ;
     2165//            if (regridTsys) newtsys[ii] += oldtsys[j] * abs( yr - zl ) ;
     2166//            wsum += abs( yr - zl ) ;
     2167//          }
     2168//          newflag[ii] = newflag[ii] && oldflag[j] ;
     2169//           }
     2170//           else {
     2171//          if (!oldflag[j]) {
     2172//            newspec[ii] += oldspec[j] * abs( dnu ) ;
     2173//            if (regridTsys) newtsys[ii] += oldtsys[j] * abs( dnu ) ;
     2174//            wsum += abs( dnu ) ;
     2175//          }
     2176//          newflag[ii] = newflag[ii] && oldflag[j] ;
     2177//             ichan = j ;
     2178//             break ;
     2179//           }
     2180//         }
     2181//         else if ( yl > zr ) {
     2182//           if ( yr >= zr ) {
     2183//          if (!oldflag[j]) {
     2184//            newspec[ii] += oldspec[j] * abs( yr - yl ) ;
     2185//            if (regridTsys) newtsys[ii] += oldtsys[j] * abs( yr - yl ) ;
     2186//            wsum += abs( yr - yl ) ;
     2187//          }
     2188//          newflag[ii] = newflag[ii] && oldflag[j] ;
     2189//           }
     2190//           else {
     2191//          if (!oldflag[j]) {
     2192//            newspec[ii] += oldspec[j] * abs( zr - yl ) ;
     2193//            if (regridTsys) newtsys[ii] += oldtsys[j] * abs( zr - yl ) ;
     2194//            wsum += abs( zr - yl ) ;
     2195//          }
     2196//          newflag[ii] = newflag[ii] && oldflag[j] ;
     2197//             ichan = j ;
     2198//             break ;
     2199//           }
     2200//         }
     2201//         else {
     2202//           ichan = j - 1 ;
     2203//           break ;
     2204//         }
     2205//       }
     2206//       if ( wsum != 0.0 ) {
     2207//         newspec[ii] /= wsum ;
     2208//      if (regridTsys) newtsys[ii] /= wsum ;
     2209//       }
     2210//       wsum = 0.0 ;
    21002211//     }
    2101 //     frac = ( wedge - pile ) / olddnu ;
    2102 //     wsum += ( 1.0 - frac ) * olddnu ;
    2103 //     newspec[ichan] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
    2104 //     newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
    2105 //     //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << ( 1.0 - frac ) << endl ;
    2106 //     //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
    2107 //     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
    2108 //     newspec[ichan] /= wsum ;
    2109 //     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << endl ;
    21102212//   }
    2111 
    2112 //   /***
    2113 //    * ichan = nChan-1
    2114 //    ***/
    2115 //   // NOTE: Assumed that all spectra have the same bandwidth
    2116 //   pile += dnu ;
    2117 //   newspec[nChan-1] += frac * olddnu * oldspec[refChan] ;
    2118 //   newflag[nChan-1] = newflag[nChan-1] || oldflag[refChan] ;
    2119 //   //ofs << "channel " << refChan << " is partly included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
    2120 //   refChan++ ;
    2121 //   wedge += olddnu ;
    2122 //   wsum = frac * olddnu ;
    2123 //   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
    2124 //   for ( int jchan = refChan ; jchan < oldsize ; jchan++ ) {
    2125 //     newspec[nChan-1] += olddnu * oldspec[jchan] ;
    2126 //     newflag[nChan-1] = newflag[nChan-1] || oldflag[jchan] ;
    2127 //     wsum += olddnu ;
    2128 //     //ofs << "channel " << jchan << " is included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
    2129 //     //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
    2130 //   }
    2131 //   //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
    2132 //   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
    2133 //   newspec[nChan-1] /= wsum ;
    2134 //   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << endl ;
    2135 
    2136 //   // ofs.close() ;
     2213// //   //ofs << "olddnu = " << olddnu << ", dnu = " << dnu << endl ;
     2214// //   pile += dnu ;
     2215// //   wedge = olddnu * ( refChan + 1 ) ;
     2216// //   while ( wedge < pile ) {
     2217// //     newspec[0] += olddnu * oldspec[refChan] ;
     2218// //     newflag[0] = newflag[0] || oldflag[refChan] ;
     2219// //     //ofs << "channel " << refChan << " is included in new channel 0" << endl ;
     2220// //     refChan++ ;
     2221// //     wedge += olddnu ;
     2222// //     wsum += olddnu ;
     2223// //     //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
     2224// //   }
     2225// //   frac = ( wedge - pile ) / olddnu ;
     2226// //   wsum += ( 1.0 - frac ) * olddnu ;
     2227// //   newspec[0] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
     2228// //   newflag[0] = newflag[0] || oldflag[refChan] ;
     2229// //   //ofs << "channel " << refChan << " is partly included in new channel 0" << " with fraction of " << ( 1.0 - frac ) << endl ;
     2230// //   //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
     2231// //   newspec[0] /= wsum ;
     2232// //   //ofs << "newspec[0] = " << newspec[0] << endl ;
     2233// //   //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
     2234
     2235// //   /***
     2236// //    * ichan = 1 - nChan-2
     2237// //    ***/
     2238// //   for ( int ichan = 1 ; ichan < nChan - 1 ; ichan++ ) {
     2239// //     pile += dnu ;
     2240// //     newspec[ichan] += frac * olddnu * oldspec[refChan] ;
     2241// //     newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
     2242// //     //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << frac << endl ;
     2243// //     refChan++ ;
     2244// //     wedge += olddnu ;
     2245// //     wsum = frac * olddnu ;
     2246// //     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
     2247// //     while ( wedge < pile ) {
     2248// //       newspec[ichan] += olddnu * oldspec[refChan] ;
     2249// //       newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
     2250// //       //ofs << "channel " << refChan << " is included in new channel " << ichan << endl ;
     2251// //       refChan++ ;
     2252// //       wedge += olddnu ;
     2253// //       wsum += olddnu ;
     2254// //       //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
     2255// //     }
     2256// //     frac = ( wedge - pile ) / olddnu ;
     2257// //     wsum += ( 1.0 - frac ) * olddnu ;
     2258// //     newspec[ichan] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
     2259// //     newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
     2260// //     //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << ( 1.0 - frac ) << endl ;
     2261// //     //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
     2262// //     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
     2263// //     newspec[ichan] /= wsum ;
     2264// //     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << endl ;
     2265// //   }
     2266
     2267// //   /***
     2268// //    * ichan = nChan-1
     2269// //    ***/
     2270// //   // NOTE: Assumed that all spectra have the same bandwidth
     2271// //   pile += dnu ;
     2272// //   newspec[nChan-1] += frac * olddnu * oldspec[refChan] ;
     2273// //   newflag[nChan-1] = newflag[nChan-1] || oldflag[refChan] ;
     2274// //   //ofs << "channel " << refChan << " is partly included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
     2275// //   refChan++ ;
     2276// //   wedge += olddnu ;
     2277// //   wsum = frac * olddnu ;
     2278// //   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
     2279// //   for ( int jchan = refChan ; jchan < oldsize ; jchan++ ) {
     2280// //     newspec[nChan-1] += olddnu * oldspec[jchan] ;
     2281// //     newflag[nChan-1] = newflag[nChan-1] || oldflag[jchan] ;
     2282// //     wsum += olddnu ;
     2283// //     //ofs << "channel " << jchan << " is included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
     2284// //     //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
     2285// //   }
     2286// //   //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
     2287// //   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
     2288// //   newspec[nChan-1] /= wsum ;
     2289// //   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << endl ;
     2290
     2291// //   // ofs.close() ;
    21372292
    21382293  specCol_.put( irow, newspec ) ;
    21392294  flagsCol_.put( irow, newflag ) ;
     2295  if (regridTsys) tsysCol_.put( irow, newtsys );
    21402296
    21412297  return ;
    … …  
    26972853  }
    26982854
    2699   addAuxWaveNumbers(addNWaves, rejectNWaves, nWaves);
     2855  addAuxWaveNumbers(whichrow, addNWaves, rejectNWaves, nWaves);
    27002856}
    27012857
    … …  
    27792935}
    27802936
    2781 void Scantable::addAuxWaveNumbers(const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, std::vector<int>& nWaves)
    2782 {
     2937void Scantable::addAuxWaveNumbers(const int whichrow, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, std::vector<int>& nWaves)
     2938{
     2939  std::vector<int> tempAddNWaves, tempRejectNWaves;
    27832940  for (uInt i = 0; i < addNWaves.size(); ++i) {
     2941    tempAddNWaves.push_back(addNWaves[i]);
     2942  }
     2943  if ((tempAddNWaves.size() == 2) && (tempAddNWaves[1] == -999)) {
     2944    setWaveNumberListUptoNyquistFreq(whichrow, tempAddNWaves);
     2945  }
     2946
     2947  for (uInt i = 0; i < rejectNWaves.size(); ++i) {
     2948    tempRejectNWaves.push_back(rejectNWaves[i]);
     2949  }
     2950  if ((tempRejectNWaves.size() == 2) && (tempRejectNWaves[1] == -999)) {
     2951    setWaveNumberListUptoNyquistFreq(whichrow, tempRejectNWaves);
     2952  }
     2953
     2954  for (uInt i = 0; i < tempAddNWaves.size(); ++i) {
    27842955    bool found = false;
    27852956    for (uInt j = 0; j < nWaves.size(); ++j) {
    2786       if (nWaves[j] == addNWaves[i]) {
     2957      if (nWaves[j] == tempAddNWaves[i]) {
    27872958        found = true;
    27882959        break;
    27892960      }
    27902961    }
    2791     if (!found) nWaves.push_back(addNWaves[i]);
    2792   }
    2793 
    2794   for (uInt i = 0; i < rejectNWaves.size(); ++i) {
     2962    if (!found) nWaves.push_back(tempAddNWaves[i]);
     2963  }
     2964
     2965  for (uInt i = 0; i < tempRejectNWaves.size(); ++i) {
    27952966    for (std::vector<int>::iterator j = nWaves.begin(); j != nWaves.end(); ) {
    2796       if (*j == rejectNWaves[i]) {
     2967      if (*j == tempRejectNWaves[i]) {
    27972968        j = nWaves.erase(j);
    27982969      } else {
    … …  
    28052976    sort(nWaves.begin(), nWaves.end());
    28062977    unique(nWaves.begin(), nWaves.end());
     2978  }
     2979}
     2980
     2981void Scantable::setWaveNumberListUptoNyquistFreq(const int whichrow, std::vector<int>& nWaves)
     2982{
     2983  if ((nWaves.size() == 2)&&(nWaves[1] == -999)) {
     2984    int val = nWaves[0];
     2985    int nyquistFreq = nchan(getIF(whichrow))/2+1;
     2986    nWaves.clear();
     2987    if (val > nyquistFreq) {  // for safety, at least nWaves contains a constant; CAS-3759
     2988      nWaves.push_back(0);
     2989    }
     2990    while (val <= nyquistFreq) {
     2991      nWaves.push_back(val);
     2992      val++;
     2993    }
    28072994  }
    28082995}
    … …  
    28443031
    28453032      //FOR DEBUGGING------------
     3033      /*
    28463034      if (whichrow < 0) {// == nRow -1) {
    28473035        cout << "+++ i=" << setw(3) << whichrow << ", IF=" << setw(2) << getIF(whichrow);
    … …  
    28553043        cout << flush;
    28563044      }
     3045      */
    28573046      //-------------------------
    28583047
    … …  
    31893378  std::vector<bool> mask = getMask(whichrow);
    31903379  uInt maskSize = mask.size();
    3191   if (maskSize != inMask.size()) {
    3192     throw(AipsError("mask sizes are not the same."));
    3193   }
    3194   for (uInt i = 0; i < maskSize; ++i) {
    3195     mask[i] = mask[i] && inMask[i];
     3380  if (inMask.size() != 0) {
     3381    if (maskSize != inMask.size()) {
     3382      throw(AipsError("mask sizes are not the same."));
     3383    }
     3384    for (uInt i = 0; i < maskSize; ++i) {
     3385      mask[i] = mask[i] && inMask[i];
     3386    }
    31963387  }
    31973388
    … …  
    35153706std::string Scantable::getMaskRangeList(const std::vector<bool>& mask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, bool verbose)
    35163707{
    3517   if (mask.size() < 2) {
    3518     throw(AipsError("The mask elements should be > 1"));
     3708  if (mask.size() <= 0) {
     3709    throw(AipsError("The mask elements should be > 0"));
    35193710  }
    35203711  int IF = getIF(whichrow);
    … …  
    35493740std::vector<int> Scantable::getMaskEdgeIndices(const std::vector<bool>& mask)
    35503741{
    3551   if (mask.size() < 2) {
    3552     throw(AipsError("The mask elements should be > 1"));
     3742  if (mask.size() <= 0) {
     3743    throw(AipsError("The mask elements should be > 0"));
    35533744  }
    35543745
Note: See TracChangeset for help on using the changeset viewer.