Line data    Source code

       1             : //# WTerm.cc: implementation of WTerm
       2             : //# Copyright (C) 2007
       3             : //# Associated Universities, Inc. Washington DC, USA.
       4             : //#
       5             : //# This library is free software; you can redistribute it and/or modify it
       6             : //# under the terms of the GNU Library General Public License as published by
       7             : //# the Free Software Foundation; either version 2 of the License, or (at your
       8             : //# option) any later version.
       9             : //#
      10             : //# This library is distributed in the hope that it will be useful, but WITHOUT
      11             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      12             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
      13             : //# License for more details.
      14             : //#
      15             : //# You should have received a copy of the GNU Library General Public License
      16             : //# along with this library; if not, write to the Free Software Foundation,
      17             : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
      18             : //#
      19             : //# Correspondence concerning AIPS++ should be adressed as follows:
      20             : //#        Internet email: casa-feedback@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             : //#
      27             : //# $Id$
      28             : #include <synthesis/TransformMachines2/WTerm.h>
      29             : #include <synthesis/TransformMachines/SynthesisMath.h>
      30             : #ifdef _OPENMP
      31             : #include <omp.h>
      32             : #endif
      33             : 
      34             : 
      35             : using namespace casacore;
      36             : namespace casa { //# NAMESPACE CASA - BEGIN
      37             :   using namespace refim;
      38          30 :   void WTerm::applySky(Matrix<Complex>& screen, 
      39             :                        const Vector<Double>& sampling,
      40             :                        const Double wValue,
      41             :                        const Int inner)
      42             :   {
      43          30 :     Int convSize = screen.shape()(0);
      44          30 :     Double twoPiW=2.0*C::pi*Double(wValue);
      45             : #ifdef _OPENMP
      46          30 :     Int Nth=max(omp_get_max_threads()-2,1);
      47             : #endif
      48             : 
      49          30 :     if (!isNoOp())
      50             :       {
      51       61470 :         for (Int iy=-inner/2;iy<inner/2;iy++) 
      52             :           {
      53       61440 :             Double m=sampling(1)*Double(iy);
      54       61440 :             Double msq=m*m;
      55             : #ifdef _OPENMP
      56             : //#pragma omp parallel default(none) firstprivate(msq,iy) shared(screen, sampling, twoPiW,convSize) num_threads(Nth)
      57       61440 : #pragma omp parallel firstprivate(msq,iy) shared(twoPiW,convSize) num_threads(Nth)
      58             : #endif
      59             :     {
      60             : #ifdef _OPENMP
      61             : #pragma omp for
      62             : #endif
      63             :             for (Int ix=-inner/2;ix<inner/2;ix++) 
      64             :               {
      65             :                 Double l=sampling(0)*Double(ix);
      66             :                 Double rsq=l*l+msq;
      67             :                 if(rsq<1.0) 
      68             :                   {
      69             :                     Double phase=twoPiW*(sqrt(1.0-rsq)-1.0), sp, cp;
      70             :                     SINCOS(phase, sp, cp);
      71             :                     screen(ix+convSize/2,iy+convSize/2)*=Complex(cp, sp);
      72             :                   }
      73             :               }
      74             :     }
      75             :           }
      76             :       }
      77          30 :   }
      78             : 
      79          64 :   void WTerm::applySky(Matrix<Complex>& screen, 
      80             :                        const Int wPixel, 
      81             :                        const Vector<Double>& sampling,
      82             :                        const Double wScale,
      83             :                        const Int inner)
      84             :   {
      85             :     //UNUSED: Int convSize = screen.shape()(0);
      86          64 :     if(wPixel>0) 
      87             :       {
      88          30 :         Double wValue=(wPixel*wPixel)/wScale;
      89          30 :         applySky(screen, sampling, wValue, inner);
      90             :       }
      91          64 :   }
      92             : };