Changeset 1610 for branches/alma/src


Ignore:
Timestamp:
08/03/09 11:34:53 (15 years ago)
Author:
Takeshi Nakazato
Message:

New Development: No

JIRA Issue: No

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Update of spectral regrigging algorithm.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/alma/src/Scantable.cpp

    r1603 r1610  
    13081308  int oldsize = abcissa.size() ;
    13091309  double olddnu = abcissa[1] - abcissa[0] ;
    1310   int refChan = 0 ;
    1311   double frac = 0.0 ;
    1312   double wedge = 0.0 ;
    1313   double pile = 0.0 ;
     1310  //int refChan = 0 ;
     1311  //double frac = 0.0 ;
     1312  //double wedge = 0.0 ;
     1313  //double pile = 0.0 ;
     1314  int ichan = 0 ;
    13141315  double wsum = 0.0 ;
    1315   /***
    1316    * ichan = 0
    1317    ***/
    1318   //ofs << "olddnu = " << olddnu << ", dnu = " << dnu << endl ;
    1319   pile += dnu ;
    1320   wedge = olddnu * ( refChan + 1 ) ;
    1321   while ( wedge < pile ) {
    1322     newspec[0] += olddnu * oldspec[refChan] ;
    1323     newflag[0] = newflag[0] || oldflag[refChan] ;
    1324     //ofs << "channel " << refChan << " is included in new channel 0" << endl ;
    1325     refChan++ ;
    1326     wedge += olddnu ;
    1327     wsum += olddnu ;
    1328     //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
    1329   }
    1330   frac = ( wedge - pile ) / olddnu ;
    1331   wsum += ( 1.0 - frac ) * olddnu ;
    1332   newspec[0] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
    1333   newflag[0] = newflag[0] || oldflag[refChan] ;
    1334   //ofs << "channel " << refChan << " is partly included in new channel 0" << " with fraction of " << ( 1.0 - frac ) << endl ;
    1335   //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
    1336   newspec[0] /= wsum ;
    1337   //ofs << "newspec[0] = " << newspec[0] << endl ;
    1338   //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
    1339 
    1340   /***
    1341    * ichan = 1 - nChan-2
    1342    ***/
    1343   for ( int ichan = 1 ; ichan < nChan - 1 ; ichan++ ) {
    1344     pile += dnu ;
    1345     newspec[ichan] += frac * olddnu * oldspec[refChan] ;
    1346     newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
    1347     //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << frac << endl ;
    1348     refChan++ ;
    1349     wedge += olddnu ;
    1350     wsum = frac * olddnu ;
    1351     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
    1352     while ( wedge < pile ) {
    1353       newspec[ichan] += olddnu * oldspec[refChan] ;
    1354       newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
    1355       //ofs << "channel " << refChan << " is included in new channel " << ichan << endl ;
    1356       refChan++ ;
    1357       wedge += olddnu ;
    1358       wsum += olddnu ;
    1359       //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
    1360     }
    1361     frac = ( wedge - pile ) / olddnu ;
    1362     wsum += ( 1.0 - frac ) * olddnu ;
    1363     newspec[ichan] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
    1364     newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
    1365     //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << ( 1.0 - frac ) << endl ;
    1366     //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
    1367     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
    1368     newspec[ichan] /= wsum ;
    1369     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << endl ;
    1370   }
    1371 
    1372   /***
    1373    * ichan = nChan-1
    1374    ***/
    1375   // NOTE: Assumed that all spectra have the same bandwidth
    1376   pile += dnu ;
    1377   newspec[nChan-1] += frac * olddnu * oldspec[refChan] ;
    1378   newflag[nChan-1] = newflag[nChan-1] || oldflag[refChan] ;
    1379   //ofs << "channel " << refChan << " is partly included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
    1380   refChan++ ;
    1381   wedge += olddnu ;
    1382   wsum = frac * olddnu ;
    1383   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
    1384   for ( int jchan = refChan ; jchan < oldsize ; jchan++ ) {
    1385     newspec[nChan-1] += olddnu * oldspec[jchan] ;
    1386     newflag[nChan-1] = newflag[nChan-1] || oldflag[jchan] ;
    1387     wsum += olddnu ;
    1388     //ofs << "channel " << jchan << " is included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
    1389     //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
    1390   }
    1391   //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
    1392   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
    1393   newspec[nChan-1] /= wsum ;
    1394   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << endl ;
     1316  Vector<Float> z( nChan ) ;
     1317  z[0] = abcissa[0] - 0.5 * olddnu + 0.5 * dnu ;
     1318  for ( int ii = 1 ; ii < nChan ; ii++ )
     1319    z[ii] = z[ii-1] + dnu ;
     1320  Vector<Float> zi( nChan+1 ) ;
     1321  Vector<Float> yi( oldsize + 1 ) ;
     1322  zi[0] = z[0] - 0.5 * dnu ;
     1323  zi[1] = z[0] + 0.5 * dnu ;
     1324  for ( int ii = 2 ; ii < nChan ; ii++ )
     1325    zi[ii] = zi[ii-1] + dnu ;
     1326  zi[nChan] = z[nChan-1] + 0.5 * dnu ;
     1327  yi[0] = abcissa[0] - 0.5 * olddnu ;
     1328  yi[1] = abcissa[1] + 0.5 * olddnu ;
     1329  for ( int ii = 2 ; ii < oldsize ; ii++ )
     1330    yi[ii] = abcissa[ii-1] + olddnu ;
     1331  yi[oldsize] = abcissa[oldsize-1] + 0.5 * olddnu ;
     1332  if ( dnu > 0.0 ) {
     1333    for ( int ii = 0 ; ii < nChan ; ii++ ) {
     1334      double zl = zi[ii] ;
     1335      double zr = zi[ii+1] ;
     1336      for ( int j = ichan ; j < oldsize ; j++ ) {
     1337        double yl = yi[j] ;
     1338        double yr = yi[j+1] ;
     1339        if ( yl <= zl ) {
     1340          if ( yr <= zl ) {
     1341            continue ;
     1342          }
     1343          else if ( yr <= zr ) {
     1344            newspec[ii] += oldspec[j] * ( yr - zl ) ;
     1345            newflag[ii] = newflag[ii] || oldflag[j] ;
     1346            wsum += ( yr - zl ) ;
     1347          }
     1348          else {
     1349            newspec[ii] += oldspec[j] * dnu ;
     1350            newflag[ii] = newflag[ii] || oldflag[j] ;
     1351            wsum += dnu ;
     1352            ichan = j ;
     1353            break ;
     1354          }
     1355        }
     1356        else if ( yl < zr ) {
     1357          if ( yr <= zr ) {
     1358              newspec[ii] += oldspec[j] * ( yr - yl ) ;
     1359              newflag[ii] = newflag[ii] || oldflag[j] ;
     1360              wsum += ( yr - yl ) ;
     1361          }
     1362          else {
     1363            newspec[ii] += oldspec[j] * ( zr - yl ) ;
     1364            newflag[ii] = newflag[ii] || oldflag[j] ;
     1365            wsum += ( zr - yl ) ;
     1366            ichan = j ;
     1367            break ;
     1368          }
     1369        }
     1370        else {
     1371          ichan = j - 1 ;
     1372          break ;
     1373        }
     1374      }
     1375      newspec[ii] /= wsum ;
     1376      wsum = 0.0 ;
     1377    }
     1378  }
     1379  else if ( dnu < 0.0 ) {
     1380    for ( int ii = 0 ; ii < nChan ; ii++ ) {
     1381      double zl = zi[ii] ;
     1382      double zr = zi[ii+1] ;
     1383      for ( int j = ichan ; j < oldsize ; j++ ) {
     1384        double yl = yi[j] ;
     1385        double yr = yi[j+1] ;
     1386        if ( yl >= zl ) {
     1387          if ( yr >= zl ) {
     1388            continue ;
     1389          }
     1390          else if ( yr >= zr ) {
     1391            newspec[ii] += oldspec[j] * abs( yr - zl ) ;
     1392            newflag[ii] = newflag[ii] || oldflag[j] ;
     1393            wsum += abs( yr - zl ) ;
     1394          }
     1395          else {
     1396            newspec[ii] += oldspec[j] * abs( dnu ) ;
     1397            newflag[ii] = newflag[ii] || oldflag[j] ;
     1398            wsum += abs( dnu ) ;
     1399            ichan = j ;
     1400            break ;
     1401          }
     1402        }
     1403        else if ( yl > zr ) {
     1404          if ( yr >= zr ) {
     1405            newspec[ii] += oldspec[j] * abs( yr - yl ) ;
     1406            newflag[ii] = newflag[ii] || oldflag[j] ;
     1407            wsum += abs( yr - yl ) ;
     1408          }
     1409          else {
     1410            newspec[ii] += oldspec[j] * abs( zr - yl ) ;
     1411            newflag[ii] = newflag[ii] || oldflag[j] ;
     1412            wsum += abs( zr - yl ) ;
     1413            ichan = j ;
     1414            break ;
     1415          }
     1416        }
     1417        else {
     1418          ichan = j - 1 ;
     1419          break ;
     1420        }
     1421      }
     1422      newspec[ii] /= wsum ;
     1423      wsum = 0.0 ;
     1424    }
     1425  }
     1426//    * ichan = 0
     1427//    ***/
     1428//   //ofs << "olddnu = " << olddnu << ", dnu = " << dnu << endl ;
     1429//   pile += dnu ;
     1430//   wedge = olddnu * ( refChan + 1 ) ;
     1431//   while ( wedge < pile ) {
     1432//     newspec[0] += olddnu * oldspec[refChan] ;
     1433//     newflag[0] = newflag[0] || oldflag[refChan] ;
     1434//     //ofs << "channel " << refChan << " is included in new channel 0" << endl ;
     1435//     refChan++ ;
     1436//     wedge += olddnu ;
     1437//     wsum += olddnu ;
     1438//     //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
     1439//   }
     1440//   frac = ( wedge - pile ) / olddnu ;
     1441//   wsum += ( 1.0 - frac ) * olddnu ;
     1442//   newspec[0] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
     1443//   newflag[0] = newflag[0] || oldflag[refChan] ;
     1444//   //ofs << "channel " << refChan << " is partly included in new channel 0" << " with fraction of " << ( 1.0 - frac ) << endl ;
     1445//   //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
     1446//   newspec[0] /= wsum ;
     1447//   //ofs << "newspec[0] = " << newspec[0] << endl ;
     1448//   //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
     1449
     1450//   /***
     1451//    * ichan = 1 - nChan-2
     1452//    ***/
     1453//   for ( int ichan = 1 ; ichan < nChan - 1 ; ichan++ ) {
     1454//     pile += dnu ;
     1455//     newspec[ichan] += frac * olddnu * oldspec[refChan] ;
     1456//     newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
     1457//     //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << frac << endl ;
     1458//     refChan++ ;
     1459//     wedge += olddnu ;
     1460//     wsum = frac * olddnu ;
     1461//     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
     1462//     while ( wedge < pile ) {
     1463//       newspec[ichan] += olddnu * oldspec[refChan] ;
     1464//       newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
     1465//       //ofs << "channel " << refChan << " is included in new channel " << ichan << endl ;
     1466//       refChan++ ;
     1467//       wedge += olddnu ;
     1468//       wsum += olddnu ;
     1469//       //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
     1470//     }
     1471//     frac = ( wedge - pile ) / olddnu ;
     1472//     wsum += ( 1.0 - frac ) * olddnu ;
     1473//     newspec[ichan] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
     1474//     newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
     1475//     //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << ( 1.0 - frac ) << endl ;
     1476//     //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
     1477//     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
     1478//     newspec[ichan] /= wsum ;
     1479//     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << endl ;
     1480//   }
     1481
     1482//   /***
     1483//    * ichan = nChan-1
     1484//    ***/
     1485//   // NOTE: Assumed that all spectra have the same bandwidth
     1486//   pile += dnu ;
     1487//   newspec[nChan-1] += frac * olddnu * oldspec[refChan] ;
     1488//   newflag[nChan-1] = newflag[nChan-1] || oldflag[refChan] ;
     1489//   //ofs << "channel " << refChan << " is partly included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
     1490//   refChan++ ;
     1491//   wedge += olddnu ;
     1492//   wsum = frac * olddnu ;
     1493//   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
     1494//   for ( int jchan = refChan ; jchan < oldsize ; jchan++ ) {
     1495//     newspec[nChan-1] += olddnu * oldspec[jchan] ;
     1496//     newflag[nChan-1] = newflag[nChan-1] || oldflag[jchan] ;
     1497//     wsum += olddnu ;
     1498//     //ofs << "channel " << jchan << " is included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
     1499//     //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
     1500//   }
     1501//   //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
     1502//   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
     1503//   newspec[nChan-1] /= wsum ;
     1504//   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << endl ;
    13951505 
    1396   specCol_.put( irow, newspec ) ;
    1397   flagsCol_.put( irow, newflag ) ;
    1398 
    1399   // ofs.close() ;
     1506//   specCol_.put( irow, newspec ) ;
     1507//   flagsCol_.put( irow, newflag ) ;
     1508
     1509//   // ofs.close() ;
     1510
    14001511
    14011512  return ;
Note: See TracChangeset for help on using the changeset viewer.