Line data    Source code

       1             : //# VisBuffer.tcc 
       2             : //# Copyright (C) 2011
       3             : //# Associated Universities, Inc. Washington DC, USA.
       4             : //#
       5             : //# This library is free software; you can redistribute it and/or modify
       6             : //# it under the terms of the GNU General Public License as published by
       7             : //# the Free Software Foundation; either version 2 of the License, or
       8             : //# (at your option) any later version.
       9             : //#
      10             : //# This library is distributed in the hope that it will be useful,
      11             : //# but WITHOUT ANY WARRANTY; without even the implied warranty of
      12             : //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      13             : //# GNU General Public License for more details.
      14             : //# 
      15             : //# You should have received a copy of the GNU General Public License
      16             : //# along with this library; if not, write to the Free Software
      17             : //# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
      18             : //#
      19             : //# Correspondence concerning AIPS++ should be addressed as follows:
      20             : //#        Internet email: aips2-request@nrao.edu.
      21             : //#        Postal address: AIPS++ Project Office
      22             : //#                        National Radio Astronomy Observatory
      23             : //#                        520 Edgemont Road
      24             : //#                        Charlottesville, VA 22903-2475 USA
      25             : //#
      26             : #include <msvis/MSVis/VisBuffer.h>
      27             : #include <casacore/casa/Arrays/Cube.h>
      28             : #include <cmath>
      29             : 
      30             : namespace casa { //# NAMESPACE CASA - BEGIN
      31             : template<class T>
      32          36 : void VisBuffer::chanAveVisCube(casacore::Cube<T>& data, casacore::Int nChanOut)
      33             : {
      34             :     using casacore::operator*;
      35          36 :   casacore::IPosition csh(data.shape());
      36          36 :   casacore::Int nChan0 = csh(1);
      37             : 
      38          36 :   if(nChan0 < nChanOut)
      39             :     // It's possible that data has already been averaged.  I could try
      40             :     // refilling data if I knew which column to use, but I don't.
      41             :     // Chuck it to the caller.
      42             :     throw(casacore::AipsError("Can't average " + casacore::String(nChan0) + " channels to " +
      43           0 :                     casacore::String(nChanOut) + " channels!"));
      44             : 
      45          36 :   csh(1) = nChanOut;
      46             : 
      47          36 :   casacore::Vector<casacore::Int>& chans(channel());
      48          36 :   casacore::Bool areShifting = true;
      49          36 :   if(chans.nelements() > 0 && chans[0] == 0)
      50           0 :     areShifting = false;
      51          36 :   if(nChan0 == nChanOut && !areShifting)
      52           0 :     return;                     // No-op.
      53             : 
      54          72 :   casacore::Cube<T> newCube(csh);
      55          36 :   newCube = T(0.0);
      56          36 :   casacore::Int nCor = nCorr();
      57          36 :   casacore::Int ichan(0);
      58             : 
      59          36 :   const casacore::Bool doSpWt(visIter_p->existsWeightSpectrum());
      60             : 
      61             :   // Make sure weightSpectrum() is unaveraged.
      62          36 :   if(doSpWt && (areShifting || weightSpectrum().shape()(1) < nChan0))
      63           0 :     fillWeightSpectrum();
      64             : 
      65          72 :   casacore::Vector<casacore::Double> totwt(nCor);
      66       19368 :   for(casacore::Int row = 0; row < nRow(); ++row){
      67       19332 :     if(!flagRow()(row)){
      68       19332 :       ichan = 0;
      69       96660 :       for(casacore::Int ochan = 0; ochan < nChanOut; ++ochan){
      70       77328 :         totwt = 0;
      71      154656 :         while(chans[ichan] >= chanAveBounds_p(ochan, 0) &&
      72      154656 :               chans[ichan] <= chanAveBounds_p(ochan, 1) &&
      73             :               ichan < nChan0){
      74      386640 :           for(casacore::Int icor = 0; icor < nCor; ++icor){
      75      309312 :             if(!flagCube()(icor, ichan, row)){
      76      302076 :               casacore::Double wt = 1.0;
      77             :               
      78      302076 :               if(doSpWt){ 
      79           0 :                 wt = weightSpectrum()(icor, ichan, row);
      80           0 :                 newCube(icor, ochan, row) += wt * data(icor, ichan, row);
      81             :               }
      82             :               else      // Mathematically equiv. but cheaper.
      83      302076 :                 newCube(icor, ochan, row) += data(icor, ichan, row);
      84             : 
      85      302076 :               totwt[icor] += wt;
      86             :             }
      87             :           }
      88       77328 :           ++ichan;
      89             :         }
      90      386640 :         for(casacore::Int icor = 0; icor < nCor; ++icor){
      91      309312 :           if(totwt[icor] > 0.0)
      92             :             //newCube(icor, ochan, row) *= T(1.0 / totwt[icor]);
      93      302076 :             newCube(icor, ochan, row) *= 1.0 / totwt[icor];
      94             :           // else...don't do any flag setting yet, flag needs to stay pristine
      95             :           // for now.  It will be updated by chanAveFlagCube().
      96             :         }
      97             :       }
      98             :     }
      99             :   }
     100          36 :   data.reference(newCube);        // Install averaged info
     101             : }
     102             : 
     103             : template<class T>
     104           0 : void VisBuffer::chanAccCube(casacore::Cube<T>& cube, casacore::Int nChanOut)
     105             : {
     106           0 :   casacore::IPosition csh(cube.shape());
     107           0 :   casacore::Int nChan0 = csh(1);
     108           0 :   csh(1) = nChanOut;
     109             : 
     110           0 :   if(nChan0 < nChanOut)
     111             :     // It's possible that cube has already been squeezed.  I could try
     112             :     // refilling data if I knew which column to use, but I don't.
     113             :     // Chuck it to the caller.
     114             :     throw(casacore::AipsError("Can't accumulate " + casacore::String(nChan0) + " channels to " +
     115           0 :                     casacore::String(nChanOut) + " channels!"));
     116           0 :   if(nChan0 == nChanOut)
     117           0 :     return;                     // No-op.
     118             : 
     119           0 :   casacore::Vector<casacore::Int>& chans(channel());
     120           0 :   casacore::Cube<T> newCube(csh);
     121           0 :   newCube = T(0.0);
     122           0 :   casacore::Int nCor = nCorr();
     123           0 :   casacore::Int ichan(0);
     124             : 
     125           0 :   for(casacore::Int row = 0; row < nRow(); ++row){
     126           0 :     if(!flagRow()(row)){
     127           0 :       ichan = 0;
     128           0 :       for(casacore::Int ochan = 0; ochan < nChanOut; ++ochan){
     129           0 :         while(chans[ichan] >= chanAveBounds_p(ochan, 0) &&
     130           0 :               chans[ichan] <= chanAveBounds_p(ochan, 1) &&
     131             :               ichan < nChan0){
     132           0 :           for(casacore::Int icor = 0; icor < nCor; ++icor)
     133           0 :             if(!flagCube()(icor, ichan, row))
     134           0 :               newCube(icor, ochan, row) += cube(icor, ichan, row);
     135           0 :           ++ichan;
     136             :         }
     137             :       }
     138             :     }
     139             :   }
     140           0 :   cube.reference(newCube);        // Install accumulated cube
     141             : }
     142             : 
     143             : } //# NAMESPACE CASA - END
     144             : 
     145             :