- Timestamp:
- 05/29/12 11:24:15 (13 years ago)
- Location:
- branches/hpc33/src
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/hpc33/src/MSFiller.cpp
r2549 r2550 304 304 // set dummy epoch 305 305 mf.set( currentTime ) ; 306 307 // initialize dirType308 //dirType = MDirection::N_Types ;309 306 310 307 // … … 957 954 pointingDirection.xyPlane( idx-1 ), pointingDirection.xyPlane( idx ) ) ; 958 955 } 959 //mf.resetEpoch( currentTime ) ;960 956 mf.set( currentTime ) ; 961 957 Quantum< Vector<Double> > tmp( d.column( 0 ), Unit( "rad" ) ) ; 962 958 if ( dirType != MDirection::J2000 ) { 963 //MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;964 959 dir = toj2000( tmp ).getAngle( "rad" ).getValue() ; 965 960 } … … 968 963 } 969 964 if ( dirType != MDirection::AZELGEO ) { 970 //MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ;971 //MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;972 965 azel = toazel( tmp ).getAngle( "rad" ).getValue() ; 973 966 } … … 981 974 { 982 975 dir = sourceDir.getAngle( "rad" ).getValue() ; 983 mf.resetEpoch( currentTime ) ; 984 //MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ; 976 mf.set( currentTime ) ; 985 977 azel = toazel( Quantum< Vector<Double> >( dir, Unit("rad") ) ).getAngle( "rad" ).getValue() ; 986 978 if ( dirType != MDirection::J2000 ) { 987 //MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;988 979 dir = toj2000( Quantum< Vector<Double> >( dir, Unit("rad") ) ).getAngle( "rad" ).getValue() ; 989 980 } -
branches/hpc33/src/STIdxIter.cpp
r2548 r2550 16 16 prod_m.resize( nfield_m ) ; 17 17 idx_m.resize( nfield_m ) ; 18 // prod_m[0] = 1 ;19 // idx_m[0] = 0 ;20 // for ( uInt i = 1 ; i < nfield_m ; i++ ) {21 // prod_m[i] = shape[i-1] * prod_m[i-1] ;22 // idx_m[i] = 0 ;23 // }24 // maxiter_m = prod_m[nfield_m-1] * shape[nfield_m-1] ;25 18 prod_m[nfield_m-1] = 1 ; 26 19 idx_m[nfield_m-1] = 0 ; … … 42 35 niter_m++ ; 43 36 uInt x = niter_m ; 44 // for ( Int i = nfield_m-1 ; i >= 0 ; i-- ) {45 37 for ( Int i = 0 ; i < nfield_m ; i++ ) { 46 38 idx_m[i] = x / prod_m[i] ; … … 136 128 // ... 137 129 storage_m.resize( arr_m.nelements()+nrow_m ) ; 138 // uInt *p = storage_m.storage() + arr_m.nelements() ; 130 // initialize constant temporary storage 139 131 uInt *p = storage_m.storage() ; 140 132 for ( uInt i = 0 ; i < nrow_m ; i++ ) { … … 142 134 p++ ; 143 135 } 136 // len_m layout 137 // len[0]: length of temporary storage that incidates whole rows (constant) 138 // len[1]: number of matched row for column 0 selection 139 // len[2]: number of matched row for column 1 selection 140 // ... 144 141 len_m.resize( ncol_m+1 ) ; 145 142 p = len_m.storage() ; … … 149 146 // cout << "len[" << i << "]=" << len_m[i] << endl ; 150 147 } 148 // skip_m layout 149 // skip_m[0]: True if column 0 is filled by unique value 150 // skip_m[1]: True if column 1 is filled by unique value 151 // ... 152 skip_m.resize( ncol_m ) ; 153 for ( uInt i = 0 ; i < ncol_m ; i++ ) { 154 skip_m[i] = (Bool)(idxlist_m[i].size()==1) ; 155 // cout << "skip_m[" << i << "]=" << skip_m[i] << endl ; 156 } 151 157 prev_m = iter_m->current() ; 152 158 p = prev_m.storage() ; … … 156 162 } 157 163 // cout << "prev_m=" << Vector<uInt>(IPosition(1,ncol_m),prev_m.storage()) << endl ; 158 skip_m.resize( ncol_m ) ;159 for ( uInt i = 0 ; i < ncol_m ; i++ ) {160 skip_m[i] = (Bool)(idxlist_m[i].size()==1) ;161 // cout << "skip_m[" << i << "]=" << skip_m[i] << endl ;162 }163 164 } 164 165 … … 175 176 // cout << "v=" << Vector<uInt>(IPosition(1,v.nelements()),v.storage(),SHARE) << endl ; 176 177 // cout << "c=" << c << endl ; 177 // if ( c < 0 ) {178 178 if ( c > ncol_m-1 ) { 179 // pos_m[0] = len_m[0] ;180 // uInt offset = 0 ;181 // while( offset < ncol_m && skip_m[offset] )182 // offset++ ;183 179 pos_m[0] = len_m[ncol_m] ; 184 180 Int offset = ncol_m - 1 ; … … 189 185 return Vector<uInt>( pos_m, storage_m.storage()+offset*nrow_m, policy ) ; 190 186 } 191 // Int offset = c + 1 ;192 // while( offset < ncol_m && skip_m[offset] )193 // offset++ ;194 187 Int offset = c - 1 ; 195 188 while( offset >= 0 && skip_m[offset] ) … … 200 193 // cout << "len_m[c+1]=" << len_m[c+1] << endl ; 201 194 // cout << "base=" << Vector<uInt>(IPosition(1,len_m[c+1]),base,SHARE) << endl ; 202 // for ( Int i = c ; i >= 0 ; i-- ) { 203 // base = updateStorage( i, base, idxlist_m[i][v[i]] ) ; 204 for ( Int i = c+1 ; i <= ncol_m ; i++ ) { 205 base = updateStorage( i, base, idxlist_m[i-1][v[i-1]] ) ; 195 for ( Int i = c ; i < ncol_m ; i++ ) { 196 base = updateStorage( i, base, idxlist_m[i][v[i]] ) ; 206 197 // cout << "len_m[" << i << "]=" << len_m[i] << endl ; 207 198 // cout << "base=" << Vector<uInt>(IPosition(1,len_m[i]),base,SHARE) << endl ; 208 199 } 209 // pos_m[0] = len_m[0] ;210 200 pos_m[0] = len_m[ncol_m] ; 211 201 // cout << "pos_m=" << pos_m << endl ; … … 216 206 Int ArrayIndexIteratorAcc::isChanged( Block<uInt> &idx ) 217 207 { 218 // Int i = ncol_m - 1 ;219 // while( i >= 0 && idx[i] == prev_m[i] )220 // i-- ;221 208 Int i = 0 ; 222 209 while( i < ncol_m && idx[i] == prev_m[i] ) … … 229 216 uInt &v ) 230 217 { 218 // CAUTION: 219 // indexes for storage_m and len_m differ from index for skip_m by 1 220 // (skip_m[0] corresponds to storage_m[1] and len_m[1]) since there 221 // is additional temporary storage at first segment in storage_m 231 222 uInt *p ; 232 // if ( skip_m[icol] ) { 233 if ( skip_m[icol-1] ) { 223 if ( skip_m[icol] ) { 234 224 // skip update, just update len_m[icol] and pass appropriate pointer 235 225 // cout << "skip " << icol << endl ; 236 226 p = base ; 237 // len_m[icol] = len_m[icol+1] ; 238 len_m[icol] = len_m[icol-1] ; 227 len_m[icol+1] = len_m[icol] ; 239 228 } 240 229 else { 241 230 // cout << "update " << icol << endl ; 242 231 uInt len = 0 ; 243 p = storage_m.storage() + icol* nrow_m ;232 p = storage_m.storage() + (icol+1) * nrow_m ; 244 233 uInt *work = p ; 245 234 Bool b ; … … 247 236 long offset = 0 ; 248 237 // warr_p points a first element of (icol)-th column 249 // const uInt *warr_p = arr_p + icol * nrow_m ; 250 const uInt *warr_p = arr_p + (icol-1) * nrow_m ; 251 // for ( uInt i = 0 ; i < len_m[icol+1] ; i++ ) { 252 for ( uInt i = 0 ; i < len_m[icol-1] ; i++ ) { 238 const uInt *warr_p = arr_p + icol * nrow_m ; 239 for ( uInt i = 0 ; i < len_m[icol] ; i++ ) { 253 240 // increment warr_p by (*(base)-*(base-1)) 254 241 warr_p += *base - offset ; … … 267 254 } 268 255 arr_m.freeStorage( arr_p, b ) ; 269 len_m[icol ] = len ;256 len_m[icol+1] = len ; 270 257 } 271 258 return p ; … … 464 451 { 465 452 convertArray( out, in ) ; 466 //uInt len = in.nelements() ;467 //Vector<Int> tmp = in.copy() ;468 //uInt n = genSort( tmp, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates ) ;469 //val.resize( n ) ;470 //for ( uInt i = 0 ; i < n ; i++ ) {471 // val[i] = tmp[i] ;472 // cout << "val[" << i << "]=" << val[i] << endl ;473 //}474 //out.resize( len ) ;475 //uInt idx = -1 ;476 //cout << "idx=" << idx << endl ;477 //int idx2 = (Int)idx ;478 //cout << "idx2=" << idx2 << endl ;479 //idx = (uInt) idx2 ;480 //cout << "idx=" << idx << endl ;481 //if ( n == 1 ) {482 // cout << "n=1" << endl ;483 // out = 0 ;484 //}485 //else if ( n == 2 ) {486 // cout << "n=2" << endl ;487 // for ( uInt i = 0 ; i < len ; i++ ) {488 // out[i] = (in[i] == val[0]) ? 0 : 1 ;489 // }490 //}491 //else {492 // cout << "n=" << n << endl ;493 // map<Int,uInt> m ;494 // for( uInt i = 0 ; i < n ; i++ )495 // m[val[i]] = i ;496 // for ( uInt i = 0 ; i < len ; i++ ) {497 //for ( uInt j = 0 ; j < n ; j++ ) {498 // if ( in[i] == val[j] ) {499 // out[i] = j ;500 // break ;501 // }502 //}503 // out[i] = m[in[i]] ;504 // }505 // }506 453 } 507 454 … … 518 465 // cout << "val[" << i << "]=" << val[i] << endl ; 519 466 } 520 //out.resize( len ) ;521 467 if ( n == 1 ) { 522 468 //cout << "n=1" << endl ; … … 535 481 m[val[i]] = i ; 536 482 for ( uInt i = 0 ; i < len ; i++ ) { 537 //for ( uInt j = 0 ; j < n ; j++ ) {538 // if ( in[i] == val[j] ) {539 // out[i] = j ;540 // break ;541 // }542 //}543 483 out[i] = m[in[i]] ; 544 484 } … … 548 488 Int STIdxIterExAcc::getSrcType() 549 489 { 550 //if ( srctype_m.nelements() > 0 )551 490 if ( srctypeid_m >= 0 ) 552 //return srctype_m[iter_m->current()[srctypeid_m]] ;553 491 return (Int)(iter_m->current()[srctypeid_m]) ; 554 492 else -
branches/hpc33/src/STMath.cpp
r2546 r2550 184 184 tout.addRow(); 185 185 // t4 = mathutil::gettimeofday_sec() ; 186 //TableCopy::copyRows(tout, subt, outrowCount, 0, 1);187 //TableCopy::copyRows(tout, baset, outrowCount, rows[0], 1);188 186 // skip to copy SPECTRA, FLAGTRA, and TSYS since those heavy columns are 189 187 // overwritten in the following process … … 235 233 // outrowCount += rowNum ; 236 234 237 // merge loop235 // merge loop 238 236 uInt i = outrowCount ; 239 // for (uInt i=0; i < tout.nrow(); ++i) { 240 241 // in[0] is already selected by TableItertor 237 // in[0] is already selected by iterator 242 238 specCol.attach(baset,"SPECTRA"); 243 239 flagCol.attach(baset,"FLAGTRA"); … … 248 244 Vector<uChar> flag; 249 245 Double inter,time; 250 // for (uInt k = 0; k < subt.nrow(); ++k ) { 246 251 247 for (uInt l = 0; l < rows.nelements(); ++l ) { 252 248 uInt k = rows[l] ; … … 273 269 // in[0] is already selected by TableIterator so that index is 274 270 // started from 1 275 //for ( int j=0; j < int(in.size()); ++j ) {276 271 for ( int j=1; j < int(in.size()); ++j ) { 277 272 const Table& tin = in[j]->table(); … … 362 357 const Vector<Bool>& msk = acc.getMask(); 363 358 if ( allEQ(msk, False) ) { 364 //uint n = rowstodelete.nelements();365 //rowstodelete.resize(n+1, True);366 //rowstodelete[n] = i;367 359 rowstodelB[nrowdel] = i ; 368 360 nrowdel++ ; … … 408 400 } 409 401 410 //if (rowstodelete.nelements() > 0) {411 402 if ( nrowdel > 0 ) { 412 403 Vector<uInt> rowstodelete( IPosition(1,nrowdel), rowstodelB.storage(), SHARE ) ; … … 5315 5306 } 5316 5307 5317 void STMath::calibrateALMA( CountedPtr<Scantable>& on,5318 CountedPtr<Scantable>& off )5319 {5320 if ( on->nrow() == 0 )5321 return ;5322 // string reftime = on->getTime( index ) ;5323 ROScalarColumn<Double> timeCol( off->table(), "TIME" ) ;5324 Vector<Double> timeVec = timeCol.getColumn() ;5325 //ROTableColumn timeCol( on->table(), "TIME" ) ;5326 timeCol.attach( on->table(), "TIME" ) ;5327 ArrayColumn<Float> tsysCol( on->table(), "TSYS" ) ;5328 vector<float> sp( on->nchan( on->getIF(0) ) ) ;5329 unsigned int spsize = sp.size() ;5330 for ( int index = 0 ; index < on->nrow() ; index++ ) {5331 double reftime = timeCol.asdouble(index) ;5332 Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ;5333 vector<float> spec = on->getSpectrum( index ) ;5334 Vector<Float> tsys = tsysCol( index ) ;5335 // ALMA Calibration5336 //5337 // Ta* = Tsys * ( ON - OFF ) / OFF5338 //5339 // 2010/01/07 Takeshi Nakazato5340 unsigned int tsyssize = tsys.nelements() ;5341 for ( unsigned int j = 0 ; j < sp.size() ; j++ ) {5342 float tscale = 0.0 ;5343 if ( tsyssize == spsize )5344 tscale = tsys[j] ;5345 else5346 tscale = tsys[0] ;5347 float v = tscale * ( spec[j] - spoff[j] ) / spoff[j] ;5348 sp[j] = v ;5349 }5350 on->setSpectrum( sp, index ) ;5351 }5352 }5353 5354 void STMath::calibrateALMA( CountedPtr<Scantable>& on,5355 CountedPtr<Scantable>& off,5356 Vector<uInt>& rows )5357 {5358 // if rows is empty, just return5359 if ( rows.nelements() == 0 )5360 return ;5361 // string reftime = on->getTime( index ) ;5362 ROScalarColumn<Double> timeCol( off->table(), "TIME" ) ;5363 Vector<Double> timeVec = timeCol.getColumn() ;5364 timeCol.attach( on->table(), "TIME" ) ;5365 ArrayColumn<Float> tsysCol( on->table(), "TSYS" ) ;5366 vector<float> sp( on->nchan( on->getIF(rows[0]) ) ) ;5367 unsigned int spsize = sp.size() ;5368 // I know that the data is contiguous5369 const uInt *p = rows.data() ;5370 for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {5371 double reftime = timeCol.asdouble(*p) ;5372 Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ;5373 vector<float> spec = on->getSpectrum( *p ) ;5374 Vector<Float> tsys = tsysCol( *p ) ;5375 // ALMA Calibration5376 //5377 // Ta* = Tsys * ( ON - OFF ) / OFF5378 //5379 // 2010/01/07 Takeshi Nakazato5380 unsigned int tsyssize = tsys.nelements() ;5381 for ( unsigned int j = 0 ; j < sp.size() ; j++ ) {5382 float tscale = 0.0 ;5383 if ( tsyssize == spsize )5384 tscale = tsys[j] ;5385 else5386 tscale = tsys[0] ;5387 float v = tscale * ( spec[j] - spoff[j] ) / spoff[j] ;5388 sp[j] = v ;5389 }5390 on->setSpectrum( sp, *p ) ;5391 p++ ;5392 }5393 }5394 5395 5308 void STMath::calibrateALMA( CountedPtr<Scantable>& out, 5396 5309 const CountedPtr<Scantable>& on, … … 5409 5322 ROArrayColumn<Float> tsysCol( on->table(), "TSYS" ) ; 5410 5323 unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ; 5411 //Vector<Float> spoff( spsize ) ;5412 //Vector<Float> spec( spsize ) ;5413 5324 // I know that the data is contiguous 5414 5325 const uInt *p = rows.data() ; … … 5417 5328 Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ; 5418 5329 Vector<Float> spec = on->specCol_( *p ) ; 5419 //getSpectrumFromTime2( spoff, reftime, timeVec, off, "linear" ) ;5420 //on->specCol_.get( *p, spec ) ;5421 5330 Vector<Float> tsys = tsysCol( *p ) ; 5422 5331 // ALMA Calibration -
branches/hpc33/src/STMath.h
r2546 r2550 422 422 casa::CountedPtr<Scantable>& off, 423 423 int index ) ; 424 void calibrateALMA( casa::CountedPtr<Scantable>& on,425 casa::CountedPtr<Scantable>& off ) ;426 void calibrateALMA( casa::CountedPtr<Scantable>& on,427 casa::CountedPtr<Scantable>& off,428 casa::Vector<casa::uInt>& rows ) ;429 424 void calibrateALMA( casa::CountedPtr<Scantable>& out, 430 425 const casa::CountedPtr<Scantable>& on,
Note:
See TracChangeset
for help on using the changeset viewer.