- Timestamp:
- 05/29/12 11:52:20 (13 years ago)
- Location:
- branches/hpc33/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/hpc33/src/STMath.cpp
r2550 r2551 4711 4711 } 4712 4712 4713 Vector<Float> STMath::getSpectrumFromTime( double reftime, 4714 const Vector<Double> &timeVec, 4715 const Matrix<Float>& spectra, 4716 string mode ) 4717 { 4718 LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ; 4719 Vector<Float> sp ; 4720 uInt ncol = spectra.ncolumn() ; 4721 4722 if ( ncol == 0 ) { 4723 os << LogIO::SEVERE << "No spectra in the input scantable. Return empty spectrum." << LogIO::POST ; 4724 return sp ; 4725 } 4726 else if ( ncol == 1 ) { 4727 //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ; 4728 //s->specCol_.get( 0, sp ) ; 4729 sp.reference( spectra.column( 0 ) ) ; 4730 return sp ; 4731 } 4732 else { 4733 vector<int> idx = getRowIdFromTime( reftime, timeVec ) ; 4734 if ( mode == "before" ) { 4735 int id = -1 ; 4736 if ( idx[0] != -1 ) { 4737 id = idx[0] ; 4738 } 4739 else if ( idx[1] != -1 ) { 4740 os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ; 4741 id = idx[1] ; 4742 } 4743 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ; 4744 //s->specCol_.get( id, sp ) ; 4745 //sp = s->getSpectrum( id ) ; 4746 sp.reference( spectra.column( id ) ) ; 4747 } 4748 else if ( mode == "after" ) { 4749 int id = -1 ; 4750 if ( idx[1] != -1 ) { 4751 id = idx[1] ; 4752 } 4753 else if ( idx[0] != -1 ) { 4754 os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ; 4755 id = idx[1] ; 4756 } 4757 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ; 4758 //s->specCol_.get( id, sp ) ; 4759 //sp = s->getSpectrum( id ) ; 4760 sp.reference( spectra.column( id ) ) ; 4761 } 4762 else if ( mode == "nearest" ) { 4763 int id = -1 ; 4764 if ( idx[0] == -1 ) { 4765 id = idx[1] ; 4766 } 4767 else if ( idx[1] == -1 ) { 4768 id = idx[0] ; 4769 } 4770 else if ( idx[0] == idx[1] ) { 4771 id = idx[0] ; 4772 } 4773 else { 4774 double t0 = timeVec[idx[0]] ; 4775 double t1 = timeVec[idx[1]] ; 4776 double tref = reftime ; 4777 if ( abs( t0 - tref ) > abs( t1 - tref ) ) { 4778 id = idx[1] ; 4779 } 4780 else { 4781 id = idx[0] ; 4782 } 4783 } 4784 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ; 4785 //s->specCol_.get( id, sp ) ; 4786 //sp = s->getSpectrum( id ) ; 4787 sp.reference( spectra.column( id ) ) ; 4788 } 4789 else if ( mode == "linear" ) { 4790 if ( idx[0] == -1 ) { 4791 // use after 4792 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ; 4793 int id = idx[1] ; 4794 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ; 4795 //s->specCol_.get( id, sp ) ; 4796 //sp = s->getSpectrum( id ) ; 4797 sp.reference( spectra.column( id ) ) ; 4798 } 4799 else if ( idx[1] == -1 ) { 4800 // use before 4801 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ; 4802 int id = idx[0] ; 4803 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ; 4804 //s->specCol_.get( id, sp ) ; 4805 //sp = s->getSpectrum( id ) ; 4806 sp.reference( spectra.column( id ) ) ; 4807 } 4808 else if ( idx[0] == idx[1] ) { 4809 // use before 4810 //os << "No need to interporate." << LogIO::POST ; 4811 int id = idx[0] ; 4812 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ; 4813 //s->specCol_.get( id, sp ) ; 4814 //sp = s->getSpectrum( id ) ; 4815 sp.reference( spectra.column( id ) ) ; 4816 } 4817 else { 4818 // do interpolation 4819 //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ; 4820 double t0 = timeVec[idx[0]] ; 4821 double t1 = timeVec[idx[1]] ; 4822 double tref = reftime ; 4823 //s->specCol_.get( idx[0], sp ) ; 4824 //sp = s->getSpectrum( idx[0] ) ; 4825 //vector<float> sp1 = s->getSpectrum( idx[1] ) ; 4826 sp = spectra.column( idx[0] ).copy() ; 4827 //Vector<Float> sp1 = s->specCol_( idx[1] ) ; 4828 Vector<Float> sp1( spectra.column( idx[1] ) ) ; 4829 double tfactor = ( tref - t0 ) / ( t1 - t0 ) ; 4830 for ( unsigned int i = 0 ; i < sp.size() ; i++ ) { 4831 sp[i] = ( sp1[i] - sp[i] ) * tfactor + sp[i] ; 4832 } 4833 } 4834 } 4835 else { 4836 os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ; 4837 } 4838 return sp ; 4839 } 4840 } 4841 4713 4842 void STMath::getSpectrumFromTime2( Vector<Float> &sp, 4714 4843 double reftime, … … 4895 5024 } 4896 5025 4897 vector<int> STMath::getRowIdFromTime( double reftime, Vector<Double> &t )5026 vector<int> STMath::getRowIdFromTime( double reftime, const Vector<Double> &t ) 4898 5027 { 4899 5028 // double reft = reftime ; -
branches/hpc33/src/STMath.h
r2550 r2551 404 404 vector<float> getSpectrumFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode = "before" ) ; 405 405 casa::Vector<casa::Float> getSpectrumFromTime( double reftime, casa::Vector<casa::Double> &v, const casa::CountedPtr<Scantable>& s, string mode = "before" ) ; 406 casa::Vector<casa::Float> getSpectrumFromTime( double reftime, const casa::Vector<casa::Double> &v, const casa::Matrix<casa::Float>& sp, string mode = "before" ) ; 406 407 void getSpectrumFromTime2( casa::Vector<casa::Float> &sp, double reftime, casa::Vector<casa::Double> &v, const casa::CountedPtr<Scantable>& s, string mode = "before" ) ; 407 408 vector<float> getTcalFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ; 408 409 vector<float> getTsysFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ; 409 410 vector<int> getRowIdFromTime( string reftime, casa::CountedPtr<Scantable>& s ) ; 410 vector<int> getRowIdFromTime( double reftime, c asa::Vector<casa::Double>& t ) ;411 vector<int> getRowIdFromTime( double reftime, const casa::Vector<casa::Double>& t ) ; 411 412 412 413 // Chopper-Wheel type calibration
Note:
See TracChangeset
for help on using the changeset viewer.