Changeset 1819 for trunk/src/Scantable.cpp
- Timestamp:
- 08/02/10 17:28:20 (14 years ago)
- Location:
- trunk/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src
-
Property
svn:mergeinfo
set to
/branches/alma/src merged eligible
-
Property
svn:mergeinfo
set to
-
trunk/src/Scantable.cpp
r1743 r1819 11 11 // 12 12 #include <map> 13 #include <fstream> 13 14 14 15 #include <casa/aips.h> … … 24 25 #include <casa/Arrays/Vector.h> 25 26 #include <casa/Arrays/VectorSTLIterator.h> 27 #include <casa/Arrays/Slice.h> 26 28 #include <casa/BasicMath/Math.h> 27 29 #include <casa/BasicSL/Constants.h> … … 29 31 #include <casa/Containers/RecordField.h> 30 32 #include <casa/Utilities/GenSort.h> 33 #include <casa/Logging/LogIO.h> 31 34 32 35 #include <tables/Tables/TableParse.h> … … 106 109 { 107 110 initFactories(); 111 108 112 Table tab(name, Table::Update); 109 113 uInt version = tab.keywordSet().asuInt("VERSION"); … … 121 125 attach(); 122 126 } 127 /* 128 Scantable::Scantable(const std::string& name, Table::TableType ttype) : 129 type_(ttype) 130 { 131 initFactories(); 132 Table tab(name, Table::Update); 133 uInt version = tab.keywordSet().asuInt("VERSION"); 134 if (version != version_) { 135 throw(AipsError("Unsupported version of ASAP file.")); 136 } 137 if ( type_ == Table::Memory ) { 138 table_ = tab.copyToMemoryTable(generateName()); 139 } else { 140 table_ = tab; 141 } 142 143 attachSubtables(); 144 originalTable_ = table_; 145 attach(); 146 } 147 */ 123 148 124 149 Scantable::Scantable( const Scantable& other, bool clear ) … … 199 224 td.addColumn(ScalarColumnDesc<uInt>("FREQ_ID")); 200 225 td.addColumn(ScalarColumnDesc<uInt>("MOLECULE_ID")); 201 td.addColumn(ScalarColumnDesc<Int>("REFBEAMNO")); 226 227 ScalarColumnDesc<Int> refbeamnoColumn("REFBEAMNO"); 228 refbeamnoColumn.setDefault(Int(-1)); 229 td.addColumn(refbeamnoColumn); 230 231 ScalarColumnDesc<uInt> flagrowColumn("FLAGROW"); 232 flagrowColumn.setDefault(uInt(0)); 233 td.addColumn(flagrowColumn); 202 234 203 235 td.addColumn(ScalarColumnDesc<Double>("TIME")); … … 257 289 originalTable_ = table_; 258 290 } 259 260 291 261 292 void Scantable::attach() … … 285 316 mfocusidCol_.attach(table_, "FOCUS_ID"); 286 317 mmolidCol_.attach(table_, "MOLECULE_ID"); 318 319 //Add auxiliary column for row-based flagging (CAS-1433 Wataru Kawasaki) 320 attachAuxColumnDef(flagrowCol_, "FLAGROW", 0); 321 322 } 323 324 template<class T, class T2> 325 void Scantable::attachAuxColumnDef(ScalarColumn<T>& col, 326 const String& colName, 327 const T2& defValue) 328 { 329 try { 330 col.attach(table_, colName); 331 } catch (TableError& err) { 332 String errMesg = err.getMesg(); 333 if (errMesg == "Table column " + colName + " is unknown") { 334 table_.addColumn(ScalarColumnDesc<T>(colName)); 335 col.attach(table_, colName); 336 col.fillColumn(static_cast<T>(defValue)); 337 } else { 338 throw; 339 } 340 } catch (...) { 341 throw; 342 } 343 } 344 345 template<class T, class T2> 346 void Scantable::attachAuxColumnDef(ArrayColumn<T>& col, 347 const String& colName, 348 const Array<T2>& defValue) 349 { 350 try { 351 col.attach(table_, colName); 352 } catch (TableError& err) { 353 String errMesg = err.getMesg(); 354 if (errMesg == "Table column " + colName + " is unknown") { 355 table_.addColumn(ArrayColumnDesc<T>(colName)); 356 col.attach(table_, colName); 357 358 int size = 0; 359 ArrayIterator<T2>& it = defValue.begin(); 360 while (it != defValue.end()) { 361 ++size; 362 ++it; 363 } 364 IPosition ip(1, size); 365 Array<T>& arr(ip); 366 for (int i = 0; i < size; ++i) 367 arr[i] = static_cast<T>(defValue[i]); 368 369 col.fillColumn(arr); 370 } else { 371 throw; 372 } 373 } catch (...) { 374 throw; 375 } 287 376 } 288 377 … … 627 716 } 628 717 718 void Scantable::clip(const Float uthres, const Float dthres, bool clipoutside, bool unflag) 719 { 720 for (uInt i=0; i<table_.nrow(); ++i) { 721 Vector<uChar> flgs = flagsCol_(i); 722 srchChannelsToClip(i, uthres, dthres, clipoutside, unflag, flgs); 723 flagsCol_.put(i, flgs); 724 } 725 } 726 727 std::vector<bool> Scantable::getClipMask(int whichrow, const Float uthres, const Float dthres, bool clipoutside, bool unflag) 728 { 729 Vector<uChar> flags; 730 flagsCol_.get(uInt(whichrow), flags); 731 srchChannelsToClip(uInt(whichrow), uthres, dthres, clipoutside, unflag, flags); 732 Vector<Bool> bflag(flags.shape()); 733 convertArray(bflag, flags); 734 //bflag = !bflag; 735 736 std::vector<bool> mask; 737 bflag.tovector(mask); 738 return mask; 739 } 740 741 void Scantable::srchChannelsToClip(uInt whichrow, const Float uthres, const Float dthres, bool clipoutside, bool unflag, 742 Vector<uChar> flgs) 743 { 744 Vector<Float> spcs = specCol_(whichrow); 745 uInt nchannel = nchan(); 746 if (spcs.nelements() != nchannel) { 747 throw(AipsError("Data has incorrect number of channels")); 748 } 749 uChar userflag = 1 << 7; 750 if (unflag) { 751 userflag = 0 << 7; 752 } 753 if (clipoutside) { 754 for (uInt j = 0; j < nchannel; ++j) { 755 Float spc = spcs(j); 756 if ((spc >= uthres) || (spc <= dthres)) { 757 flgs(j) = userflag; 758 } 759 } 760 } else { 761 for (uInt j = 0; j < nchannel; ++j) { 762 Float spc = spcs(j); 763 if ((spc < uthres) && (spc > dthres)) { 764 flgs(j) = userflag; 765 } 766 } 767 } 768 } 769 629 770 void Scantable::flag(const std::vector<bool>& msk, bool unflag) 630 771 { … … 673 814 flagsCol_.put(i, flgs); 674 815 } 816 } 817 818 void Scantable::flagRow(const std::vector<uInt>& rows, bool unflag) 819 { 820 if ( selector_.empty() && (rows.size() == table_.nrow()) ) 821 throw(AipsError("Trying to flag whole scantable.")); 822 823 uInt rowflag = (unflag ? 0 : 1); 824 std::vector<uInt>::const_iterator it; 825 for (it = rows.begin(); it != rows.end(); ++it) 826 flagrowCol_.put(*it, rowflag); 675 827 } 676 828 … … 793 945 table_.keywordSet().get("FluxUnit", tmp); 794 946 oss << setw(15) << "Flux Unit:" << tmp << endl; 795 Vector<Double> vec(moleculeTable_.getRestFrequencies()); 947 //Vector<Double> vec(moleculeTable_.getRestFrequencies()); 948 int nid = moleculeTable_.nrow(); 949 Bool firstline = True; 796 950 oss << setw(15) << "Rest Freqs:"; 797 if (vec.nelements() > 0) { 798 oss << setprecision(10) << vec << " [Hz]" << endl; 799 } else { 800 oss << "none" << endl; 951 for (int i=0; i<nid; i++) { 952 Table t = table_(table_.col("MOLECULE_ID") == i); 953 if (t.nrow() > 0) { 954 Vector<Double> vec(moleculeTable_.getRestFrequency(i)); 955 if (vec.nelements() > 0) { 956 if (firstline) { 957 oss << setprecision(10) << vec << " [Hz]" << endl; 958 firstline=False; 959 } 960 else{ 961 oss << setw(15)<<" " << setprecision(10) << vec << " [Hz]" << endl; 962 } 963 } else { 964 oss << "none" << endl; 965 } 966 } 801 967 } 802 968 … … 901 1067 const MDirection& md = getDirection(whichrow); 902 1068 const MEpoch& me = timeCol_(whichrow); 903 Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow)); 1069 //Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow)); 1070 Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow)); 904 1071 return freqTable_.getSpectralCoordinate(md, mp, me, rf, 905 1072 mfreqidCol_(whichrow)); … … 975 1142 const MDirection& md = getDirection(whichrow); 976 1143 const MEpoch& me = timeCol_(whichrow); 977 const Double& rf = mmolidCol_(whichrow); 1144 //const Double& rf = mmolidCol_(whichrow); 1145 const Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow)); 978 1146 SpectralCoordinate spc = 979 1147 freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow)); … … 992 1160 } 993 1161 994 void Scantable::setRestFrequencies( double rf, const std::string& name, 1162 /** 1163 void asap::Scantable::setRestFrequencies( double rf, const std::string& name, 995 1164 const std::string& unit ) 1165 **/ 1166 void Scantable::setRestFrequencies( vector<double> rf, const vector<std::string>& name, 1167 const std::string& unit ) 1168 996 1169 { 997 1170 ///@todo lookup in line table to fill in name and formattedname 998 1171 Unit u(unit); 999 Quantum<Double> urf(rf, u); 1000 uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), name, ""); 1172 //Quantum<Double> urf(rf, u); 1173 Quantum<Vector<Double> >urf(rf, u); 1174 Vector<String> formattedname(0); 1175 //cerr<<"Scantable::setRestFrequnecies="<<urf<<endl; 1176 1177 //uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), name, ""); 1178 uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), mathutil::toVectorString(name), formattedname); 1001 1179 TableVector<uInt> tabvec(table_, "MOLECULE_ID"); 1002 1180 tabvec = id; 1003 1181 } 1004 1182 1005 void Scantable::setRestFrequencies( const std::string& name ) 1183 /** 1184 void asap::Scantable::setRestFrequencies( const std::string& name ) 1006 1185 { 1007 1186 throw(AipsError("setRestFrequencies( const std::string& name ) NYI")); 1187 ///@todo implement 1188 } 1189 **/ 1190 void Scantable::setRestFrequencies( const vector<std::string>& name ) 1191 { 1192 throw(AipsError("setRestFrequencies( const vector<std::string>& name ) NYI")); 1008 1193 ///@todo implement 1009 1194 } … … 1044 1229 void Scantable::addFit( const STFitEntry& fit, int row ) 1045 1230 { 1046 cout << mfitidCol_(uInt(row)) << endl; 1231 //cout << mfitidCol_(uInt(row)) << endl; 1232 LogIO os( LogOrigin( "Scantable", "addFit()", WHERE ) ) ; 1233 os << mfitidCol_(uInt(row)) << LogIO::POST ; 1047 1234 uInt id = fitTable_.addEntry(fit, mfitidCol_(uInt(row))); 1048 1235 mfitidCol_.put(uInt(row), id); … … 1059 1246 } 1060 1247 1061 std::string Scantable::getAntennaName() const1248 String Scantable::getAntennaName() const 1062 1249 { 1063 1250 String out; … … 1078 1265 Table subt = t( t.col("SCAN") == scanlist[i]+1 ); 1079 1266 if (subt.nrow()==0) { 1080 cerr <<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<endl; 1267 //cerr <<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<endl; 1268 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ; 1269 os <<LogIO::WARN<<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<LogIO::POST; 1081 1270 ret = 1; 1082 1271 break; … … 1090 1279 Table subt2 = t( t.col("SCAN") == scanlist[i+1]+1 ); 1091 1280 if ( subt2.nrow() == 0) { 1092 cerr<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<endl; 1281 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ; 1282 1283 //cerr<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<endl; 1284 os<<LogIO::WARN<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<LogIO::POST; 1093 1285 ret = 1; 1094 1286 break; … … 1100 1292 if (scan1seqn == 1 && scan2seqn == 2) { 1101 1293 if (laston1 == laston2) { 1102 cerr<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl; 1294 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ; 1295 //cerr<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl; 1296 os<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST; 1103 1297 i +=1; 1104 1298 } 1105 1299 else { 1106 cerr<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl; 1300 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ; 1301 //cerr<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl; 1302 os<<LogIO::WARN<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST; 1107 1303 } 1108 1304 } 1109 1305 else if (scan1seqn==2 && scan2seqn == 1) { 1110 1306 if (laston1 == laston2) { 1111 cerr<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<endl; 1307 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ; 1308 //cerr<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<endl; 1309 os<<LogIO::WARN<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<LogIO::POST; 1112 1310 ret = 1; 1113 1311 break; … … 1115 1313 } 1116 1314 else { 1117 cerr<<"The other scan for "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<endl; 1315 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ; 1316 //cerr<<"The other scan for "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<endl; 1317 os<<LogIO::WARN<<"The other scan for "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<LogIO::POST; 1118 1318 ret = 1; 1119 1319 break; … … 1122 1322 } 1123 1323 else { 1124 cerr<<"The scan does not appear to be standard obsevation."<<endl; 1324 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ; 1325 //cerr<<"The scan does not appear to be standard obsevation."<<endl; 1326 os<<LogIO::WARN<<"The scan does not appear to be standard obsevation."<<LogIO::POST; 1125 1327 } 1126 1328 //if ( i >= nscan ) break; … … 1128 1330 } 1129 1331 else { 1130 cerr<<"No reference to GBT_GO table."<<endl; 1332 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ; 1333 //cerr<<"No reference to GBT_GO table."<<endl; 1334 os<<LogIO::WARN<<"No reference to GBT_GO table."<<LogIO::POST; 1131 1335 ret = 1; 1132 1336 } … … 1142 1346 } 1143 1347 1348 void asap::Scantable::reshapeSpectrum( int nmin, int nmax ) 1349 throw( casa::AipsError ) 1350 { 1351 // assumed that all rows have same nChan 1352 Vector<Float> arr = specCol_( 0 ) ; 1353 int nChan = arr.nelements() ; 1354 1355 // if nmin < 0 or nmax < 0, nothing to do 1356 if ( nmin < 0 ) { 1357 throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ; 1358 } 1359 if ( nmax < 0 ) { 1360 throw( casa::indexError<int>( nmax, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ; 1361 } 1362 1363 // if nmin > nmax, exchange values 1364 if ( nmin > nmax ) { 1365 int tmp = nmax ; 1366 nmax = nmin ; 1367 nmin = tmp ; 1368 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ; 1369 os << "Swap values. Applied range is [" 1370 << nmin << ", " << nmax << "]" << LogIO::POST ; 1371 } 1372 1373 // if nmin exceeds nChan, nothing to do 1374 if ( nmin >= nChan ) { 1375 throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Specified minimum exceeds nChan." ) ) ; 1376 } 1377 1378 // if nmax exceeds nChan, reset nmax to nChan 1379 if ( nmax >= nChan ) { 1380 if ( nmin == 0 ) { 1381 // nothing to do 1382 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ; 1383 os << "Whole range is selected. Nothing to do." << LogIO::POST ; 1384 return ; 1385 } 1386 else { 1387 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ; 1388 os << "Specified maximum exceeds nChan. Applied range is [" 1389 << nmin << ", " << nChan-1 << "]." << LogIO::POST ; 1390 nmax = nChan - 1 ; 1391 } 1392 } 1393 1394 // reshape specCol_ and flagCol_ 1395 for ( int irow = 0 ; irow < nrow() ; irow++ ) { 1396 reshapeSpectrum( nmin, nmax, irow ) ; 1397 } 1398 1399 // update FREQUENCIES subtable 1400 Double refpix ; 1401 Double refval ; 1402 Double increment ; 1403 int freqnrow = freqTable_.table().nrow() ; 1404 Vector<uInt> oldId( freqnrow ) ; 1405 Vector<uInt> newId( freqnrow ) ; 1406 for ( int irow = 0 ; irow < freqnrow ; irow++ ) { 1407 freqTable_.getEntry( refpix, refval, increment, irow ) ; 1408 /*** 1409 * need to shift refpix to nmin 1410 * note that channel nmin in old index will be channel 0 in new one 1411 ***/ 1412 refval = refval - ( refpix - nmin ) * increment ; 1413 refpix = 0 ; 1414 freqTable_.setEntry( refpix, refval, increment, irow ) ; 1415 } 1416 1417 // update nchan 1418 int newsize = nmax - nmin + 1 ; 1419 table_.rwKeywordSet().define( "nChan", newsize ) ; 1420 1421 // update bandwidth 1422 // assumed all spectra in the scantable have same bandwidth 1423 table_.rwKeywordSet().define( "Bandwidth", increment * newsize ) ; 1424 1425 return ; 1426 } 1427 1428 void asap::Scantable::reshapeSpectrum( int nmin, int nmax, int irow ) 1429 { 1430 // reshape specCol_ and flagCol_ 1431 Vector<Float> oldspec = specCol_( irow ) ; 1432 Vector<uChar> oldflag = flagsCol_( irow ) ; 1433 uInt newsize = nmax - nmin + 1 ; 1434 specCol_.put( irow, oldspec( Slice( nmin, newsize, 1 ) ) ) ; 1435 flagsCol_.put( irow, oldflag( Slice( nmin, newsize, 1 ) ) ) ; 1436 1437 return ; 1438 } 1439 1440 void asap::Scantable::regridChannel( int nChan, double dnu ) 1441 { 1442 LogIO os( LogOrigin( "Scantable", "regridChannel()", WHERE ) ) ; 1443 os << "Regrid abcissa with channel number " << nChan << " and spectral resoultion " << dnu << "Hz." << LogIO::POST ; 1444 // assumed that all rows have same nChan 1445 Vector<Float> arr = specCol_( 0 ) ; 1446 int oldsize = arr.nelements() ; 1447 1448 // if oldsize == nChan, nothing to do 1449 if ( oldsize == nChan ) { 1450 os << "Specified channel number is same as current one. Nothing to do." << LogIO::POST ; 1451 return ; 1452 } 1453 1454 // if oldChan < nChan, unphysical operation 1455 if ( oldsize < nChan ) { 1456 os << "Unphysical operation. Nothing to do." << LogIO::POST ; 1457 return ; 1458 } 1459 1460 // change channel number for specCol_ and flagCol_ 1461 Vector<Float> newspec( nChan, 0 ) ; 1462 Vector<uChar> newflag( nChan, false ) ; 1463 vector<string> coordinfo = getCoordInfo() ; 1464 string oldinfo = coordinfo[0] ; 1465 coordinfo[0] = "Hz" ; 1466 setCoordInfo( coordinfo ) ; 1467 for ( int irow = 0 ; irow < nrow() ; irow++ ) { 1468 regridChannel( nChan, dnu, irow ) ; 1469 } 1470 coordinfo[0] = oldinfo ; 1471 setCoordInfo( coordinfo ) ; 1472 1473 1474 // NOTE: this method does not update metadata such as 1475 // FREQUENCIES subtable, nChan, Bandwidth, etc. 1476 1477 return ; 1478 } 1479 1480 void asap::Scantable::regridChannel( int nChan, double dnu, int irow ) 1481 { 1482 // logging 1483 //ofstream ofs( "average.log", std::ios::out | std::ios::app ) ; 1484 //ofs << "IFNO = " << getIF( irow ) << " irow = " << irow << endl ; 1485 1486 Vector<Float> oldspec = specCol_( irow ) ; 1487 Vector<uChar> oldflag = flagsCol_( irow ) ; 1488 Vector<Float> newspec( nChan, 0 ) ; 1489 Vector<uChar> newflag( nChan, false ) ; 1490 1491 // regrid 1492 vector<double> abcissa = getAbcissa( irow ) ; 1493 int oldsize = abcissa.size() ; 1494 double olddnu = abcissa[1] - abcissa[0] ; 1495 //int refChan = 0 ; 1496 //double frac = 0.0 ; 1497 //double wedge = 0.0 ; 1498 //double pile = 0.0 ; 1499 int ichan = 0 ; 1500 double wsum = 0.0 ; 1501 Vector<Float> z( nChan ) ; 1502 z[0] = abcissa[0] - 0.5 * olddnu + 0.5 * dnu ; 1503 for ( int ii = 1 ; ii < nChan ; ii++ ) 1504 z[ii] = z[ii-1] + dnu ; 1505 Vector<Float> zi( nChan+1 ) ; 1506 Vector<Float> yi( oldsize + 1 ) ; 1507 zi[0] = z[0] - 0.5 * dnu ; 1508 zi[1] = z[0] + 0.5 * dnu ; 1509 for ( int ii = 2 ; ii < nChan ; ii++ ) 1510 zi[ii] = zi[ii-1] + dnu ; 1511 zi[nChan] = z[nChan-1] + 0.5 * dnu ; 1512 yi[0] = abcissa[0] - 0.5 * olddnu ; 1513 yi[1] = abcissa[1] + 0.5 * olddnu ; 1514 for ( int ii = 2 ; ii < oldsize ; ii++ ) 1515 yi[ii] = abcissa[ii-1] + olddnu ; 1516 yi[oldsize] = abcissa[oldsize-1] + 0.5 * olddnu ; 1517 if ( dnu > 0.0 ) { 1518 for ( int ii = 0 ; ii < nChan ; ii++ ) { 1519 double zl = zi[ii] ; 1520 double zr = zi[ii+1] ; 1521 for ( int j = ichan ; j < oldsize ; j++ ) { 1522 double yl = yi[j] ; 1523 double yr = yi[j+1] ; 1524 if ( yl <= zl ) { 1525 if ( yr <= zl ) { 1526 continue ; 1527 } 1528 else if ( yr <= zr ) { 1529 newspec[ii] += oldspec[j] * ( yr - zl ) ; 1530 newflag[ii] = newflag[ii] || oldflag[j] ; 1531 wsum += ( yr - zl ) ; 1532 } 1533 else { 1534 newspec[ii] += oldspec[j] * dnu ; 1535 newflag[ii] = newflag[ii] || oldflag[j] ; 1536 wsum += dnu ; 1537 ichan = j ; 1538 break ; 1539 } 1540 } 1541 else if ( yl < zr ) { 1542 if ( yr <= zr ) { 1543 newspec[ii] += oldspec[j] * ( yr - yl ) ; 1544 newflag[ii] = newflag[ii] || oldflag[j] ; 1545 wsum += ( yr - yl ) ; 1546 } 1547 else { 1548 newspec[ii] += oldspec[j] * ( zr - yl ) ; 1549 newflag[ii] = newflag[ii] || oldflag[j] ; 1550 wsum += ( zr - yl ) ; 1551 ichan = j ; 1552 break ; 1553 } 1554 } 1555 else { 1556 ichan = j - 1 ; 1557 break ; 1558 } 1559 } 1560 newspec[ii] /= wsum ; 1561 wsum = 0.0 ; 1562 } 1563 } 1564 else if ( dnu < 0.0 ) { 1565 for ( int ii = 0 ; ii < nChan ; ii++ ) { 1566 double zl = zi[ii] ; 1567 double zr = zi[ii+1] ; 1568 for ( int j = ichan ; j < oldsize ; j++ ) { 1569 double yl = yi[j] ; 1570 double yr = yi[j+1] ; 1571 if ( yl >= zl ) { 1572 if ( yr >= zl ) { 1573 continue ; 1574 } 1575 else if ( yr >= zr ) { 1576 newspec[ii] += oldspec[j] * abs( yr - zl ) ; 1577 newflag[ii] = newflag[ii] || oldflag[j] ; 1578 wsum += abs( yr - zl ) ; 1579 } 1580 else { 1581 newspec[ii] += oldspec[j] * abs( dnu ) ; 1582 newflag[ii] = newflag[ii] || oldflag[j] ; 1583 wsum += abs( dnu ) ; 1584 ichan = j ; 1585 break ; 1586 } 1587 } 1588 else if ( yl > zr ) { 1589 if ( yr >= zr ) { 1590 newspec[ii] += oldspec[j] * abs( yr - yl ) ; 1591 newflag[ii] = newflag[ii] || oldflag[j] ; 1592 wsum += abs( yr - yl ) ; 1593 } 1594 else { 1595 newspec[ii] += oldspec[j] * abs( zr - yl ) ; 1596 newflag[ii] = newflag[ii] || oldflag[j] ; 1597 wsum += abs( zr - yl ) ; 1598 ichan = j ; 1599 break ; 1600 } 1601 } 1602 else { 1603 ichan = j - 1 ; 1604 break ; 1605 } 1606 } 1607 newspec[ii] /= wsum ; 1608 wsum = 0.0 ; 1609 } 1610 } 1611 // * ichan = 0 1612 // ***/ 1613 // //ofs << "olddnu = " << olddnu << ", dnu = " << dnu << endl ; 1614 // pile += dnu ; 1615 // wedge = olddnu * ( refChan + 1 ) ; 1616 // while ( wedge < pile ) { 1617 // newspec[0] += olddnu * oldspec[refChan] ; 1618 // newflag[0] = newflag[0] || oldflag[refChan] ; 1619 // //ofs << "channel " << refChan << " is included in new channel 0" << endl ; 1620 // refChan++ ; 1621 // wedge += olddnu ; 1622 // wsum += olddnu ; 1623 // //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ; 1624 // } 1625 // frac = ( wedge - pile ) / olddnu ; 1626 // wsum += ( 1.0 - frac ) * olddnu ; 1627 // newspec[0] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ; 1628 // newflag[0] = newflag[0] || oldflag[refChan] ; 1629 // //ofs << "channel " << refChan << " is partly included in new channel 0" << " with fraction of " << ( 1.0 - frac ) << endl ; 1630 // //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ; 1631 // newspec[0] /= wsum ; 1632 // //ofs << "newspec[0] = " << newspec[0] << endl ; 1633 // //ofs << "wedge = " << wedge << ", pile = " << pile << endl ; 1634 1635 // /*** 1636 // * ichan = 1 - nChan-2 1637 // ***/ 1638 // for ( int ichan = 1 ; ichan < nChan - 1 ; ichan++ ) { 1639 // pile += dnu ; 1640 // newspec[ichan] += frac * olddnu * oldspec[refChan] ; 1641 // newflag[ichan] = newflag[ichan] || oldflag[refChan] ; 1642 // //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << frac << endl ; 1643 // refChan++ ; 1644 // wedge += olddnu ; 1645 // wsum = frac * olddnu ; 1646 // //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ; 1647 // while ( wedge < pile ) { 1648 // newspec[ichan] += olddnu * oldspec[refChan] ; 1649 // newflag[ichan] = newflag[ichan] || oldflag[refChan] ; 1650 // //ofs << "channel " << refChan << " is included in new channel " << ichan << endl ; 1651 // refChan++ ; 1652 // wedge += olddnu ; 1653 // wsum += olddnu ; 1654 // //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ; 1655 // } 1656 // frac = ( wedge - pile ) / olddnu ; 1657 // wsum += ( 1.0 - frac ) * olddnu ; 1658 // newspec[ichan] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ; 1659 // newflag[ichan] = newflag[ichan] || oldflag[refChan] ; 1660 // //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << ( 1.0 - frac ) << endl ; 1661 // //ofs << "wedge = " << wedge << ", pile = " << pile << endl ; 1662 // //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ; 1663 // newspec[ichan] /= wsum ; 1664 // //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << endl ; 1665 // } 1666 1667 // /*** 1668 // * ichan = nChan-1 1669 // ***/ 1670 // // NOTE: Assumed that all spectra have the same bandwidth 1671 // pile += dnu ; 1672 // newspec[nChan-1] += frac * olddnu * oldspec[refChan] ; 1673 // newflag[nChan-1] = newflag[nChan-1] || oldflag[refChan] ; 1674 // //ofs << "channel " << refChan << " is partly included in new channel " << nChan-1 << " with fraction of " << frac << endl ; 1675 // refChan++ ; 1676 // wedge += olddnu ; 1677 // wsum = frac * olddnu ; 1678 // //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ; 1679 // for ( int jchan = refChan ; jchan < oldsize ; jchan++ ) { 1680 // newspec[nChan-1] += olddnu * oldspec[jchan] ; 1681 // newflag[nChan-1] = newflag[nChan-1] || oldflag[jchan] ; 1682 // wsum += olddnu ; 1683 // //ofs << "channel " << jchan << " is included in new channel " << nChan-1 << " with fraction of " << frac << endl ; 1684 // //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ; 1685 // } 1686 // //ofs << "wedge = " << wedge << ", pile = " << pile << endl ; 1687 // //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ; 1688 // newspec[nChan-1] /= wsum ; 1689 // //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << endl ; 1690 1691 // specCol_.put( irow, newspec ) ; 1692 // flagsCol_.put( irow, newflag ) ; 1693 1694 // // ofs.close() ; 1695 1696 1697 return ; 1698 } 1699 1144 1700 std::vector<float> Scantable::getWeather(int whichrow) const 1145 1701 { … … 1154 1710 1155 1711 } 1156 1712 //namespace asap
Note: See TracChangeset
for help on using the changeset viewer.