- Timestamp:
- 06/15/12 15:35:25 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/hpc33/src/STMath.cpp
r2567 r2571 3955 3955 masks, 3956 3956 "TINT" ) ; 3957 3957 3958 // on scan 3958 3959 // t0 = mathutil::gettimeofday_sec() ; … … 5115 5116 string weight ) 5116 5117 { 5117 Table ttab = s->table() ; 5118 ROScalarColumn<Double> *timeCol = new ROScalarColumn<Double>( ttab, "TIME" ) ; 5119 uInt nrow = timeCol->nrow() ; 5120 Vector<Double> timeSep = timeCol->getColumn() ; 5121 Vector<uInt> scannoOrg = s->scanCol_.getColumn() ; 5122 delete timeCol ; 5123 for ( uInt i = nrow-2 ; i > 0 ; i-- ) { 5124 timeSep[i] -= timeSep[i-1] ; 5125 } 5126 Vector<Double> interval = s->integrCol_.getColumn() ; 5127 interval /= 86400.0 ; 5128 uInt *newscan = new uInt[nrow] ; 5129 Vector<uInt> newscanno( IPosition(1,nrow), newscan, TAKE_OVER ) ; 5130 uInt *p = newscan ; 5131 uInt newid = 0 ; 5132 for ( uInt i = 0 ; i < nrow - 1 ; i++ ) { 5133 *p = newid ; 5134 double gap = 2.0 * timeSep[i+1] / ( interval[i] + interval[i+1] ) ; 5118 // prepare output table 5119 bool insitu = insitu_ ; 5120 insitu_ = false ; 5121 CountedPtr<Scantable> a = getScantable( s, true ) ; 5122 insitu_ = insitu ; 5123 Table &atab = a->table() ; 5124 ScalarColumn<Double> timeColOut( atab, "TIME" ) ; 5125 5126 if ( s->nrow() == 0 ) 5127 return a ; 5128 5129 // setup RowAccumulator 5130 WeightType wtype = stringToWeight( weight ) ; 5131 RowAccumulator acc( wtype ) ; 5132 Vector<Bool> cmask( mask ) ; 5133 acc.setUserMask( cmask ) ; 5134 5135 vector<string> cols( 1, "POLNO" ) ; 5136 STIdxIterAcc iter( s, cols ) ; 5137 5138 uInt nchan = s->nchan(s->getIF(0)) ; 5139 Table ttab = s->table() ; 5140 ROScalarColumn<Double> *timeCol = new ROScalarColumn<Double>( ttab, "TIME" ) ; 5141 Vector<Double> timeVec = timeCol->getColumn() ; 5142 delete timeCol ; 5143 Vector<Double> interval = s->integrCol_.getColumn() ; 5144 uInt nrow = timeVec.nelements() ; 5145 uInt outrow = 0 ; 5146 5147 Vector<uChar> flag( nchan ) ; 5148 Vector<Bool> bflag( nchan ) ; 5149 Vector<Float> spec( nchan ) ; 5150 Vector<Float> tsys( nchan ) ; 5151 5152 while( !iter.pastEnd() ) { 5153 5154 Vector<uInt> rows = iter.getRows( SHARE ) ; 5155 uInt len = rows.nelements() ; 5156 5157 Vector<Double> timeSep( len-1 ) ; 5158 for ( uInt i = 0 ; i < len-1 ; i++ ) { 5159 timeSep[i] = timeVec[rows[i+1]] - timeVec[rows[i]] ; 5160 } 5161 5162 uInt irow ; 5163 uInt jrow ; 5164 for ( uInt i = 0 ; i < len-1 ; i++ ) { 5165 irow = rows[i] ; 5166 jrow = rows[i+1] ; 5167 // accumulate data 5168 s->flagsCol_.get( irow, flag ) ; 5169 convertArray( bflag, flag ) ; 5170 s->specCol_.get( irow, spec ) ; 5171 tsys.assign( s->tsysCol_( irow ) ) ; 5172 if ( !allEQ(bflag,True) ) 5173 acc.add( spec, !bflag, tsys, interval[irow], timeVec[irow] ) ; 5174 double gap = 2.0 * 86400.0 * timeSep[i] / ( interval[jrow] + interval[irow] ) ; 5135 5175 //cout << "gap[" << i << "]=" << setw(5) << gap << endl ; 5136 5176 if ( gap > 1.1 ) { 5137 // cout << "detected gap between " << i << " and " << i+1 << endl ; 5138 newid++ ; 5139 } 5140 p++ ; 5141 } 5142 *p = newid ; 5143 s->scanCol_.putColumn( newscanno ) ; 5144 // double t2 = mathutil::gettimeofday_sec() ; 5145 vector< CountedPtr<Scantable> > tmp( 1, s ) ; 5146 CountedPtr<Scantable> a = average( tmp, mask, weight, "SCAN" ) ; 5147 s->scanCol_.putColumn( scannoOrg ) ; 5148 // double t3 = mathutil::gettimeofday_sec() ; 5149 // cout << "a.nrow = " << a->nrow() << endl ; 5150 // t1 = mathutil::gettimeofday_sec() ; 5151 // cout << " elapsed time for average(): " << t3-t2 << " sec" << endl ; 5152 return a ; 5153 } 5177 //cout << "detected gap between " << i << " and " << i+1 << endl ; 5178 // put data to output table 5179 // reset RowAccumulator 5180 if ( acc.state() ) { 5181 atab.addRow() ; 5182 copyRows( atab, ttab, outrow, irow, 1, False, False, False ) ; 5183 acc.replaceNaN() ; 5184 const Vector<Bool> &msk = acc.getMask() ; 5185 convertArray( flag, !msk ) ; 5186 for (uInt k = 0; k < nchan; ++k) { 5187 uChar userFlag = 1 << 7; 5188 if (msk[k]==True) userFlag = 0 << 7; 5189 flag(k) = userFlag; 5190 } 5191 a->flagsCol_.put( outrow, flag ) ; 5192 a->specCol_.put( outrow, acc.getSpectrum() ) ; 5193 a->tsysCol_.put( outrow, acc.getTsys() ) ; 5194 a->integrCol_.put( outrow, acc.getInterval() ) ; 5195 timeColOut.put( outrow, acc.getTime() ) ; 5196 a->cycleCol_.put( outrow, 0 ) ; 5197 } 5198 acc.reset() ; 5199 outrow++ ; 5200 } 5201 } 5202 5203 // accumulate and add last data 5204 irow = rows[len-1] ; 5205 s->flagsCol_.get( irow, flag ) ; 5206 convertArray( bflag, flag ) ; 5207 s->specCol_.get( irow, spec ) ; 5208 tsys.assign( s->tsysCol_( irow ) ) ; 5209 if (!allEQ(bflag,True) ) 5210 acc.add( spec, !bflag, tsys, interval[irow], timeVec[irow] ) ; 5211 if ( acc.state() ) { 5212 atab.addRow() ; 5213 copyRows( atab, ttab, outrow, irow, 1, False, False, False ) ; 5214 acc.replaceNaN() ; 5215 const Vector<Bool> &msk = acc.getMask() ; 5216 convertArray( flag, !msk ) ; 5217 for (uInt k = 0; k < nchan; ++k) { 5218 uChar userFlag = 1 << 7; 5219 if (msk[k]==True) userFlag = 0 << 7; 5220 flag(k) = userFlag; 5221 } 5222 a->flagsCol_.put( outrow, flag ) ; 5223 a->specCol_.put( outrow, acc.getSpectrum() ) ; 5224 a->tsysCol_.put( outrow, acc.getTsys() ) ; 5225 a->integrCol_.put( outrow, acc.getInterval() ) ; 5226 timeColOut.put( outrow, acc.getTime() ) ; 5227 a->cycleCol_.put( outrow, 0 ) ; 5228 } 5229 acc.reset() ; 5230 outrow++ ; 5231 5232 iter.next() ; 5233 } 5234 5235 return a ; 5236 }
Note:
See TracChangeset
for help on using the changeset viewer.