LCOV - code coverage report
Current view: top level - synthesis/TransformMachines2 - MosaicFT.cc (source / functions) Hit Total Coverage
Test: casacpp_coverage.info Lines: 0 944 0.0 %
Date: 2024-10-04 16:51:10 Functions: 0 33 0.0 %

          Line data    Source code
       1             : //# MosaicFT.cc: Implementation of MosaicFT class
       2             : //# Copyright (C) 2002-2016
       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 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 General Public
      13             : //# 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 Foundation,
      17             : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
      18             : //#
      19             : //# Correspondence concerning AIPS++ should be addressed 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             : //# $Id$
      27             : 
      28             : #include <casacore/casa/Quanta/UnitMap.h>
      29             : #include <casacore/casa/Quanta/MVTime.h>
      30             : #include <casacore/casa/Quanta/UnitVal.h>
      31             : #include <casacore/measures/Measures/Stokes.h>
      32             : #include <casacore/measures/Measures/UVWMachine.h>
      33             : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
      34             : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
      35             : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
      36             : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
      37             : #include <casacore/coordinates/Coordinates/Projection.h>
      38             : #include <casacore/ms/MeasurementSets/MSColumns.h>
      39             : #include <casacore/casa/BasicSL/Constants.h>
      40             : #include <casacore/scimath/Mathematics/FFTServer.h>
      41             : #include <synthesis/TransformMachines2/MosaicFT.h>
      42             : #include <synthesis/TransformMachines2/SimplePBConvFunc.h>
      43             : #include <synthesis/TransformMachines2/HetArrayConvFunc.h>
      44             : #include <synthesis/TransformMachines/PBMath.h>
      45             : #include <synthesis/TransformMachines2/VPSkyJones.h>
      46             : #include <casacore/scimath/Mathematics/RigidVector.h>
      47             : #include <msvis/MSVis/StokesVector.h>
      48             : #include <synthesis/TransformMachines/StokesImageUtil.h>
      49             : #include <msvis/MSVis/VisBuffer2.h>
      50             : #include <msvis/MSVis/VisibilityIterator2.h>
      51             : #include <casacore/images/Images/ImageInterface.h>
      52             : #include <casacore/images/Images/PagedImage.h>
      53             : #include <casacore/images/Images/SubImage.h>
      54             : #include <casacore/images/Regions/ImageRegion.h>
      55             : #include <casacore/images/Regions/WCBox.h>
      56             : #include <casacore/casa/Containers/Block.h>
      57             : #include <casacore/casa/Containers/Record.h>
      58             : #include <casacore/casa/Arrays/ArrayLogical.h>
      59             : #include <casacore/casa/Arrays/ArrayMath.h>
      60             : #include <casacore/casa/Arrays/Array.h>
      61             : #include <casacore/casa/Arrays/MaskedArray.h>
      62             : #include <casacore/casa/Arrays/Vector.h>
      63             : #include <casacore/casa/Arrays/Slice.h>
      64             : #include <casacore/casa/Arrays/Matrix.h>
      65             : #include <casacore/casa/Arrays/Cube.h>
      66             : #include <casacore/casa/Arrays/MatrixIter.h>
      67             : #include <casacore/casa/BasicSL/String.h>
      68             : #include <casacore/casa/Utilities/Assert.h>
      69             : #include <casacore/casa/Exceptions/Error.h>
      70             : #include <casacore/lattices/Lattices/ArrayLattice.h>
      71             : #include <casacore/lattices/Lattices/SubLattice.h>
      72             : #include <casacore/lattices/LRegions/LCBox.h>
      73             : #include <casacore/lattices/LEL/LatticeExpr.h>
      74             : #include <casacore/lattices/Lattices/LatticeCache.h>
      75             : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
      76             : #include <casacore/lattices/Lattices/LatticeIterator.h>
      77             : #include <casacore/lattices/Lattices/LatticeStepper.h>
      78             : #include <casacore/casa/Utilities/CompositeNumber.h>
      79             : #include <casacore/casa/OS/Timer.h>
      80             : #include <casacore/casa/OS/HostInfo.h>
      81             : #include <sstream>
      82             : #ifdef _OPENMP
      83             : #include <omp.h>
      84             : #endif
      85             : using namespace casacore;
      86             : namespace casa { //# NAMESPACE CASA - BEGIN
      87             : namespace refim {//# namespace for imaging refactor
      88             : using namespace casacore;
      89             : using namespace casa;
      90             : using namespace casacore;
      91             : using namespace casa::refim;
      92             : 
      93           0 :   MosaicFT::MosaicFT(SkyJones* sj, MPosition mloc, String stokes,
      94             :                    Long icachesize, Int itilesize, 
      95           0 :                      Bool usezero, Bool useDoublePrec, Bool useConjConvFunc, Bool usePointing)
      96           0 :   : FTMachine(), sj_p(sj),
      97           0 :     imageCache(0),  cachesize(icachesize), tilesize(itilesize), gridder(nullptr),
      98           0 :     isTiled(false),
      99           0 :     maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
     100           0 :     mspc(0), msac(0), pointingToImage(0), usezero_p(usezero), convSampling(1),
     101           0 :     skyCoverage_p( ), machineName_p("MosaicFT"), stokes_p(stokes), useConjConvFunc_p(useConjConvFunc), usePointingTable_p(usePointing),timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0)
     102             : {
     103           0 :   convSize=0;
     104           0 :   lastIndex_p=0;
     105           0 :   doneWeightImage_p=false;
     106           0 :   convWeightImage_p=0;
     107           0 :   pbConvFunc_p=new SimplePBConvFunc();
     108             :     
     109           0 :   mLocation_p=mloc;
     110           0 :   useDoubleGrid_p=useDoublePrec;  
     111             :   // We should get rid of the ms dependence in the constructor
     112             :   // not used
     113           0 : }
     114             : 
     115           0 : MosaicFT::MosaicFT(const RecordInterface& stateRec)
     116           0 :   : FTMachine()
     117             : {
     118             :   // Construct from the input state record
     119           0 :   String error;
     120           0 :   if (!fromRecord(error, stateRec)) {
     121           0 :     throw (AipsError("Failed to create MosaicFT: " + error));
     122             :   };
     123           0 : }
     124             : 
     125             : //---------------------------------------------------------------------- 
     126           0 : MosaicFT& MosaicFT::operator=(const MosaicFT& other)
     127             : {
     128           0 :   if(this!=&other) {
     129             : 
     130             :     //Do the base parameters
     131           0 :     FTMachine::operator=(other);
     132             :      
     133           0 :     convSampling=other.convSampling;
     134           0 :     sj_p=other.sj_p;
     135           0 :     imageCache=other.imageCache;
     136           0 :     cachesize=other.cachesize;
     137           0 :     tilesize=other.tilesize;
     138           0 :     isTiled=other.isTiled;
     139             :     //lattice=other.lattice;
     140           0 :     lattice=0;
     141             :     // arrayLattice=other.arrayLattice;
     142             :     // weightLattice=other.weightLattice;
     143             :     //if(arrayLattice) delete arrayLattice;
     144           0 :     arrayLattice=0;
     145             :     //if(weightLattice) delete weightLattice;
     146           0 :     weightLattice=0;
     147           0 :     maxAbsData=other.maxAbsData;
     148           0 :     centerLoc=other.centerLoc;
     149           0 :     offsetLoc=other.offsetLoc;
     150           0 :     pointingToImage=other.pointingToImage;
     151           0 :     usezero_p=other.usezero_p;
     152           0 :     doneWeightImage_p=other.doneWeightImage_p;
     153           0 :     pbConvFunc_p=other.pbConvFunc_p;
     154           0 :     stokes_p=other.stokes_p;
     155           0 :     if(!other.skyCoverage_p.null())
     156           0 :       skyCoverage_p=other.skyCoverage_p;
     157             :     else
     158           0 :       skyCoverage_p=0;
     159           0 :     if(other.convWeightImage_p !=0)
     160           0 :       convWeightImage_p=(TempImage<Complex> *)other.convWeightImage_p->cloneII();
     161             :     else
     162           0 :       convWeightImage_p=0;
     163           0 :     if(other.gridder==0)
     164           0 :       gridder.reset(nullptr);
     165             :     else{
     166           0 :       uvScale=other.uvScale;
     167           0 :       uvOffset=other.uvOffset;
     168           0 :       gridder.reset(new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
     169           0 :                                                          uvScale, uvOffset,
     170           0 :                                                          "SF"));
     171             :           
     172             :     }
     173           0 :     useConjConvFunc_p=other.useConjConvFunc_p;
     174           0 :     usePointingTable_p=other.usePointingTable_p;
     175           0 :     timemass_p=other.timemass_p;
     176           0 :     timegrid_p=other.timegrid_p;
     177           0 :     timedegrid_p=other.timedegrid_p;
     178             :   };
     179           0 :   return *this;
     180             : };
     181             : 
     182             : //----------------------------------------------------------------------
     183           0 :   MosaicFT::MosaicFT(const MosaicFT& other): FTMachine(),machineName_p("MosaicFT")
     184             : {
     185           0 :   operator=(other);
     186           0 : }
     187             : 
     188             : //----------------------------------------------------------------------
     189             : //void MosaicFT::setSharingFT(MosaicFT& otherFT){
     190             : //  otherFT_p=&otherFT;
     191             : //}
     192           0 :   void MosaicFT::init(const vi::VisBuffer2& /*vb*/) {
     193             :   
     194             :   /* if((image->shape().product())>cachesize) {
     195             :     isTiled=true;
     196             :   }
     197             :   else {
     198             :     isTiled=false;
     199             :   }
     200             :   */
     201             :   //For now only isTiled false works.
     202           0 :   isTiled=false;
     203           0 :   nx    = image->shape()(0);
     204           0 :   ny    = image->shape()(1);
     205           0 :   npol  = image->shape()(2);
     206           0 :   nchan = image->shape()(3);
     207             : 
     208             : 
     209             :   //  if(skyCoverage_p==0){
     210             :   //    Double memoryMB=HostInfo::memoryTotal()/1024.0/(20.0);
     211             :   //    skyCoverage_p=new TempImage<Float> (IPosition(4,nx,ny,1,1),
     212             :   //                                    image->coordinates(), memoryMB);
     213             :     //Setting it to zero
     214             : //   (*skyCoverage_p)-=(*skyCoverage_p);
     215             : //  }
     216           0 :   sumWeight.resize(npol, nchan);
     217             :   
     218           0 :   convSupport=0;
     219             : 
     220           0 :   uvScale.resize(2);
     221           0 :   uvScale=0.0;
     222           0 :   uvScale(0)=Float(nx)*image->coordinates().increment()(0); 
     223           0 :   uvScale(1)=Float(ny)*image->coordinates().increment()(1); 
     224             :     
     225           0 :   uvOffset.resize(2);
     226           0 :   uvOffset(0)=nx/2;
     227           0 :   uvOffset(1)=ny/2;
     228             :   
     229           0 :   gridder.reset(new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
     230           0 :                                                      uvScale, uvOffset,
     231           0 :                                                      "SF"));
     232             : 
     233             :   // Set up image cache needed for gridding. 
     234           0 :   if(imageCache) delete imageCache; imageCache=0;
     235             :   /*
     236             :   if(isTiled) {
     237             :     Float tileOverlap=0.5;
     238             :     tilesize=min(256,tilesize);
     239             :     IPosition tileShape=IPosition(4,tilesize,tilesize,npol,nchan);
     240             :     Vector<Float> tileOverlapVec(4);
     241             :     tileOverlapVec=0.0;
     242             :     tileOverlapVec(0)=tileOverlap;
     243             :     tileOverlapVec(1)=tileOverlap;
     244             :     Int tmpCacheVal=static_cast<Int>(cachesize);
     245             :     imageCache=new LatticeCache <Complex> (*image, tmpCacheVal, tileShape, 
     246             :                                            tileOverlapVec,
     247             :                                            (tileOverlap>0.0));
     248             :     
     249             :   }
     250             :   */
     251           0 : }
     252             : 
     253             : // This is nasty, we should use CountedPointers here.
     254           0 : MosaicFT::~MosaicFT() {
     255           0 :   if(imageCache) delete imageCache; imageCache=0;
     256             :   //  if(arrayLattice) delete arrayLattice; arrayLattice=0;
     257           0 : }
     258             : 
     259             : 
     260           0 : void MosaicFT::setConvFunc(CountedPtr<SimplePBConvFunc>& pbconvFunc){
     261             : 
     262             : 
     263           0 :   pbConvFunc_p=pbconvFunc;
     264             :   
     265             : 
     266           0 : }
     267             : 
     268           0 : CountedPtr<SimplePBConvFunc>& MosaicFT::getConvFunc(){
     269           0 :   return pbConvFunc_p;
     270             : }
     271             : 
     272           0 : void MosaicFT::findConvFunction(const ImageInterface<Complex>& iimage,
     273             :                                 const vi::VisBuffer2& vb, const Matrix<Double>& /*rotateduvw*/) {
     274             :   
     275             :   
     276             :   //oversample if image is small
     277             :   //But not more than 5000 pixels
     278           0 :   convSampling=(max(nx, ny) < 50) ? 100: Int(ceil(5000.0/max(nx, ny)));
     279           0 :   if(convSampling <1) 
     280           0 :     convSampling=1;
     281           0 :   if(pbConvFunc_p.null())
     282           0 :     pbConvFunc_p=new SimplePBConvFunc();
     283           0 :   if(sj_p)
     284           0 :     pbConvFunc_p->setSkyJones(sj_p.get());
     285             :   ////TEST for HetArray only for now
     286           0 :   if(pbConvFunc_p->name()=="HetArrayConvFunc"){
     287           0 :     if(convSampling <10) 
     288           0 :       convSampling=10;
     289           0 :     AipsrcValue<Int>::find (convSampling, "mosaic.oversampling", 10);
     290             :   }
     291           0 :   if((pbConvFunc_p->getVBUtil()).null()){
     292           0 :     if(vbutil_p.null()){
     293           0 :         vbutil_p=new VisBufferUtil(vb);
     294             :     }
     295           0 :     pbConvFunc_p->setVBUtil(vbutil_p);
     296             :   }
     297             :   //cerr << "NELEMS " << interpVisFreq_p.nelements() << "  lsr " << lsrFreq_p.nelements() << endl;
     298           0 :   pbConvFunc_p->findConvFunction(iimage, vb, convSampling, interpVisFreq_p, convFunc, weightConvFunc_p, convSizePlanes_p, convSupportPlanes_p,
     299           0 :                                  convPolMap_p, convChanMap_p, convRowMap_p, (useConjConvFunc_p && !toVis_p), MVDirection(-(movingDirShift_p.getAngle())), fixMovingSource_p);
     300             : 
     301             :   // cerr << "MAX of convFunc " << max(abs(convFunc)) << endl;
     302             :   //For now only use one size and support
     303           0 :   if(convSizePlanes_p.nelements() ==0)
     304           0 :     convSize=0;
     305             :   else
     306           0 :     convSize=max(convSizePlanes_p);
     307           0 :   if(convSupportPlanes_p.nelements() ==0)
     308           0 :     convSupport=0;
     309             :   else
     310           0 :     convSupport=max(convSupportPlanes_p);
     311             :                                  
     312           0 : }
     313             : 
     314           0 : void MosaicFT::initializeToVis(ImageInterface<Complex>& iimage,
     315             :                                const vi::VisBuffer2& vb)
     316             : {
     317           0 :   image=&iimage;
     318           0 :   toVis_p=true;
     319           0 :   ok();
     320             :   
     321             :   //  if(convSize==0) {
     322           0 :     init(vb);
     323             :     
     324             :     //  }
     325             :   
     326             :   // Initialize the maps for polarization and channel. These maps
     327             :   // translate visibility indices into image indices
     328           0 :   initMaps(vb);
     329           0 :   pbConvFunc_p->setVBUtil(vbutil_p);
     330           0 :   pbConvFunc_p->setUsePointing(usePointingTable_p);
     331             :  //make sure we rotate the first field too
     332           0 :   lastFieldId_p=-1;
     333           0 :   phaseShifter_p=new UVWMachine(*uvwMachine_p);
     334             :   //This is needed here as we need to know the grid correction before FFTing 
     335           0 :   findConvFunction(*image, vb, vb.uvw());
     336             :   
     337           0 :   prepGridForDegrid();
     338             : 
     339           0 : }
     340             : 
     341             : 
     342           0 : void MosaicFT::prepGridForDegrid(){
     343             : 
     344             :   //For now isTiled=false
     345           0 :   isTiled=false;
     346           0 :   nx    = image->shape()(0);
     347           0 :   ny    = image->shape()(1);
     348           0 :   npol  = image->shape()(2);
     349           0 :   nchan = image->shape()(3);
     350             : 
     351           0 :   IPosition gridShape(4, nx, ny, npol, nchan);
     352           0 :   griddedData.resize(gridShape);
     353           0 :   griddedData=Complex(0.0);
     354             :   
     355           0 :   IPosition stride(4, 1);
     356           0 :   IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
     357           0 :                 (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
     358           0 :   IPosition trc(blc+image->shape()-stride);
     359             :     
     360           0 :   IPosition start(4, 0);
     361           0 :   griddedData(blc, trc) = image->getSlice(start, image->shape());
     362             :   
     363           0 :   image->clearCache();
     364             :     //if(arrayLattice) delete arrayLattice; arrayLattice=0;
     365           0 :   arrayLattice = new ArrayLattice<Complex>(griddedData);
     366           0 :   lattice=arrayLattice;
     367             :    {///UnDo the grid correction
     368           0 :       Int inx = lattice->shape()(0);
     369             :       //Int iny = lattice->shape()(1);
     370           0 :       Vector<Complex> correction(inx);
     371           0 :       correction=Complex(1.0, 0.0);
     372             : 
     373             :       // Int npixCorr= max(nx,ny);
     374           0 :       Vector<Float> sincConvX(nx);
     375           0 :       for (Int ix=0;ix<nx;ix++) {
     376           0 :         Float x=C::pi*Float(ix-nx/2)/(Float(nx)*Float(convSampling));
     377           0 :         if(ix==nx/2) {
     378           0 :           sincConvX(ix)=1.0;
     379             :         }
     380             :         else {
     381           0 :           sincConvX(ix)=sin(x)/x;
     382             :         }
     383             :       }
     384           0 :       Vector<Float> sincConvY(ny);
     385           0 :       for (Int ix=0;ix<ny;ix++) {
     386           0 :         Float x=C::pi*Float(ix-ny/2)/(Float(ny)*Float(convSampling));
     387           0 :         if(ix==ny/2) {
     388           0 :           sincConvY(ix)=1.0;
     389             :         }
     390             :         else {
     391           0 :           sincConvY(ix)=sin(x)/x;
     392             :         }
     393             :       }
     394             :  
     395           0 :        IPosition cursorShape(4, inx, 1, 1, 1);
     396           0 :       IPosition axisPath(4, 0, 1, 2, 3);
     397           0 :       LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
     398           0 :       LatticeIterator<Complex> lix(*lattice, lsx);
     399           0 :       for(lix.reset();!lix.atEnd();lix++) {
     400             :         //Int pol=lix.position()(2);
     401             :         //Int chan=lix.position()(3);
     402             :         
     403           0 :         Int iy=lix.position()(1);
     404             :         //gridder->correctX1D(correction,iy);
     405           0 :         for (Int ix=0;ix<nx;ix++) {
     406           0 :           correction(ix)=(sincConvX(ix)*sincConvY(iy));
     407             :         }
     408           0 :         lix.rwVectorCursor()*=correction;
     409             :         
     410             :       }
     411           0 :   }
     412             :     
     413             :   
     414           0 :   logIO() << LogIO::DEBUGGING << "Starting FFT of image" << LogIO::POST;
     415             :    // Now do the FFT2D in place
     416           0 :   ft_p.c2cFFT(*lattice);
     417             :   ///////////////////////
     418             :   /*{
     419             :     CoordinateSystem ftCoords(image->coordinates());
     420             :     Int directionIndex=ftCoords.findCoordinate(Coordinate::DIRECTION);
     421             :     DirectionCoordinate dc=ftCoords.directionCoordinate(directionIndex);
     422             :     Vector<Bool> axes(2); axes(0)=true;axes(1)=true;
     423             :     Vector<Int> shape(2); shape(0)=griddedData.shape()(0) ;shape(1)=griddedData.shape()(1);
     424             :     Coordinate* ftdc=dc.makeFourierCoordinate(axes,shape);
     425             :     ftCoords.replaceCoordinate(*ftdc, directionIndex);
     426             :     delete ftdc; ftdc=0;
     427             :     PagedImage<Float> thisScreen(griddedData.shape(), ftCoords, String("MODEL_GRID_VIS"));
     428             :     thisScreen.put(amplitude(griddedData));
     429             :     }*/
     430             :   ////////////////////////
     431           0 :   logIO() << LogIO::DEBUGGING << "Finished FFT" << LogIO::POST;
     432             : 
     433             : 
     434             : 
     435           0 : }
     436             : 
     437             : 
     438           0 : void MosaicFT::finalizeToVis()
     439             : {
     440           0 :   logIO() << LogOrigin("MosaicFT", "finalizeToVis")  << LogIO::NORMAL;
     441           0 :   logIO()<< LogIO::NORMAL2 << "Time degrid " << timedegrid_p << LogIO::POST;
     442           0 :   timedegrid_p=0.0;
     443             :   
     444           0 :   if(!arrayLattice.null()) arrayLattice=0;
     445           0 :   if(!lattice.null()) lattice=0;
     446           0 :   griddedData.resize();
     447             : 
     448             :   /*
     449             :   if(isTiled) {
     450             :     
     451             :     logIO() << LogOrigin("MosaicFT", "finalizeToVis")  << LogIO::NORMAL;
     452             :     
     453             :     AlwaysAssert(imageCache, AipsError);
     454             :     AlwaysAssert(image, AipsError);
     455             :     ostringstream o;
     456             :     imageCache->flush();
     457             :     imageCache->showCacheStatistics(o);
     458             :     logIO() << o.str() << LogIO::POST;
     459             :   }
     460             :   */
     461           0 :   if(pointingToImage)
     462           0 :     delete pointingToImage;
     463           0 :   pointingToImage=0;
     464           0 : }
     465             : 
     466             : 
     467             : // Initialize the FFT to the Sky. Here we have to setup and initialize the
     468             : // grid. 
     469             : 
     470             : 
     471             : 
     472             : 
     473             : 
     474           0 : void MosaicFT::initializeToSky(ImageInterface<Complex>& iimage,
     475             :                                Matrix<Float>& weight,
     476             :                                const vi::VisBuffer2& vb)
     477             : {
     478             :   // image always points to the image
     479           0 :   image=&iimage;
     480           0 :   toVis_p=False;
     481             :   //  if(convSize==0) {
     482           0 :     init(vb);
     483             :     
     484             :     //  }
     485             :   
     486             :   // Initialize the maps for polarization and channel. These maps
     487             :   // translate visibility indices into image indices
     488           0 :   initMaps(vb);
     489           0 :   pbConvFunc_p->setVBUtil(vbutil_p);
     490           0 :   pbConvFunc_p->setUsePointing(usePointingTable_p);
     491             :   //make sure we rotate the first field too
     492           0 :   lastFieldId_p=-1;
     493           0 :   phaseShifter_p=new UVWMachine(*uvwMachine_p);
     494             :   //findConvFunction(*image, vb);
     495             :   /*if((image->shape().product())>cachesize) {
     496             :     isTiled=true;
     497             :   }
     498             :   else {
     499             :     isTiled=false;
     500             :   }
     501             :   */
     502             :   //For now isTiled has to be false
     503           0 :   isTiled=false;
     504           0 :   nx    = image->shape()(0);
     505           0 :   ny    = image->shape()(1);
     506           0 :   npol  = image->shape()(2);
     507           0 :   nchan = image->shape()(3);
     508             : 
     509           0 :   sumWeight=0.0;
     510           0 :   weight.resize(sumWeight.shape());
     511           0 :   weight=0.0;
     512             :   
     513           0 :   image->clearCache();
     514             :   // Initialize for in memory or to disk gridding. lattice will
     515             :   // point to the appropriate Lattice, either the ArrayLattice for
     516             :   // in memory gridding or to the image for to disk gridding.
     517             :   /*if(isTiled) {
     518             :     imageCache->flush();
     519             :     image->set(Complex(0.0));
     520             :     lattice=CountedPtr<Lattice<Complex> >(image, false);
     521             :     if( !doneWeightImage_p && (convWeightImage_p==0)){
     522             :       
     523             :       convWeightImage_p=new  TempImage<Complex> (iimage.shape(), 
     524             :                                                  iimage.coordinates());
     525             : 
     526             : 
     527             : 
     528             : 
     529             :       convWeightImage_p->set(Complex(0.0));
     530             :       weightLattice=convWeightImage_p;
     531             : 
     532             :     }
     533             :     }
     534             :     else*/ 
     535             :   {
     536           0 :     IPosition gridShape(4, nx, ny, npol, nchan);
     537           0 :     if(!useDoubleGrid_p) {
     538           0 :         griddedData.resize(gridShape);
     539           0 :         griddedData=Complex(0.0);
     540             :       }
     541             :     else {
     542           0 :       griddedData2.resize(gridShape);
     543           0 :       griddedData2=DComplex(0.0);
     544             :     }
     545             :     //if(arrayLattice) delete arrayLattice; arrayLattice=0;
     546             :     //arrayLattice = new ArrayLattice<Complex>(griddedData);
     547             :     //lattice=arrayLattice;
     548             :       
     549           0 :     if( !doneWeightImage_p && (convWeightImage_p==0)){
     550             :      
     551             :       
     552             :  
     553           0 :       convWeightImage_p=new  TempImage<Complex> (iimage.shape(), 
     554           0 :                                                  iimage.coordinates());
     555           0 :       griddedWeight.resize(gridShape);
     556             :       /*IPosition stride(4, 1);
     557             :       IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
     558             :                     (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
     559             :       IPosition trc(blc+image->shape()-stride);
     560             :       
     561             :       griddedWeight(blc, trc).set(Complex(0.0));
     562             :       */
     563           0 :       if(useDoubleGrid_p){
     564           0 :         griddedWeight2.resize(gridShape);
     565           0 :         griddedWeight2=DComplex(0.0);
     566             :       }
     567             :       else{
     568           0 :         griddedWeight=Complex(0.0);
     569             :       }
     570             :       //if(weightLattice) delete weightLattice; weightLattice=0;
     571           0 :       weightLattice = new ArrayLattice<Complex>(griddedWeight);
     572             : 
     573             :     }
     574             : 
     575           0 :   }
     576             : 
     577             :   //cerr << "initializetosky lastfield " << lastFieldId_p << endl;
     578             :   // AlwaysAssert(lattice, AipsError);
     579             :   
     580           0 : }
     581             : 
     582           0 : void MosaicFT::reset(){
     583             :   //call the base class reset
     584           0 :   FTMachine::reset();
     585           0 :   doneWeightImage_p=false;
     586           0 :   convWeightImage_p=nullptr;
     587           0 :   pbConvFunc_p->reset();
     588           0 : }
     589             : 
     590           0 : void MosaicFT::finalizeToSky()
     591             : {
     592           0 :   logIO() << LogOrigin("MosaicFT", "finalizeToSky")  << LogIO::NORMAL;
     593           0 :   logIO() << LogIO::NORMAL2 << "time to massage data " << timemass_p << LogIO::POST;
     594           0 :   logIO() << LogIO::NORMAL2<< "time gridding " << timegrid_p << LogIO::POST;
     595           0 :    timemass_p=0.0;
     596           0 :    timegrid_p=0.0;
     597             :   // Now we flush the cache and report statistics
     598             :   // For memory based, we don't write anything out yet.
     599             :   /*if(isTiled) {
     600             :     logIO() << LogOrigin("MosaicFT", "finalizeToSky")  << LogIO::NORMAL;
     601             :     
     602             :     AlwaysAssert(image, AipsError);
     603             :     AlwaysAssert(imageCache, AipsError);
     604             :     imageCache->flush();
     605             :     ostringstream o;
     606             :     imageCache->showCacheStatistics(o);
     607             :     logIO() << o.str() << LogIO::POST;
     608             :   }
     609             :   */
     610             :   
     611             :   
     612           0 :   if(!doneWeightImage_p){
     613           0 :     if(useDoubleGrid_p){
     614           0 :       convertArray(griddedWeight, griddedWeight2);
     615             :       //Don't need the double-prec grid anymore...
     616           0 :       griddedWeight2.resize();
     617             :     }
     618           0 :     ft_p.c2cFFT(*weightLattice, false);
     619             :     //Get the stokes right
     620           0 :     CoordinateSystem coords=convWeightImage_p->coordinates();
     621           0 :     Int stokesIndex=coords.findCoordinate(Coordinate::STOKES);
     622           0 :     Int npol=1;
     623           0 :     Vector<Int> whichStokes(npol);
     624           0 :     if(stokes_p=="I" || stokes_p=="RR" || stokes_p=="LL" ||stokes_p=="XX" 
     625           0 :        || stokes_p=="YY" || stokes_p=="Q" || stokes_p=="U" || stokes_p=="V"){
     626           0 :       npol=1;
     627           0 :       whichStokes(0)=Stokes::type(stokes_p);
     628             :        // if single plane Q U or V are used...the weight should be the I weight
     629             :       //if(stokes_p=="Q" || stokes_p=="U" || stokes_p=="V")
     630             :       //whichStokes(0)=Stokes::type("I");
     631             :     }
     632           0 :     else if(stokes_p=="IV"){
     633           0 :       npol=2;
     634           0 :       whichStokes.resize(2);
     635           0 :       whichStokes(0)=Stokes::I;
     636           0 :       whichStokes(1)=Stokes::V;
     637             :     }
     638           0 :     else if(stokes_p=="QU"){
     639           0 :       npol=2;
     640           0 :       whichStokes.resize(2);
     641           0 :       whichStokes(0)=Stokes::Q;
     642           0 :       whichStokes(1)=Stokes::U;
     643             :     }
     644           0 :     else if(stokes_p=="RRLL"){
     645           0 :       npol=2;
     646           0 :       whichStokes.resize(2);
     647           0 :       whichStokes(0)=Stokes::RR;
     648           0 :       whichStokes(1)=Stokes::LL;
     649             :     }   
     650           0 :     else if(stokes_p=="XXYY"){
     651           0 :       npol=2;
     652           0 :       whichStokes.resize(2);
     653           0 :       whichStokes(0)=Stokes::XX;
     654           0 :       whichStokes(1)=Stokes::YY;
     655             :     }  
     656           0 :     else if(stokes_p=="IQU"){
     657           0 :       npol=3;
     658           0 :       whichStokes.resize(3);
     659           0 :       whichStokes(0)=Stokes::I;
     660           0 :       whichStokes(1)=Stokes::Q;
     661           0 :       whichStokes(2)=Stokes::U;
     662             :     }
     663           0 :     else if(stokes_p=="IQUV"){
     664           0 :       npol=4;
     665           0 :       whichStokes.resize(4);
     666           0 :       whichStokes(0)=Stokes::I;
     667           0 :       whichStokes(1)=Stokes::Q;
     668           0 :       whichStokes(2)=Stokes::U;
     669           0 :       whichStokes(3)=Stokes::V;
     670             :     } 
     671             :     
     672           0 :     StokesCoordinate newStokesCoord(whichStokes);
     673           0 :     coords.replaceCoordinate(newStokesCoord, stokesIndex);
     674           0 :     IPosition imshp=convWeightImage_p->shape();
     675           0 :     imshp(2)=npol;
     676             : 
     677             : 
     678           0 :     skyCoverage_p=new TempImage<Float> (imshp, coords,1.0);
     679           0 :     IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
     680           0 :                     (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
     681           0 :     IPosition stride(4, 1);
     682           0 :     IPosition trc(blc+image->shape()-stride);
     683             :     
     684             :     // Do the copy
     685           0 :     IPosition start(4, 0);
     686           0 :     convWeightImage_p->put(griddedWeight(blc, trc));
     687           0 :     StokesImageUtil::ToStokesPSF(*skyCoverage_p, *convWeightImage_p);
     688           0 :     if(npol>1){
     689             :       // only the I get it right Q and U or V may end up with zero depending 
     690             :       // if RR or XX
     691           0 :       blc(0)=0; blc(1)=0; blc(3)=0;blc(2)=0;
     692           0 :       trc=skyCoverage_p->shape()-stride;
     693           0 :       trc(2)=0;
     694           0 :       SubImage<Float> isubim(*skyCoverage_p, Slicer(blc, trc, Slicer::endIsLast));
     695           0 :       for (Int k=1; k < npol; ++k){
     696           0 :         blc(2)=k; trc(2)=k;
     697           0 :         SubImage<Float> quvsubim(*skyCoverage_p, Slicer(blc, trc, Slicer::endIsLast), true);
     698           0 :         quvsubim.copyData(isubim);
     699           0 :       }
     700             : 
     701           0 :     }
     702             :     //Store this image in the pbconvfunc object as
     703             :     //it can be used for rescaling or shared by other ftmachines that use
     704             :     //this pbconvfunc
     705           0 :     pbConvFunc_p->setWeightImage(skyCoverage_p);
     706           0 :     if(convWeightImage_p) delete convWeightImage_p;
     707           0 :     convWeightImage_p=0;
     708           0 :     doneWeightImage_p=true;
     709             : 
     710             :     /*
     711             :     if(0){
     712             :       PagedImage<Float> thisScreen(skyCoverage_p->shape(), 
     713             :                                    skyCoverage_p->coordinates(), "Screen");
     714             :       thisScreen.copyData(*skyCoverage_p);
     715             :     }
     716             :     */
     717             : 
     718           0 :   }
     719             : 
     720           0 :   if(!weightLattice.null()) weightLattice=0;
     721           0 :   griddedWeight.resize();
     722           0 :   if(pointingToImage) delete pointingToImage; pointingToImage=0;
     723           0 : }
     724             : 
     725           0 : void MosaicFT::setWeightImage(CountedPtr<ImageInterface<Float> >& wgtimage){
     726           0 :   IPosition shp=wgtimage->shape();
     727           0 :   CoordinateSystem cs=wgtimage->coordinates();
     728           0 :   CountedPtr<TempImage<Float> > wgtim=new TempImage<Float>(shp, cs);
     729           0 :   wgtim->copyData(*(wgtimage));
     730           0 :   skyCoverage_p=wgtim;
     731           0 :   Record rec=skyCoverage_p->miscInfo();
     732             :   //For mosaicFTNew it has the nx*ny factor already in
     733           0 :   rec.define("isscaled", True);
     734           0 :   skyCoverage_p->setMiscInfo(rec);
     735             :   //cerr << "IN SET " << max(wgtimage->get()) << endl;
     736           0 :   pbConvFunc_p->setWeightImage(skyCoverage_p);
     737           0 :   doneWeightImage_p=true;
     738           0 : }
     739             : 
     740           0 : Array<Complex>* MosaicFT::getDataPointer(const IPosition& centerLoc2D,
     741             :                                          Bool readonly) {
     742             :   Array<Complex>* result;
     743             :   // Is tiled: get tiles and set up offsets
     744           0 :   centerLoc(0)=centerLoc2D(0);
     745           0 :   centerLoc(1)=centerLoc2D(1);
     746           0 :   result=&imageCache->tile(offsetLoc,centerLoc, readonly);
     747           0 :   gridder->setOffset(IPosition(2, offsetLoc(0), offsetLoc(1)));
     748           0 :   return result;
     749             : }
     750             : 
     751             : #define NEED_UNDERSCORES
     752             : #if defined(NEED_UNDERSCORES)
     753             : #define sectgmoss2 sectgmoss2_
     754             : #define gmoss2 gmoss2_
     755             : #define sectgmosd2 sectgmosd2_
     756             : #define gmosd2 gmosd2_
     757             : #define sectdmos2 sectdmos2_
     758             : #define dmos2 dmos2_
     759             : #define gmoswgtd gmoswgtd_
     760             : #define gmoswgts gmoswgts_
     761             : #define locuvw locuvw_
     762             : #endif
     763             : 
     764             : extern "C" { 
     765             :   void locuvw(const Double*, const Double*, const Double*, const Int*, const Double*, const Double*, const Int*, 
     766             :               Int*, Int*, Complex*, const Int*, const Int*, const Double*);
     767             :   void gmoswgtd(const Int*/*nvispol*/, const Int*/*nvischan*/,
     768             :                 const Int*/*flag*/, const Int*/*rflag*/, const Float*/*weight*/, const Int*/*nrow*/, 
     769             :                 const Int*/*nx*/, const Int*/*ny*/, const Int*/*npol*/, const Int*/*nchan*/, 
     770             :                 const Int*/*support*/, const Int*/*convsize*/, const Int*/*sampling*/, 
     771             :                 const Int*/*chanmap*/, const Int*/*polmap*/,
     772             :                 DComplex* /*weightgrid*/, Double* /*sumwt*/, const Complex*/*convweight*/, const Int*/*convplanemap*/, 
     773             :                 const Int*/*convchanmap*/,  const Int*/*convpolmap*/, 
     774             :                 const Int*/*nconvplane*/, const Int*/*nconvchan*/, const Int*/*nconvpol*/, const Int*/*rbeg*/, 
     775             :                 const Int*/*rend*/, const Int*/*loc*/, const Int*/*off*/, const Complex*/*phasor*/);
     776             :  void gmoswgts(const Int*/*nvispol*/, const Int*/*nvischan*/,
     777             :                 const Int*/*flag*/, const Int*/*rflag*/, const Float*/*weight*/, const Int*/*nrow*/, 
     778             :                 const Int*/*nx*/, const Int*/*ny*/, const Int*/*npol*/, const Int*/*nchan*/, 
     779             :                 const Int*/*support*/, const Int*/*convsize*/, const Int*/*sampling*/, 
     780             :                 const Int*/*chanmap*/, const Int*/*polmap*/,
     781             :                Complex* /*weightgrid*/, Double*/*sumwt*/, const Complex*/*convweight*/, const Int*/*convplanemap*/, 
     782             :                 const Int*/*convchanmap*/,  const Int*/*convpolmap*/, 
     783             :                 const Int*/*nconvplane*/, const Int*/*nconvchan*/, const Int*/*nconvpol*/, const Int*/*rbeg*/, 
     784             :                 const Int*/*rend*/, const Int*/*loc*/, const Int*/*off*/, const Complex*/*phasor*/);
     785             :   void sectgmosd2(const Complex* /*values*/,
     786             :                   Int* /*nvispol*/, Int* /*nvischan*/,
     787             :                   Int* /*dopsf*/, const Int* /*flag*/, const Int* /*rflag*/, const Float* /*weight*/,
     788             :                   Int* /* nrow*/, DComplex* /*grid*/, Int* /*nx*/, Int* /*ny*/, Int * /*npol*/, Int * /*nchan  */,
     789             :                   Int*/*support*/, Int*/*convsize*/, Int*/*sampling*/, const Complex*/*convfunc*/,
     790             :                   const Int*/*chanmap*/, const Int*/*polmap*/,
     791             :                   Double*/*sumwgt*/, const Int*/*convplanemap*/,
     792             :                   const Int*/*convchanmap*/, const Int*/*convpolmap*/, 
     793             :                   Int*/*nconvplane*/, Int*/*nconvchan*/, Int* /*nconvpol*/,
     794             :                   const Int*/*x0*/,const Int*/*y0*/, const Int*/*nxsub*/, const Int*/*nysub*/, const Int*/*rbeg*/, 
     795             :                   const Int* /*rend*/, const Int*/*loc*/, const Int* /*off*/, const Complex*/*phasor*/);     
     796             : 
     797             :  void sectgmoss2(const Complex* /*values*/,
     798             :                   Int* /*nvispol*/, Int* /*nvischan*/,
     799             :                   Int* /*dopsf*/, const Int* /*flag*/, const Int* /*rflag*/, const Float* /*weight*/,
     800             :                   Int* /* nrow*/, Complex* /*grid*/, Int* /*nx*/, Int* /*ny*/, Int * /*npol*/, Int * /*nchan  */,
     801             :                   Int*/*support*/, Int*/*convsize*/, Int*/*sampling*/, const Complex*/*convfunc*/,
     802             :                   const Int*/*chanmap*/, const Int*/*polmap*/,
     803             :                   Double*/*sumwgt*/, const Int*/*convplanemap*/,
     804             :                   const Int*/*convchanmap*/, const Int*/*convpolmap*/, 
     805             :                   Int*/*nconvplane*/, Int*/*nconvchan*/, Int* /*nconvpol*/,
     806             :                   const Int*/*x0*/,const Int*/*y0*/, const Int*/*nxsub*/, const Int*/*nysub*/, const Int*/*rbeg*/, 
     807             :                   const Int* /*rend*/, const Int*/*loc*/, const Int* /*off*/, const Complex*/*phasor*/);     
     808             : 
     809             : 
     810             :   void gmosd2(const Double*,
     811             :               Double*,
     812             :               const Complex*,
     813             :               Int*,
     814             :               Int*,
     815             :               Int*,
     816             :               const Int*,
     817             :               const Int*,
     818             :               const Float*,
     819             :               Int*,
     820             :               Int*,
     821             :               Double*,
     822             :               Double*,
     823             :               DComplex*,
     824             :               Int*,
     825             :               Int*,
     826             :               Int *,
     827             :               Int *,
     828             :               const Double*,
     829             :               const Double*,
     830             :               Int*,
     831             :               Int*,
     832             :               Int*,
     833             :               const Complex*,
     834             :               Int*,
     835             :               Int*,
     836             :               Double*,
     837             :               DComplex*,
     838             :               Complex*,
     839             :               Int*,
     840             :               Int*,
     841             :               Int*,Int*, Int*, Int*, Int*);
     842             :   /*  void gmoss(const Double*,
     843             :               Double*,
     844             :               const Complex*,
     845             :               Int*,
     846             :               Int*,
     847             :               Int*,
     848             :               const Int*,
     849             :               const Int*,
     850             :               const Float*,
     851             :               Int*,
     852             :               Int*,
     853             :               Double*,
     854             :               Double*,
     855             :               Complex*,
     856             :               Int*,
     857             :               Int*,
     858             :               Int *,
     859             :               Int *,
     860             :               const Double*,
     861             :               const Double*,
     862             :               Int*,
     863             :               Int*,
     864             :               Int*,
     865             :               const Complex*,
     866             :               Int*,
     867             :               Int*,
     868             :               Double*,
     869             :               Complex*,
     870             :               Complex*,
     871             :               Int*,
     872             :               Int*,
     873             :               Int*);
     874             :   */
     875             :     void gmoss2(const Double*,
     876             :               Double*,
     877             :               const Complex*,
     878             :               Int*,
     879             :               Int*,
     880             :               Int*,
     881             :               const Int*,
     882             :               const Int*,
     883             :               const Float*,
     884             :               Int*, //10
     885             :               Int*,
     886             :               Double*,
     887             :               Double*,
     888             :               Complex*,
     889             :               Int*,
     890             :               Int*,
     891             :               Int *,
     892             :               Int *,
     893             :               const Double*,
     894             :               const Double*, //20
     895             :               Int*,
     896             :               Int*,
     897             :               Int*,
     898             :               const Complex*,
     899             :               Int*,
     900             :               Int*,
     901             :               Double*,
     902             :               Complex*,
     903             :               Complex*,
     904             :               Int*, //30
     905             :               Int*,
     906             :               Int*, Int*, Int*, Int*, Int*);
     907             : 
     908             :   /* void dmos(const Double*,
     909             :               Double*,
     910             :               Complex*,
     911             :               Int*,
     912             :               Int*,
     913             :               const Int*,
     914             :               const Int*,
     915             :               Int*,
     916             :               Int*,
     917             :             Double*, //10
     918             :               Double*,
     919             :               const Complex*,
     920             :               Int*,
     921             :               Int*,
     922             :               Int *,
     923             :               Int *,
     924             :               const Double*,
     925             :               const Double*,
     926             :               Int*,
     927             :             Int*,//20
     928             :               Int*,
     929             :               const Complex*,
     930             :               Int*,
     931             :               Int*,
     932             :               Int*,
     933             :               Int*);
     934             :   */
     935             : 
     936             :   void dmos2(const Double*,
     937             :               Double*,
     938             :               Complex*,
     939             :               Int*,
     940             :               Int*,
     941             :               const Int*,
     942             :               const Int*,
     943             :               Int*,
     944             :               Int*,
     945             :               Double*,
     946             :               Double*,
     947             :               const Complex*,
     948             :               Int*,
     949             :               Int*,
     950             :               Int *,
     951             :               Int *,
     952             :               const Double*,
     953             :               const Double*,
     954             :               Int*,
     955             :               Int*,
     956             :               Int*,
     957             :               const Complex*,
     958             :               Int*,
     959             :               Int*,
     960             :               Int*,
     961             :               Int*, Int*, Int*, Int*, Int*);
     962             :   void sectdmos2(Complex*,
     963             :               Int*,
     964             :               Int*,
     965             :               const Int*,
     966             :               const Int*,
     967             :                  Int*,
     968             :               const Complex*,
     969             :               Int*,
     970             :               Int*,
     971             :               Int *,
     972             :                  Int *,
     973             :               Int*,
     974             :               Int*,
     975             :               Int*,
     976             :               const Complex*,
     977             :               const Int*,
     978             :               const Int*,
     979             :               const Int*,
     980             :               const  Int*, 
     981             :               const Int*, 
     982             :               Int*, Int*, Int*,
     983             :                  //rbeg
     984             :                  const Int*,
     985             :                  const Int*,
     986             :                  const Int*,
     987             :                  const Int*,
     988             :                  const Complex*);
     989             : 
     990             :              
     991             : 
     992             : }
     993           0 : void MosaicFT::put(const vi::VisBuffer2& vb, Int row, Bool dopsf,
     994             :                    FTMachine::Type type)
     995             : {
     996             : 
     997             : 
     998             :   
     999             :   
    1000           0 :   Timer tim;
    1001           0 :   tim.mark();
    1002             :  
    1003           0 :   matchChannel(vb);
    1004             :  
    1005             : 
    1006             :   //cerr << "CHANMAP " << chanMap << endl;
    1007             :   //No point in reading data if its not matching in frequency
    1008           0 :   if(max(chanMap)==-1)
    1009           0 :     return;
    1010             : 
    1011             :   //const Matrix<Float> *imagingweight;
    1012             :   //imagingweight=&(vb.imagingWeight());
    1013           0 :   Matrix<Float> imagingweight;
    1014           0 :   getImagingWeight(imagingweight, vb);
    1015             : 
    1016           0 :   if(dopsf) type=FTMachine::PSF;
    1017             : 
    1018           0 :   Cube<Complex> data;
    1019             :   //Fortran gridder need the flag as ints 
    1020           0 :   Cube<Int> flags;
    1021           0 :   Matrix<Float> elWeight;
    1022           0 :   interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, type);
    1023             :   
    1024             :  
    1025             : 
    1026             :   Bool iswgtCopy;
    1027             :   const Float *wgtStorage;
    1028           0 :   wgtStorage=elWeight.getStorage(iswgtCopy);
    1029             : 
    1030             : 
    1031             :   
    1032             : 
    1033             :   Bool isCopy;
    1034           0 :   const Complex *datStorage=0;
    1035             : 
    1036             :   // cerr << "dopsf " << dopsf << " isWeightCopy " << iswgtCopy << "  " << wgtStorage<< endl;
    1037           0 :   if(!dopsf)
    1038           0 :     datStorage=data.getStorage(isCopy);
    1039             :     
    1040             :   
    1041             :   // If row is -1 then we pass through all rows
    1042             :   Int startRow, endRow, nRow;
    1043           0 :   if (row==-1) {
    1044           0 :     nRow=vb.nRows();
    1045           0 :     startRow=0;
    1046           0 :     endRow=nRow-1;
    1047             :   } else {
    1048           0 :     nRow=1;
    1049           0 :     startRow=row;
    1050           0 :     endRow=row;
    1051             :   }
    1052             :   
    1053             :   // Get the uvws in a form that Fortran can use and do that
    1054             :   // necessary phase rotation. 
    1055           0 :   Matrix<Double> uvw(negateUV(vb));
    1056           0 :   Vector<Double> dphase(vb.nRows());
    1057           0 :   dphase=0.0;
    1058             :  
    1059           0 :   doUVWRotation_p=true;
    1060           0 :   girarUVW(uvw, dphase, vb);
    1061           0 :   refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
    1062             :   // This needs to be after the interp to get the interpolated channels
    1063             :   //Also has to be after rotateuvw in case tracking is on
    1064             :   //cerr << "orig " << vb.uvw().row(2) << endl;
    1065             :   //vi::VisBuffer2& vbRotuvw=const_cast<vi::VisBuffer2&>(vb);
    1066             :   //vbRotuvw.setUvw(uvw);
    1067             : 
    1068             :   // cerr << "rot  " << vbRotuvw.uvw().row(2) << endl;
    1069             : 
    1070           0 :   findConvFunction(*image, vb, uvw);
    1071             :   
    1072             :   //cerr << "Put convsup " << convSupport << " max min convFunc " << max(convFunc) << "   " << min(convFunc) << "  "  << max(weightConvFunc_p) << min(weightConvFunc_p)  << "SHP " << convFunc.shape() << "   " << weightConvFunc_p.shape() << endl;
    1073             :   //cerr << "convRowMap " << convRowMap_p  << " " << convChanMap_p << "  " << convPolMap_p << endl; 
    1074             :   //nothing to grid here as the pointing resulted in a zero support convfunc
    1075           0 :   if(convSupport <= 0)
    1076           0 :     return;
    1077             :   
    1078             :   // Get the pointing positions. This can easily consume a lot 
    1079             :   // of time thus we are for now assuming a field per 
    1080             :   // vb chunk...need to change that accordingly if we start using
    1081             :   // multiple pointings per vb.
    1082             :   //Warning 
    1083             : 
    1084             :   // Take care of translation of Bools to Integer
    1085           0 :   Int idopsf=0;
    1086           0 :   if(dopsf) idopsf=1;
    1087             :   
    1088             :   
    1089           0 :   Vector<Int> rowFlags(vb.nRows());
    1090           0 :   rowFlags=0;
    1091           0 :   rowFlags(vb.flagRow())=true;
    1092           0 :   if(!usezero_p) {
    1093           0 :     for (Int rownr=startRow; rownr<=endRow; rownr++) {
    1094           0 :       if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
    1095             :     }
    1096             :   }
    1097             :   
    1098             :   
    1099             : 
    1100             :   //cerr << "convSamp " << convSampling << " convsupp " << convSupport << " consize " << convSize << " convFunc " << convFunc.shape() << endl;
    1101             :   //TESTOO
    1102             :   /*{
    1103             :     ArrayIterator<Complex> itC(convFunc, IPosition(2,0,1));
    1104             :     ArrayIterator<Complex> itW(weightConvFunc_p, IPosition(2,0,1));
    1105             :     itC.origin();
    1106             :     itW.origin();
    1107             :     Int k=0;
    1108             :     while(!itC.pastEnd()){
    1109             :       cerr << k << "sum conv plane " << sum(itC.array()) << "  wt " << sum(itW.array()) << endl;
    1110             : 
    1111             :       itC.next();
    1112             :       itW.next();
    1113             :       ++k;
    1114             :     }
    1115             : 
    1116             :     }*/
    1117             :   //TESTOO
    1118             :   
    1119             :   //Tell the gridder to grid the weights too ...need to do that once only
    1120             :   //Int doWeightGridding=1;
    1121             :   //if(doneWeightImage_p)
    1122             :   //  doWeightGridding=-1;
    1123             :   Bool del;
    1124             :   //    IPosition s(flags.shape());
    1125           0 :   const IPosition& fs=flags.shape();
    1126             :   //cerr << "flags shape " << fs << endl;
    1127           0 :   std::vector<Int>s(fs.begin(), fs.end());
    1128           0 :   Int nvp=s[0];
    1129           0 :   Int nvc=s[1];
    1130           0 :   Int nvisrow=s[2];
    1131           0 :   Int csamp=convSampling;
    1132             :   Bool uvwcopy; 
    1133           0 :   const Double *uvwstor=uvw.getStorage(uvwcopy);
    1134             :   Bool gridcopy;
    1135             :   Bool convcopy;
    1136             :   Bool wconvcopy;
    1137           0 :   const Complex *convstor=convFunc.getStorage(convcopy);
    1138           0 :   const Complex *wconvstor=weightConvFunc_p.getStorage(wconvcopy);
    1139           0 :   Int nPolConv=convFunc.shape()[2];
    1140           0 :   Int nChanConv=convFunc.shape()[3];
    1141           0 :   Int nConvFunc=convFunc.shape()(4);
    1142             :   Bool weightcopy;
    1143             :   ////////**************************
    1144           0 :   Cube<Int> loc(2, nvc, nRow);
    1145           0 :   Cube<Int> off(2, nvc, nRow);
    1146           0 :   Matrix<Complex> phasor(nvc, nRow);
    1147             :   Bool delphase;
    1148           0 :   Complex * phasorstor=phasor.getStorage(delphase);
    1149           0 :   const Double * visfreqstor=interpVisFreq_p.getStorage(del);
    1150           0 :   const Double * scalestor=uvScale.getStorage(del);
    1151           0 :   const Double * offsetstor=uvOffset.getStorage(del);
    1152           0 :   Int * locstor=loc.getStorage(del);
    1153           0 :   Int * offstor=off.getStorage(del);
    1154           0 :   const Double *dpstor=dphase.getStorage(del);
    1155             :   Int irow;
    1156           0 :   Int nth=1;
    1157             : #ifdef _OPENMP
    1158           0 :   if(numthreads_p >0){
    1159           0 :     nth=min(numthreads_p, omp_get_max_threads());
    1160             :   }
    1161             :   else{   
    1162           0 :     nth= omp_get_max_threads();
    1163             :   }
    1164             :   //nth=min(4,nth);
    1165             : #endif
    1166           0 :   Double cinv=Double(1.0)/C::c;
    1167             :  
    1168           0 :   Int dow=0;
    1169           0 : #pragma omp parallel default(none) private(irow) firstprivate(visfreqstor, nvc, scalestor, offsetstor, csamp, phasorstor, uvwstor, locstor, offstor, dpstor, dow, cinv) shared(startRow, endRow) num_threads(nth)  
    1170             : {
    1171             : #pragma omp for
    1172             :   for (irow=startRow; irow<=endRow;irow++){
    1173             :     /*locateuvw(uvwstor,dpstor, visfreqstor, nvc, scalestor, offsetstor, csamp, 
    1174             :               locstor, 
    1175             :               offstor, phasorstor, irow, false);*/
    1176             :     locuvw(uvwstor, dpstor, visfreqstor, &nvc, scalestor, offsetstor, &csamp, locstor, offstor, phasorstor, &irow, &dow, &cinv);
    1177             :   }  
    1178             : 
    1179             :  }//end pragma parallel
    1180             : 
    1181             : 
    1182             :  
    1183           0 :  timemass_p +=tim.real();
    1184             :  Int  ixsub, iysub, icounter;
    1185           0 :  ixsub=1;
    1186           0 :  iysub=1;
    1187             :   //////***********************DEBUGGING
    1188             :   //nth=1;
    1189             :   ////////***************
    1190           0 :   if (nth >3){
    1191           0 :     ixsub=8;
    1192           0 :     iysub=8; 
    1193             :   }
    1194           0 :   else if(nth >1){
    1195           0 :      ixsub=2;
    1196           0 :      iysub=2; 
    1197             :   }
    1198           0 :   Int rbeg=startRow+1;
    1199           0 :   Int rend=endRow+1;
    1200           0 :   Block<Matrix<Double> > sumwgt(ixsub*iysub);
    1201           0 :   Vector<Double *> swgtptr(ixsub*iysub);
    1202           0 :   Vector<Bool> swgtdel(ixsub*iysub);
    1203           0 :   for (icounter=0; icounter < ixsub*iysub; ++icounter){
    1204           0 :     sumwgt[icounter].resize(sumWeight.shape());
    1205           0 :     sumwgt[icounter].set(0.0);
    1206           0 :     swgtptr[icounter]=sumwgt[icounter].getStorage(swgtdel(icounter));
    1207             :   }
    1208             :   //cerr << "done thread " << doneThreadPartition_p << "  " << ixsub*iysub << endl;
    1209           0 :    if(doneThreadPartition_p < 0){
    1210           0 :     xsect_p.resize(ixsub*iysub);
    1211           0 :     ysect_p.resize(ixsub*iysub);
    1212           0 :     nxsect_p.resize(ixsub*iysub);
    1213           0 :     nysect_p.resize(ixsub*iysub);
    1214           0 :     for (icounter=0; icounter < ixsub*iysub; ++icounter){
    1215           0 :       findGridSector(nx, ny, ixsub, iysub, 0, 0, icounter, xsect_p(icounter), ysect_p(icounter), nxsect_p(icounter), nysect_p(icounter), true);
    1216             :     }
    1217             :   }
    1218           0 :    Vector<Int> xsect, ysect, nxsect, nysect;
    1219           0 :    xsect=xsect_p; ysect=ysect_p; nxsect=nxsect_p; nysect=nysect_p;
    1220             :    //cerr << xsect.shape() << "  " << xsect << endl;
    1221           0 :   const Int* pmapstor=polMap.getStorage(del);
    1222           0 :   const Int* cmapstor=chanMap.getStorage(del);
    1223             : // Dummy sumwt for gridweight part
    1224           0 :   Matrix<Double> dumSumWeight(npol, nchan);
    1225           0 :   dumSumWeight=sumWeight;
    1226             :   Bool isDSWC;
    1227           0 :   Double *dsumwtstor=dumSumWeight.getStorage(isDSWC);
    1228           0 :   Int nc=nchan;
    1229           0 :   Int np=npol;
    1230           0 :   Int nxp=nx;
    1231           0 :   Int nyp=ny;
    1232           0 :   Int csize=convSize;
    1233           0 :   const Int * flagstor=flags.getStorage(del);
    1234           0 :   const Int * rowflagstor=rowFlags.getStorage(del);
    1235           0 :   Int csupp=convSupport;
    1236           0 :   const Int *convrowmapstor=convRowMap_p.getStorage(del);
    1237           0 :   const Int *convchanmapstor=convChanMap_p.getStorage(del);
    1238           0 :   const Int *convpolmapstor=convPolMap_p.getStorage(del);
    1239             :   ///
    1240             : 
    1241             :   
    1242             :   ////////***************************
    1243           0 :   tim.mark(); 
    1244             : 
    1245           0 :   if(useDoubleGrid_p) {
    1246           0 :     DComplex *gridstor=griddedData2.getStorage(gridcopy);
    1247             :     
    1248           0 : #pragma omp parallel default(none) private(icounter, del) firstprivate(idopsf, /*doWeightGridding,*/ datStorage, wgtStorage, flagstor, rowflagstor, convstor, wconvstor, pmapstor, cmapstor, gridstor,  csupp, nxp, nyp, np, nc,ixsub, iysub, rend, rbeg, csamp, csize, nvp, nvc, nvisrow, phasorstor, locstor, offstor, convrowmapstor, convchanmapstor, convpolmapstor, nPolConv, nChanConv, nConvFunc,xsect, ysect, nxsect, nysect) shared(swgtptr) 
    1249             :     {   
    1250             : #pragma omp for schedule(dynamic)      
    1251             :     for(icounter=0; icounter < ixsub*iysub; ++icounter){
    1252             :       Int x0=xsect(icounter);
    1253             :       Int y0=ysect(icounter);
    1254             :       Int nxsub=nxsect(icounter);
    1255             :       Int nysub=nysect(icounter);
    1256             :       
    1257             :       /*
    1258             :       ix= (icounter+1)-((icounter)/ixsub)*ixsub;
    1259             :       iy=(icounter)/ixsub+1;
    1260             :       y0=(nyp/iysub)*(iy-1)+1;
    1261             :       nysub=nyp/iysub;
    1262             :       if( iy == iysub) {
    1263             :         nysub=nyp-(nyp/iysub)*(iy-1);
    1264             :       }
    1265             :       x0=(nxp/ixsub)*(ix-1)+1;
    1266             :       nxsub=nxp/ixsub;
    1267             :       if( ix == ixsub){
    1268             :         nxsub=nxp-(nxp/ixsub)*(ix-1);
    1269             :       } 
    1270             :       */
    1271             : 
    1272             :     sectgmosd2(datStorage,
    1273             :            &nvp,
    1274             :            &nvc,
    1275             :            &idopsf,
    1276             :            flagstor,
    1277             :            rowflagstor,
    1278             :            wgtStorage,
    1279             :            &nvisrow,
    1280             :            gridstor,
    1281             :            &nxp,
    1282             :            &nyp,
    1283             :            &np,
    1284             :            &nc,
    1285             :            &csupp, 
    1286             :            &csize,
    1287             :            &csamp,
    1288             :            convstor,
    1289             :            cmapstor,
    1290             :            pmapstor,
    1291             :            swgtptr[icounter],
    1292             :            convrowmapstor,
    1293             :            convchanmapstor,
    1294             :            convpolmapstor,
    1295             :                &nConvFunc, &nChanConv, &nPolConv,
    1296             :                &x0, &y0, &nxsub, &nysub, &rbeg, &rend, locstor, offstor,
    1297             :                  phasorstor
    1298             :                );
    1299             :     }
    1300             :     }//end pragma parallel
    1301           0 :     for (icounter=0; icounter < ixsub*iysub; ++icounter){
    1302           0 :       sumwgt[icounter].putStorage(swgtptr[icounter],swgtdel[icounter]);
    1303           0 :       sumWeight=sumWeight+sumwgt[icounter];
    1304             :     }    
    1305             : 
    1306             :     //cerr << "SUMWEIG " << sumWeight << endl;
    1307           0 :     griddedData2.putStorage(gridstor, gridcopy);
    1308           0 :     if(dopsf && (nth >4))
    1309           0 :       tweakGridSector(nx, ny, ixsub, iysub);
    1310           0 :     timegrid_p+=tim.real();
    1311           0 :     tim.mark();
    1312           0 :     if(!doneWeightImage_p){
    1313             :       //This can be parallelized by making copy of the central part of the griddedWeight
    1314             :       //and adding it after dooing the gridding
    1315           0 :       DComplex *gridwgtstor=griddedWeight2.getStorage(weightcopy);
    1316           0 :       gmoswgtd(&nvp, &nvc,flagstor, rowflagstor, wgtStorage, &nvisrow, 
    1317             :                &nxp, &nyp, &np, &nc, &csupp, &csize, &csamp, 
    1318             :                cmapstor, pmapstor,
    1319             :                gridwgtstor, dsumwtstor, wconvstor, convrowmapstor, 
    1320             :                convchanmapstor,  convpolmapstor, 
    1321             :                &nConvFunc, &nChanConv, &nPolConv, &rbeg, 
    1322             :                &rend, locstor, offstor, phasorstor);
    1323           0 :       griddedWeight2.putStorage(gridwgtstor, weightcopy);
    1324             :     
    1325             :     }
    1326           0 :     timemass_p+=tim.real();
    1327             :   }
    1328             :   else {
    1329             :     //cerr << "maps "  << convChanMap_p << "   " << chanMap  << endl;
    1330             :     //cerr << "nchan " << nchan << "  nchanconv " << nChanConv << endl;
    1331           0 :     Complex *gridstor=griddedData.getStorage(gridcopy);
    1332           0 : #pragma omp parallel default(none) private(icounter, del) firstprivate(idopsf, /*doWeightGridding,*/ datStorage, wgtStorage, flagstor, rowflagstor, convstor, wconvstor, pmapstor, cmapstor, gridstor, csupp, nxp, nyp, np, nc,ixsub, iysub, rend, rbeg, csamp, csize, nvp, nvc, nvisrow, phasorstor, locstor, offstor, convrowmapstor, convchanmapstor, convpolmapstor, nPolConv, nChanConv, nConvFunc, xsect, ysect, nxsect, nysect)  shared(swgtptr) 
    1333             :     {   
    1334             : #pragma omp for schedule(dynamic)      
    1335             :       for(icounter=0; icounter < ixsub*iysub; ++icounter){
    1336             :         /*ix= (icounter+1)-((icounter)/ixsub)*ixsub;
    1337             :         iy=(icounter)/ixsub+1;
    1338             :         y0=(nyp/iysub)*(iy-1)+1;
    1339             :         nysub=nyp/iysub;
    1340             :         if( iy == iysub) {
    1341             :           nysub=nyp-(nyp/iysub)*(iy-1);
    1342             :         }
    1343             :         x0=(nxp/ixsub)*(ix-1)+1;
    1344             :         nxsub=nxp/ixsub;
    1345             :         if( ix == ixsub){
    1346             :           nxsub=nxp-(nxp/ixsub)*(ix-1);
    1347             :         } 
    1348             :         */
    1349             :         Int x0=xsect(icounter);
    1350             :         Int y0=ysect(icounter);
    1351             :         Int nxsub=nxsect(icounter);
    1352             :         Int nysub=nysect(icounter);
    1353             : 
    1354             :            sectgmoss2(datStorage,
    1355             :            &nvp,
    1356             :            &nvc,
    1357             :            &idopsf,
    1358             :            flagstor,
    1359             :            rowflagstor,
    1360             :            wgtStorage,
    1361             :            &nvisrow,
    1362             :            gridstor,
    1363             :            &nxp,
    1364             :            &nyp,
    1365             :            &np,
    1366             :            &nc,
    1367             :            &csupp, 
    1368             :            &csize,
    1369             :            &csamp,
    1370             :            convstor,
    1371             :            cmapstor,
    1372             :            pmapstor,
    1373             :            swgtptr[icounter],
    1374             :            convrowmapstor,
    1375             :            convchanmapstor,
    1376             :            convpolmapstor,
    1377             :                &nConvFunc, &nChanConv, &nPolConv,
    1378             :                &x0, &y0, &nxsub, &nysub, &rbeg, &rend, locstor, offstor,
    1379             :                  phasorstor
    1380             :                );
    1381             : 
    1382             : 
    1383             :     }
    1384             :     } //end pragma   
    1385           0 :      for (icounter=0; icounter < ixsub*iysub; ++icounter){
    1386           0 :        sumwgt[icounter].putStorage(swgtptr[icounter],swgtdel[icounter]);
    1387           0 :        sumWeight=sumWeight+sumwgt[icounter];
    1388             :     }
    1389           0 :     griddedData.putStorage(gridstor, gridcopy);
    1390           0 :     if(dopsf && (nth > 4))
    1391           0 :       tweakGridSector(nx, ny, ixsub, iysub);
    1392           0 :     timegrid_p+=tim.real();
    1393           0 :     tim.mark();
    1394           0 :     if(!doneWeightImage_p){
    1395           0 :       Complex *gridwgtstor=griddedWeight.getStorage(weightcopy);
    1396           0 :       gmoswgts(&nvp, &nvc,flagstor, rowflagstor, wgtStorage, &nvisrow, 
    1397             :                &nxp, &nyp, &np, &nc, &csupp, &csize, &csamp, 
    1398             :                cmapstor, pmapstor,
    1399             :                gridwgtstor, dsumwtstor, wconvstor, convrowmapstor, 
    1400             :                convchanmapstor,  convpolmapstor, 
    1401             :                &nConvFunc, &nChanConv, &nPolConv, &rbeg, 
    1402             :                &rend, locstor, offstor, phasorstor);
    1403           0 :       griddedWeight.putStorage(gridwgtstor, weightcopy);
    1404             :     
    1405             :     }
    1406           0 :     timemass_p+=tim.real();
    1407             :   }
    1408           0 :   convFunc.freeStorage(convstor, convcopy);
    1409           0 :   weightConvFunc_p.freeStorage(wconvstor, wconvcopy);
    1410           0 :   dumSumWeight.putStorage(dsumwtstor, isDSWC);
    1411             :   //cerr << "dumSumwe " << dumSumWeight << endl;
    1412           0 :   uvw.freeStorage(uvwstor, uvwcopy);
    1413           0 :   if(!dopsf)
    1414           0 :     data.freeStorage(datStorage, isCopy);
    1415             : 
    1416           0 :   elWeight.freeStorage(wgtStorage,iswgtCopy);
    1417             :   
    1418             : 
    1419             : 
    1420             : 
    1421           0 : }
    1422             : 
    1423           0 : void MosaicFT::gridImgWeights(const vi::VisBuffer2& vb){
    1424             : 
    1425           0 :   if(doneWeightImage_p)
    1426           0 :     return;
    1427           0 :   matchChannel(vb);
    1428             :  
    1429             :   
    1430             :   //cerr << "CHANMAP " << chanMap << endl;
    1431             :   //No point in reading data if its not matching in frequency
    1432           0 :   if(max(chanMap)==-1)
    1433           0 :     return;
    1434             : 
    1435             :   Int startRow, endRow, nRow;
    1436           0 :   nRow=vb.nRows();
    1437           0 :   startRow=0;
    1438           0 :   endRow=nRow-1;
    1439             :   
    1440             :   
    1441             :   //const Matrix<Float> *imagingweight;
    1442             :   //imagingweight=&(vb.imagingWeight());
    1443           0 :   Matrix<Float> imagingweight;
    1444           0 :   getImagingWeight(imagingweight, vb);
    1445             : 
    1446             : 
    1447           0 :   Cube<Complex> data;
    1448             :   //Fortran gridder need the flag as ints 
    1449           0 :   Cube<Int> flags;
    1450           0 :   Matrix<Float> elWeight;
    1451           0 :   interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, FTMachine::PSF);
    1452             :   
    1453             :  
    1454             : 
    1455             :   Bool iswgtCopy;
    1456             :   const Float *wgtStorage;
    1457           0 :   wgtStorage=elWeight.getStorage(iswgtCopy);
    1458             :   Bool issumWgtCopy;
    1459           0 :   Double* sumwgtstor=sumWeight.getStorage(issumWgtCopy);
    1460             : 
    1461             :   
    1462             :  
    1463             :   // Get the uvws in a form that Fortran can use and do that
    1464             :   // necessary phase rotation. 
    1465           0 :   Matrix<Double> uvw(negateUV(vb));
    1466           0 :   Vector<Double> dphase(vb.nRows());
    1467           0 :   dphase=0.0;
    1468             :  
    1469           0 :   doUVWRotation_p=true;
    1470           0 :   girarUVW(uvw, dphase, vb);
    1471           0 :   refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
    1472             :   // This needs to be after the interp to get the interpolated channels
    1473             :   //Also has to be after rotateuvw in case tracking is on
    1474             :   //vi::VisBuffer2& vbRotuvw=const_cast<vi::VisBuffer2&>(vb);
    1475             :   //vbRotuvw.setUvw(uvw);
    1476             : 
    1477           0 :   findConvFunction(*image, vb, uvw);
    1478             :   //nothing to grid here as the pointing resulted in a zero support convfunc
    1479           0 :   if(convSupport <= 0)
    1480           0 :     return;
    1481             :   
    1482             :   Bool del;
    1483             :   
    1484           0 :   const Int* pmapstor=polMap.getStorage(del);
    1485           0 :   const Int* cmapstor=chanMap.getStorage(del);
    1486             :   
    1487           0 :   Vector<Int> rowFlags(vb.nRows());
    1488           0 :   rowFlags=0;
    1489           0 :   rowFlags(vb.flagRow())=true;
    1490           0 :   if(!usezero_p) {
    1491           0 :     for (uInt rownr=0; rownr< vb.nRows(); rownr++) {
    1492           0 :       if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
    1493             :     }
    1494             :   }
    1495             : 
    1496             :   //Fortran indexing
    1497             :   
    1498           0 :   Int rbeg=1;
    1499           0 :   Int rend=vb.nRows();
    1500             : 
    1501           0 :   const Int * flagstor=flags.getStorage(del);
    1502           0 :   const Int * rowflagstor=rowFlags.getStorage(del);
    1503             : 
    1504           0 :   const Int *convrowmapstor=convRowMap_p.getStorage(del);
    1505           0 :   const Int *convchanmapstor=convChanMap_p.getStorage(del);
    1506           0 :   const Int *convpolmapstor=convPolMap_p.getStorage(del);
    1507             :   
    1508             :   //Tell the gridder to grid the weights too ...need to do that once only
    1509             :   //Int doWeightGridding=1;
    1510             :   //if(doneWeightImage_p)
    1511             :   //  doWeightGridding=-1;
    1512             :   //    IPosition s(flags.shape());
    1513           0 :   const IPosition& fs=flags.shape();
    1514             :   //cerr << "flags shape " << fs << endl;
    1515           0 :   std::vector<Int>s(fs.begin(), fs.end());
    1516           0 :   Int nvp=s[0];
    1517           0 :   Int nvc=s[1];
    1518           0 :   Int nvisrow=s[2];
    1519           0 :   Int csamp=convSampling;
    1520             :   Bool uvwcopy; 
    1521           0 :   const Double *uvwstor=uvw.getStorage(uvwcopy);
    1522             :   Bool gridcopy;
    1523             :   Bool convcopy;
    1524             :   Bool wconvcopy;
    1525           0 :   const Complex *wconvstor=weightConvFunc_p.getStorage(wconvcopy);
    1526           0 :   Int nPolConv=convFunc.shape()[2];
    1527           0 :   Int nChanConv=convFunc.shape()[3];
    1528           0 :   Int nConvFunc=convFunc.shape()(4);
    1529             :   Bool weightcopy;
    1530             :   ////////**************************
    1531           0 :   Cube<Int> loc(2, nvc, vb.nRows());
    1532           0 :   Cube<Int> off(2, nvc, vb.nRows());
    1533           0 :   Matrix<Complex> phasor(nvc, vb.nRows());
    1534             :   Bool delphase;
    1535           0 :   Complex * phasorstor=phasor.getStorage(delphase);
    1536           0 :   const Double * visfreqstor=interpVisFreq_p.getStorage(del);
    1537           0 :   const Double * scalestor=uvScale.getStorage(del);
    1538           0 :   const Double * offsetstor=uvOffset.getStorage(del);
    1539           0 :   Int * locstor=loc.getStorage(del);
    1540           0 :   Int * offstor=off.getStorage(del);
    1541           0 :   const Double *dpstor=dphase.getStorage(del);
    1542             : 
    1543             :   Int irow;
    1544           0 :   Int nth=1;
    1545             : #ifdef _OPENMP
    1546           0 :   if(numthreads_p >0){
    1547           0 :     nth=min(numthreads_p, omp_get_max_threads());
    1548             :   }
    1549             :   else{   
    1550           0 :     nth= omp_get_max_threads();
    1551             :   }
    1552             :   //nth=min(4,nth);
    1553             : #endif
    1554             : 
    1555           0 :   Double cinv=Double(1.0)/C::c;
    1556             :  
    1557           0 :   Int dow=0;
    1558             : 
    1559           0 : #pragma omp parallel default(none) private(irow) firstprivate(visfreqstor, nvc, scalestor, offsetstor, csamp, phasorstor, uvwstor, locstor, offstor, dpstor, dow, cinv) shared(startRow, endRow) num_threads(nth)  
    1560             : {
    1561             : #pragma omp for
    1562             :   for (irow=startRow; irow<=endRow;irow++){
    1563             :     /*locateuvw(uvwstor,dpstor, visfreqstor, nvc, scalestor, offsetstor, csamp, 
    1564             :               locstor, 
    1565             :               offstor, phasorstor, irow, false);*/
    1566             :     locuvw(uvwstor, dpstor, visfreqstor, &nvc, scalestor, offsetstor, &csamp, locstor, offstor, phasorstor, &irow, &dow, &cinv);
    1567             :   }  
    1568             : 
    1569             :  }//end pragma parallel
    1570             : 
    1571             : 
    1572             : 
    1573             : 
    1574           0 :   if(useDoubleGrid_p) {
    1575             :       //This can be parallelized by making copy of the central part of the griddedWeight
    1576             :       //and adding it after dooing the gridding
    1577           0 :       DComplex *gridwgtstor=griddedWeight2.getStorage(weightcopy);
    1578           0 :       gmoswgtd(&nvp, &nvc,flagstor, rowflagstor, wgtStorage, &nvisrow, 
    1579           0 :                &nx, &ny, &npol, &nchan, &convSupport, &convSize, &convSampling, 
    1580             :                cmapstor, pmapstor,
    1581             :                gridwgtstor, sumwgtstor, wconvstor, convrowmapstor, 
    1582             :                convchanmapstor,  convpolmapstor, 
    1583             :                &nConvFunc, &nChanConv, &nPolConv, &rbeg, 
    1584             :                &rend, locstor, offstor, phasorstor);
    1585           0 :       griddedWeight2.putStorage(gridwgtstor, weightcopy);
    1586             :     
    1587             :     
    1588             : 
    1589             : 
    1590             : 
    1591             :   }
    1592             :   else{
    1593           0 :     Complex *gridwgtstor=griddedWeight.getStorage(weightcopy);
    1594           0 :     gmoswgts(&nvp, &nvc,flagstor, rowflagstor, wgtStorage, &nvisrow, 
    1595           0 :              &nx, &ny, &npol, &nchan, &convSupport, &convSize, &convSampling, 
    1596             :              cmapstor, pmapstor,
    1597             :              gridwgtstor, sumwgtstor, wconvstor, convrowmapstor, 
    1598             :              convchanmapstor,  convpolmapstor, 
    1599             :              &nConvFunc, &nChanConv, &nPolConv, &rbeg, 
    1600             :              &rend, locstor, offstor, phasorstor);
    1601           0 :     griddedWeight.putStorage(gridwgtstor, weightcopy);
    1602             : 
    1603             : 
    1604             :   }
    1605           0 :   sumWeight.putStorage(sumwgtstor, issumWgtCopy); 
    1606           0 :   elWeight.freeStorage(wgtStorage,iswgtCopy);
    1607             :     
    1608           0 : }
    1609             : 
    1610           0 : void MosaicFT::get(vi::VisBuffer2& vb, Int row)
    1611             : {
    1612             :   
    1613             : 
    1614             :   
    1615             :   // If row is -1 then we pass through all rows
    1616             :   Int startRow, endRow, nRow;
    1617           0 :   if (row==-1) {
    1618           0 :     nRow=vb.nRows();
    1619           0 :     startRow=0;
    1620           0 :     endRow=nRow-1;
    1621             :     //  vb.modelVisCube()=Complex(0.0,0.0);
    1622             :   } else {
    1623           0 :     nRow=1;
    1624           0 :     startRow=row;
    1625           0 :     endRow=row;
    1626             :     //  vb.modelVisCube().xyPlane(row)=Complex(0.0,0.0);
    1627             :   }
    1628             :   
    1629             : 
    1630             :  
    1631             : 
    1632           0 :   matchChannel(vb);
    1633             :  
    1634             :   //No point in reading data if its not matching in frequency
    1635           0 :   if(max(chanMap)==-1)
    1636           0 :     return;
    1637             : 
    1638             :   // Get the uvws in a form that Fortran can use
    1639           0 :   Matrix<Double> uvw(negateUV(vb));
    1640           0 :   Vector<Double> dphase(vb.nRows());
    1641           0 :   dphase=0.0;
    1642             :  
    1643           0 :   doUVWRotation_p=true;
    1644           0 :   girarUVW(uvw, dphase, vb);
    1645           0 :   refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
    1646             :   
    1647             :   
    1648             :   
    1649             :  
    1650           0 :   Cube<Complex> data;
    1651           0 :   Cube<Int> flags;
    1652           0 :   getInterpolateArrays(vb, data, flags);
    1653             : 
    1654             :   //vi::VisBuffer2& vbRotuvw=const_cast<vi::VisBuffer2&>(vb);
    1655             :   //vbRotuvw.setUvw(uvw);
    1656             : 
    1657             :   //Need to get interpolated freqs
    1658           0 :   findConvFunction(*image, vb, uvw);
    1659             : 
    1660             :   // no valid pointing in this buffer
    1661           0 :   if(convSupport <= 0)
    1662           0 :     return;
    1663             :   Complex *datStorage;
    1664             :   Bool isCopy;
    1665           0 :   datStorage=data.getStorage(isCopy);
    1666             :   
    1667             : 
    1668           0 :   Vector<Int> rowFlags(vb.nRows());
    1669           0 :   rowFlags=0;
    1670           0 :   rowFlags(vb.flagRow())=true;
    1671           0 :   if(!usezero_p) {
    1672           0 :     for (Int rownr=startRow; rownr<=endRow; rownr++) {
    1673           0 :       if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
    1674             :     }
    1675             :   }
    1676           0 :   Int nvp=data.shape()[0];
    1677           0 :   Int nvc=data.shape()[1];
    1678           0 :   Int nvisrow=data.shape()[2];
    1679           0 :   Int csamp=convSampling;
    1680           0 :   Int csize=convSize;
    1681           0 :   Int csupp=convSupport;
    1682           0 :   Int nc=nchan;
    1683           0 :   Int np=npol;
    1684           0 :   Int nxp=nx;
    1685           0 :   Int nyp=ny;
    1686             :   Bool uvwcopy; 
    1687           0 :   const Double *uvwstor=uvw.getStorage(uvwcopy);
    1688           0 :   Int nPolConv=convFunc.shape()[2];
    1689           0 :   Int nChanConv=convFunc.shape()[3];
    1690           0 :   Int nConvFunc=convFunc.shape()(4);
    1691             :   ////////**************************
    1692           0 :   Cube<Int> loc(2, nvc, nRow);
    1693           0 :   Cube<Int> off(2, nvc, nRow);
    1694           0 :   Matrix<Complex> phasor(nvc, nRow);
    1695             :   Bool delphase;
    1696             :   Bool del;
    1697           0 :   const Int* pmapstor=polMap.getStorage(del);
    1698           0 :   const Int* cmapstor=chanMap.getStorage(del);
    1699           0 :   Complex * phasorstor=phasor.getStorage(delphase);
    1700           0 :   const Double * visfreqstor=interpVisFreq_p.getStorage(del);
    1701           0 :   const Double * scalestor=uvScale.getStorage(del);
    1702           0 :   const Double * offsetstor=uvOffset.getStorage(del);
    1703           0 :   const Int * flagstor=flags.getStorage(del);
    1704           0 :   const Int * rowflagstor=rowFlags.getStorage(del);
    1705           0 :   Int * locstor=loc.getStorage(del);
    1706           0 :   Int * offstor=off.getStorage(del);
    1707           0 :   const Double *dpstor=dphase.getStorage(del);
    1708           0 :   const Int *convrowmapstor=convRowMap_p.getStorage(del);
    1709           0 :   const Int *convchanmapstor=convChanMap_p.getStorage(del);
    1710           0 :   const Int *convpolmapstor=convPolMap_p.getStorage(del);
    1711             :   ////////***************************
    1712             : 
    1713             :   Int irow;
    1714           0 :   Int nth=1;
    1715             :  #ifdef _OPENMP
    1716           0 :   if(numthreads_p >0){
    1717           0 :     nth=min(numthreads_p, omp_get_max_threads());
    1718             :   }
    1719             :   else{   
    1720           0 :     nth= omp_get_max_threads();
    1721             :   }
    1722             :   //nth=min(4,nth);
    1723             : #endif
    1724             :  
    1725           0 :   Timer tim;
    1726           0 :   tim.mark();
    1727             : 
    1728           0 :    Int dow=0;
    1729           0 :    Double cinv=Double(1.0)/C::c;
    1730           0 : #pragma omp parallel default(none) private(irow) firstprivate(visfreqstor, nvc, scalestor, offsetstor, csamp, phasorstor, uvwstor, locstor, offstor, dpstor, dow, cinv) shared(startRow, endRow) num_threads(nth)  
    1731             : {
    1732             : #pragma omp for
    1733             :   for (irow=startRow; irow<=endRow;irow++){
    1734             :     /////////////////*locateuvw(uvwstor,dpstor, visfreqstor, nvc, scalestor, offsetstor, csamp, 
    1735             :     //    locstor, 
    1736             :                 ///////////           offstor, phasorstor, irow, false);
    1737             :     //using the fortran version which is significantly faster ...this can account for 10% less overall degridding time
    1738             :     locuvw(uvwstor, dpstor, visfreqstor, &nvc, scalestor, offsetstor, &csamp, locstor, 
    1739             :            offstor, phasorstor, &irow, &dow, &cinv);
    1740             :   }  
    1741             : 
    1742             :  }//end pragma parallel
    1743           0 :  Int rbeg=startRow+1;
    1744           0 :  Int rend=endRow+1;
    1745           0 :  Int npart=nth;
    1746             :  
    1747             :  Bool gridcopy;
    1748           0 :  const Complex *gridstor=griddedData.getStorage(gridcopy);
    1749             :  Bool convcopy;
    1750             :  ////Degridding needs the conjugate ...doing it here
    1751           0 :  Array<Complex> conjConvFunc=conj(convFunc);
    1752           0 :  const Complex *convstor=conjConvFunc.getStorage(convcopy);
    1753           0 :   Int ix=0;
    1754           0 : #pragma omp parallel default(none) private(ix, rbeg, rend) firstprivate(uvwstor, datStorage, flagstor, rowflagstor, convstor, pmapstor, cmapstor, gridstor, nxp, nyp, np, nc, csamp, csize, csupp, nvp, nvc, nvisrow, phasorstor, locstor, offstor, nPolConv, nChanConv, nConvFunc, convrowmapstor, convpolmapstor, convchanmapstor, npart)  num_threads(npart)
    1755             :   {
    1756             :     #pragma omp for schedule(dynamic) 
    1757             :     for (ix=0; ix< npart; ++ix){
    1758             :       rbeg=ix*(nvisrow/npart)+1;
    1759             :       rend=(ix != (npart-1)) ? (rbeg+(nvisrow/npart)-1) : (rbeg+(nvisrow/npart)-1+nvisrow%npart) ;
    1760             :       //cerr << "maps "  << convChanMap_p << "   " << chanMap  << endl;
    1761             :       //cerr << "nchan " << nchan << "  nchanconv " << nChanConv << " npolconv " << nPolConv << " nRowConv " << nConvFunc << endl;
    1762             :      sectdmos2(
    1763             :                datStorage,
    1764             :                &nvp,
    1765             :                &nvc,
    1766             :                flagstor,
    1767             :                rowflagstor,
    1768             :                &nvisrow,
    1769             :                gridstor,
    1770             :                &nxp,
    1771             :                &nyp,
    1772             :                &np,
    1773             :                &nc,
    1774             :                &csupp,
    1775             :                &csize,   
    1776             :                &csamp,
    1777             :                convstor,
    1778             :                cmapstor,
    1779             :                pmapstor,
    1780             :                convrowmapstor, convchanmapstor,
    1781             :                convpolmapstor,
    1782             :                &nConvFunc, &nChanConv, &nPolConv,
    1783             :                &rbeg, &rend, locstor, offstor, phasorstor
    1784             :                );
    1785             : 
    1786             : 
    1787             :     }
    1788             :   }//end pragma omp
    1789             : 
    1790             : 
    1791           0 :   data.putStorage(datStorage, isCopy);
    1792           0 :   griddedData.freeStorage(gridstor, gridcopy);
    1793           0 :   convFunc.freeStorage(convstor, convcopy);
    1794             :   
    1795           0 :    timedegrid_p+=tim.real();
    1796             : 
    1797           0 :   interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
    1798           0 : }
    1799             : 
    1800             : 
    1801             : /*
    1802             : void MosaicFT::get(VisBuffer& vb, Int row)
    1803             : {
    1804             :   
    1805             : 
    1806             :   
    1807             :   // If row is -1 then we pass through all rows
    1808             :   Int startRow, endRow, nRow;
    1809             :   if (row==-1) {
    1810             :     nRow=vb.nRow();
    1811             :     startRow=0;
    1812             :     endRow=nRow-1;
    1813             :     //  vb.modelVisCube()=Complex(0.0,0.0);
    1814             :   } else {
    1815             :     nRow=1;
    1816             :     startRow=row;
    1817             :     endRow=row;
    1818             :     //  vb.modelVisCube().xyPlane(row)=Complex(0.0,0.0);
    1819             :   }
    1820             :   
    1821             : 
    1822             :  
    1823             : 
    1824             : 
    1825             :   // Get the uvws in a form that Fortran can use
    1826             :   Matrix<Double> uvw(3, vb.uvw().nelements());
    1827             :   uvw=0.0;
    1828             :   Vector<Double> dphase(vb.uvw().nelements());
    1829             :   dphase=0.0;
    1830             :   //NEGATING to correct for an image inversion problem
    1831             :   for (Int i=startRow;i<=endRow;i++) {
    1832             :     for (Int idim=0;idim<2;idim++) uvw(idim,i)=-vb.uvw()(i)(idim);
    1833             :     uvw(2,i)=vb.uvw()(i)(2);
    1834             :   }
    1835             :   
    1836             :   doUVWRotation_p=true;
    1837             :   girarUVW(uvw, dphase, vb);
    1838             :   refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
    1839             :   
    1840             :   
    1841             :   //Check if ms has changed then cache new spw and chan selection
    1842             :   if(vb.newMS())
    1843             :     matchAllSpwChans(vb);
    1844             :   
    1845             :   //Here we redo the match or use previous match
    1846             :   
    1847             :   //Channel matching for the actual spectral window of buffer
    1848             :   if(doConversion_p[vb.spectralWindow()]){
    1849             :     matchChannel(vb.spectralWindow(), vb);
    1850             :   }
    1851             :   else{
    1852             :     chanMap.resize();
    1853             :     chanMap=multiChanMap_p[vb.spectralWindow()];
    1854             :   }
    1855             :   //No point in reading data if its not matching in frequency
    1856             :   if(max(chanMap)==-1)
    1857             :     return;
    1858             : 
    1859             :   Cube<Complex> data;
    1860             :   Cube<Int> flags;
    1861             :   getInterpolateArrays(vb, data, flags);
    1862             :   //Need to get interpolated freqs
    1863             :   findConvFunction(*image, vb);
    1864             :   // no valid pointing in this buffer
    1865             :   if(convSupport <= 0)
    1866             :     return;
    1867             :   Complex *datStorage;
    1868             :   Bool isCopy;
    1869             :   datStorage=data.getStorage(isCopy);
    1870             : 
    1871             : 
    1872             :   Vector<Int> rowFlags(vb.nRow());
    1873             :   rowFlags=0;
    1874             :   rowFlags(vb.flagRow())=true;
    1875             :   if(!usezero_p) {
    1876             :     for (Int rownr=startRow; rownr<=endRow; rownr++) {
    1877             :       if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
    1878             :     }
    1879             :   }
    1880             :   Int nPolConv=convFunc.shape()[2];
    1881             :   Int nChanConv=convFunc.shape()[3];
    1882             :   Int nConvFunc=convFunc.shape()(4);
    1883             :  
    1884             :   {
    1885             :     Bool del;
    1886             :      Bool uvwcopy; 
    1887             :      const Double *uvwstor=uvw.getStorage(uvwcopy);
    1888             :      Bool gridcopy;
    1889             :      const Complex *gridstor=griddedData.getStorage(gridcopy);
    1890             :      Bool convcopy;
    1891             :      ////Degridding needs the conjugate ...doing it here
    1892             :      Array<Complex> conjConvFunc=conj(convFunc);
    1893             :      const Complex *convstor=conjConvFunc.getStorage(convcopy);
    1894             :      //     IPosition s(data.shape());
    1895             :      const IPosition& fs=data.shape();
    1896             :      std::vector<Int> s(fs.begin(), fs.end());
    1897             :      //cerr << "maps "  << convChanMap_p << "   " << chanMap  << endl;
    1898             :      //cerr << "nchan " << nchan << "  nchanconv " << nChanConv << " npolconv " << nPolConv << " nRowConv " << nConvFunc << endl;
    1899             :      dmos2(uvwstor,
    1900             :             dphase.getStorage(del),
    1901             :             datStorage,
    1902             :             &s[0],
    1903             :             &s[1],
    1904             :             flags.getStorage(del),
    1905             :             rowFlags.getStorage(del),
    1906             :             &s[2],
    1907             :             &row,
    1908             :            uvScale.getStorage(del), //10
    1909             :             uvOffset.getStorage(del),
    1910             :             gridstor,
    1911             :             &nx,
    1912             :             &ny,
    1913             :             &npol,
    1914             :             &nchan,
    1915             :             interpVisFreq_p.getStorage(del),
    1916             :             &C::c,
    1917             :             &convSupport,
    1918             :            &convSize,   //20
    1919             :             &convSampling,
    1920             :             convstor,
    1921             :             chanMap.getStorage(del),
    1922             :             polMap.getStorage(del),
    1923             :             convRowMap_p.getStorage(del), convChanMap_p.getStorage(del),
    1924             :             convPolMap_p.getStorage(del),
    1925             :             &nConvFunc, &nChanConv, &nPolConv);
    1926             :       data.putStorage(datStorage, isCopy);
    1927             :      uvw.freeStorage(uvwstor, uvwcopy);
    1928             :      griddedData.freeStorage(gridstor, gridcopy);
    1929             :      convFunc.freeStorage(convstor, convcopy);
    1930             :   }
    1931             :   interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
    1932             : }
    1933             : 
    1934             : */
    1935             : 
    1936             : 
    1937             : // Finalize the FFT to the Sky. Here we actually do the FFT and
    1938             : // return the resulting image
    1939           0 : ImageInterface<Complex>& MosaicFT::getImage(Matrix<Float>& weights,
    1940             :                                             Bool normalize) 
    1941             : {
    1942             :   //AlwaysAssert(lattice, AipsError);
    1943           0 :   AlwaysAssert(image, AipsError);
    1944             : 
    1945           0 :    if((griddedData.nelements() ==0) && (griddedData2.nelements()==0)){
    1946           0 :     throw(AipsError("Programmer error ...request for image without right order of calls"));
    1947             :   }
    1948             : 
    1949           0 :   logIO() << LogOrigin("MosaicFT", "getImage") << LogIO::NORMAL;
    1950             :   
    1951           0 :   weights.resize(sumWeight.shape());
    1952           0 :   convertArray(weights, sumWeight);
    1953           0 :   SynthesisUtilMethods::getResource("mem peak in getImage");
    1954             :   
    1955             :   // If the weights are all zero then we cannot normalize
    1956             :   // otherwise we don't care.
    1957           0 :   if(max(weights)==0.0) {
    1958           0 :     if(normalize) {
    1959             :       logIO() << LogIO::SEVERE << "No useful data in MosaicFT: weights all zero"
    1960           0 :               << LogIO::POST;
    1961             :     }
    1962             :     else {
    1963             :       logIO() << LogIO::WARN << "No useful data in MosaicFT: weights all zero"
    1964           0 :               << LogIO::POST;
    1965             :     }
    1966             :   }
    1967             :   else {
    1968             :     
    1969             :     //const IPosition latticeShape = lattice->shape();
    1970             :     
    1971             :     logIO() << LogIO::DEBUGGING
    1972           0 :             << "Starting FFT and scaling of image" << LogIO::POST;
    1973           0 :     if(useDoubleGrid_p){
    1974           0 :       ArrayLattice<DComplex> darrayLattice(griddedData2);
    1975           0 :       ft_p.c2cFFT(darrayLattice,false);
    1976           0 :       griddedData.resize(griddedData2.shape());
    1977           0 :       convertArray(griddedData, griddedData2);
    1978             :       
    1979             :       //Don't need the double-prec grid anymore...
    1980           0 :       griddedData2.resize();
    1981           0 :       arrayLattice = new ArrayLattice<Complex>(griddedData);
    1982           0 :       lattice=arrayLattice;
    1983             : 
    1984           0 :     }
    1985             :     else{
    1986           0 :       arrayLattice = new ArrayLattice<Complex>(griddedData);
    1987           0 :       lattice=arrayLattice;
    1988           0 :       ft_p.c2cFFT(*lattice,false);
    1989             :     }
    1990             :    {////Do the grid correction
    1991           0 :       Int inx = lattice->shape()(0);
    1992             :       //Int iny = lattice->shape()(1);
    1993           0 :       Vector<Complex> correction(inx);
    1994           0 :       correction=Complex(1.0, 0.0);
    1995             : 
    1996             :       // Int npixCorr= max(nx,ny);
    1997           0 :       Vector<Float> sincConvX(nx);
    1998           0 :       for (Int ix=0;ix<nx;ix++) {
    1999           0 :         Float x=C::pi*Float(ix-nx/2)/(Float(nx)*Float(convSampling));
    2000           0 :         if(ix==nx/2) {
    2001           0 :           sincConvX(ix)=1.0;
    2002             :         }
    2003             :         else {
    2004           0 :           sincConvX(ix)=sin(x)/x;
    2005             :         }
    2006             :       }
    2007           0 :       Vector<Float> sincConvY(ny);
    2008           0 :       for (Int ix=0;ix<ny;ix++) {
    2009           0 :         Float x=C::pi*Float(ix-ny/2)/(Float(ny)*Float(convSampling));
    2010           0 :         if(ix==ny/2) {
    2011           0 :           sincConvY(ix)=1.0;
    2012             :         }
    2013             :         else {
    2014           0 :           sincConvY(ix)=sin(x)/x;
    2015             :         }
    2016             :       }
    2017             :     
    2018             : 
    2019           0 :       IPosition cursorShape(4, inx, 1, 1, 1);
    2020           0 :       IPosition axisPath(4, 0, 1, 2, 3);
    2021           0 :       LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
    2022           0 :       LatticeIterator<Complex> lix(*lattice, lsx);
    2023           0 :       for(lix.reset();!lix.atEnd();lix++) {
    2024           0 :         Int pol=lix.position()(2);
    2025           0 :         Int chan=lix.position()(3);
    2026           0 :         if(weights(pol, chan)!=0.0) {
    2027           0 :           Int iy=lix.position()(1);
    2028             :           //gridder->correctX1D(correction,iy);
    2029           0 :           for (Int ix=0;ix<nx;ix++) {
    2030           0 :             correction(ix)=1.0/(sincConvX(ix)*sincConvY(iy));
    2031             :           }
    2032           0 :           lix.rwVectorCursor()*=correction;
    2033           0 :           if(normalize) {
    2034           0 :             Complex rnorm(1.0/weights(pol,chan));
    2035           0 :             lix.rwCursor()*=rnorm;
    2036             :           }
    2037             :         }
    2038             :         else {
    2039           0 :           lix.woCursor()=0.0;
    2040             :         }
    2041             :       }
    2042           0 :     }
    2043             :      
    2044             :     
    2045             : 
    2046             :     //if(!isTiled) 
    2047             :     {
    2048           0 :       LatticeLocker lock1 (*(image), FileLocker::Write);
    2049             :       // Check the section from the image BEFORE converting to a lattice 
    2050           0 :       IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
    2051           0 :                     (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
    2052           0 :       IPosition stride(4, 1);
    2053           0 :       IPosition trc(blc+image->shape()-stride);
    2054             :       
    2055             :       // Do the copy
    2056           0 :       IPosition start(4, 0);
    2057           0 :       image->put(griddedData(blc, trc));
    2058           0 :     }
    2059             :   }
    2060           0 :   if(!arrayLattice.null()) arrayLattice=0;
    2061           0 :   if(!lattice.null()) lattice=0;
    2062           0 :    griddedData.resize();
    2063           0 :   image->clearCache();
    2064           0 :   return *image;
    2065             : }
    2066             : 
    2067             : // Get weight image
    2068           0 : void MosaicFT::getWeightImage(ImageInterface<Float>& weightImage,
    2069             :                               Matrix<Float>& weights) 
    2070             : {
    2071             :   
    2072           0 :   logIO() << LogOrigin("MosaicFT", "getWeightImage") << LogIO::NORMAL;
    2073             :   
    2074           0 :   weights.resize(sumWeight.shape());
    2075           0 :   convertArray(weights,sumWeight);
    2076             :   /*
    2077             :   weightImage.copyData((LatticeExpr<Float>) 
    2078             :                        (iif((pbConvFunc_p->getFluxScaleImage()) > (0.0), 
    2079             :                             (*skyCoverage_p),0.0)));
    2080             :   */
    2081           0 :   weightImage.copyData(*skyCoverage_p);
    2082             :   //cerr << "getWeightIm " << max(sumWeight) << "    " << max(skyCoverage_p->get()) << endl;
    2083             :   
    2084           0 :    skyCoverage_p->tempClose();
    2085             : 
    2086           0 : }
    2087             : 
    2088           0 : void MosaicFT::getFluxImage(ImageInterface<Float>& fluxImage) {
    2089             : 
    2090           0 :   if (stokes_p=="QU" || stokes_p=="IV"){
    2091           0 :     pbConvFunc_p->sliceFluxScale(2);
    2092             :   }
    2093           0 :   else if(stokes_p=="Q" || stokes_p=="U" ||  stokes_p=="V" || stokes_p=="I" ){
    2094           0 :      pbConvFunc_p->sliceFluxScale(1);
    2095             :   }
    2096           0 :    else if(stokes_p=="IQU"){
    2097           0 :        pbConvFunc_p->sliceFluxScale(3);
    2098             :    }
    2099           0 :   IPosition inShape=(pbConvFunc_p->getFluxScaleImage()).shape();
    2100           0 :   IPosition outShape=fluxImage.shape();
    2101           0 :   if(outShape==inShape){
    2102           0 :     fluxImage.copyData(pbConvFunc_p->getFluxScaleImage());
    2103             :   }
    2104           0 :   else if((outShape(0)==inShape(0)) && (outShape(1)==inShape(1)) 
    2105           0 :           && (outShape(2)==inShape(2))){
    2106             :     //case where CubeSkyEquation is chunking...copy the first pol-cube
    2107           0 :     IPosition cursorShape(4, inShape(0), inShape(1), inShape(2), 1);
    2108           0 :     IPosition axisPath(4, 0, 1, 2, 3);
    2109           0 :     LatticeStepper lsout(outShape, cursorShape, axisPath);
    2110           0 :     LatticeStepper lsin(inShape, cursorShape, axisPath);
    2111           0 :     LatticeIterator<Float> liout(fluxImage, lsout);
    2112           0 :     RO_LatticeIterator<Float> liin(pbConvFunc_p->getFluxScaleImage(), lsin);
    2113           0 :     liin.reset();
    2114           0 :     for(liout.reset();!liout.atEnd();liout++) {
    2115           0 :       if(inShape(2)==1)
    2116           0 :         liout.woMatrixCursor()=liin.matrixCursor();
    2117             :       else
    2118           0 :         liout.woCubeCursor()=liin.cubeCursor();
    2119             :     }
    2120             : 
    2121             : 
    2122           0 :   }
    2123             :   else{
    2124             :     //Should not reach here but the we're getting old
    2125           0 :     cout << "Bad case of shape mismatch in flux image shape" << endl;
    2126             :   }
    2127           0 : }
    2128             : 
    2129           0 : CountedPtr<TempImage<Float> >& MosaicFT::getConvWeightImage(){
    2130           0 :   if(!doneWeightImage_p)
    2131           0 :     finalizeToSky();
    2132           0 :   return skyCoverage_p;
    2133             : }
    2134             : 
    2135           0 : Bool MosaicFT::toRecord(String&  error,
    2136             :                         RecordInterface& outRec, Bool withImage, const String diskimage)
    2137             : {  
    2138             :   // Save the current MosaicFT object to an output state record
    2139           0 :   Bool retval = true;
    2140           0 :   if(!FTMachine::toRecord(error, outRec, withImage, diskimage))
    2141           0 :     return false;
    2142             :   
    2143           0 :   if(sj_p){
    2144           0 :     outRec.define("telescope", sj_p->telescope());
    2145             :     //cerr <<" Telescope " << sj_p->telescope() << endl;
    2146             :   }
    2147           0 :   outRec.define("uvscale", uvScale);
    2148           0 :   outRec.define("uvoffset", uvOffset);
    2149           0 :   outRec.define("cachesize", Int64(cachesize));
    2150           0 :   outRec.define("tilesize", tilesize);
    2151           0 :   outRec.define("maxabsdata", maxAbsData);
    2152           0 :   Vector<Int> center_loc(4), offset_loc(4);
    2153           0 :   for (Int k=0; k<4 ; k++){
    2154           0 :     center_loc(k)=centerLoc(k);
    2155           0 :     offset_loc(k)=offsetLoc(k);
    2156             :   }
    2157           0 :   outRec.define("centerloc", center_loc);
    2158           0 :   outRec.define("offsetloc", offset_loc);
    2159           0 :   outRec.define("usezero", usezero_p);
    2160           0 :   outRec.define("convfunc", convFunc);
    2161           0 :   outRec.define("weightconvfunc", weightConvFunc_p);
    2162           0 :   outRec.define("convsampling", convSampling);
    2163           0 :   outRec.define("convsize", convSize);
    2164           0 :   outRec.define("convsupport", convSupport);
    2165           0 :   outRec.define("convsupportplanes", convSupportPlanes_p);
    2166           0 :   outRec.define("convsizeplanes", convSizePlanes_p);
    2167           0 :   outRec.define("convRowMap",  convRowMap_p);
    2168           0 :   outRec.define("stokes", stokes_p);
    2169           0 :   outRec.define("useconjconvfunc", useConjConvFunc_p);
    2170           0 :   outRec.define("usepointingtable", usePointingTable_p);
    2171           0 :   if(!pbConvFunc_p.null()){
    2172           0 :     Record subRec;
    2173             :     //cerr << "Doing pbconvrec " << endl;
    2174           0 :     pbConvFunc_p->toRecord(subRec);
    2175           0 :     outRec.defineRecord("pbconvfunc", subRec);        
    2176           0 :   }
    2177             :   
    2178             : 
    2179           0 :   return retval;
    2180           0 : }
    2181             : 
    2182           0 : Bool MosaicFT::fromRecord(String& error,
    2183             :                           const RecordInterface& inRec)
    2184             : {
    2185           0 :   Bool retval = true;
    2186           0 :   pointingToImage=0;
    2187           0 :   doneWeightImage_p=false;
    2188           0 :   convWeightImage_p=nullptr;
    2189             :   
    2190           0 :   if(!FTMachine::fromRecord(error, inRec))
    2191           0 :     return false;
    2192           0 :   sj_p=0;
    2193           0 :   if(inRec.isDefined("telescope")){
    2194           0 :     String tel=inRec.asString("telescope");
    2195             :     PBMath::CommonPB pbtype;
    2196           0 :     Quantity freq(1e12, "Hz");// no useful band...just get default beam
    2197           0 :     String band="";
    2198           0 :     String pbname;
    2199           0 :     PBMath::whichCommonPBtoUse(tel, freq, band, pbtype, pbname);
    2200           0 :     if(pbtype != PBMath::UNKNOWN)
    2201           0 :       sj_p.reset(new VPSkyJones(tel,pbtype));
    2202           0 :   }
    2203             : 
    2204           0 :   inRec.get("name", machineName_p);
    2205           0 :   inRec.get("uvscale", uvScale);
    2206           0 :   inRec.get("uvoffset", uvOffset);
    2207           0 :   cachesize=inRec.asInt64("cachesize");
    2208           0 :   inRec.get("tilesize", tilesize);
    2209           0 :   inRec.get("maxabsdata", maxAbsData);
    2210           0 :   Vector<Int> center_loc(4), offset_loc(4);
    2211           0 :   inRec.get("centerloc", center_loc);
    2212           0 :   inRec.get("offsetloc", offset_loc);
    2213           0 :   uInt ndim4 = 4;
    2214           0 :   centerLoc=IPosition(ndim4, center_loc(0), center_loc(1), center_loc(2), 
    2215           0 :                       center_loc(3));
    2216           0 :   offsetLoc=IPosition(ndim4, offset_loc(0), offset_loc(1), offset_loc(2), 
    2217           0 :                       offset_loc(3));
    2218           0 :   imageCache=0; lattice=0; arrayLattice=0;
    2219           0 :   inRec.get("usezero", usezero_p);
    2220           0 :   inRec.get("convfunc", convFunc);
    2221           0 :   inRec.get("weightconvfunc", weightConvFunc_p);
    2222           0 :   inRec.get("convsampling", convSampling);
    2223           0 :   inRec.get("convsize", convSize);
    2224           0 :   inRec.get("convsupport", convSupport);
    2225           0 :   inRec.get("convsupportplanes", convSupportPlanes_p);
    2226           0 :   inRec.get("convsizeplanes", convSizePlanes_p);
    2227           0 :   inRec.get("convRowMap",  convRowMap_p);
    2228           0 :   inRec.get("stokes", stokes_p);
    2229           0 :   inRec.get("useconjconvfunc", useConjConvFunc_p);
    2230           0 :   inRec.get("usepointingtable", usePointingTable_p);
    2231           0 :   if(inRec.isDefined("pbconvfunc")){
    2232           0 :     Record subRec=inRec.asRecord("pbconvfunc");
    2233           0 :     String elname=subRec.asString("name");
    2234             :     // if we are predicting only ...no need to estimate fluxscale
    2235           0 :     if(elname=="HetArrayConvFunc"){
    2236             :     
    2237           0 :       pbConvFunc_p=new HetArrayConvFunc(subRec, !toVis_p);
    2238             :     }
    2239             :     else{
    2240           0 :       pbConvFunc_p=new SimplePBConvFunc(subRec, !toVis_p);
    2241           0 :       if(!sj_p)
    2242           0 :         throw(AipsError("Failed to recovermosaic FTmachine;\n If you are seeing this message when try to get model vis \n then either try to reset the model or use scratch column for now"));
    2243             :     }
    2244           0 :   }
    2245             :   else{
    2246           0 :     pbConvFunc_p=0;
    2247             :   }
    2248           0 :   gridder.reset(nullptr);
    2249           0 :    return retval;
    2250           0 : }
    2251             : 
    2252           0 : void MosaicFT::ok() {
    2253           0 :   AlwaysAssert(image, AipsError);
    2254           0 : }
    2255             : 
    2256             : // Make a plain straightforward honest-to-God image. This returns
    2257             : // a complex image, without conversion to Stokes. The representation
    2258             : // is that required for the visibilities.
    2259             : //----------------------------------------------------------------------
    2260           0 : void MosaicFT::makeImage(FTMachine::Type type, 
    2261             :                          vi::VisibilityIterator2& vi,
    2262             :                          ImageInterface<Complex>& theImage,
    2263             :                          Matrix<Float>& weight) {
    2264             :   
    2265             :   
    2266           0 :   logIO() << LogOrigin("MosaicFT", "makeImage") << LogIO::NORMAL;
    2267             :   
    2268           0 :   if(type==FTMachine::COVERAGE) {
    2269             :     logIO() << "Type COVERAGE not defined for Fourier transforms"
    2270           0 :             << LogIO::EXCEPTION;
    2271             :   }
    2272             :   
    2273             :   
    2274             : 
    2275             :   // Loop over all visibilities and pixels
    2276           0 :   vi::VisBuffer2* vb=vi.getVisBuffer();
    2277             :   
    2278             :   // Initialize put (i.e. transform to Sky) for this model
    2279           0 :   vi.origin();
    2280             :   
    2281           0 :   if(vb->polarizationFrame()==MSIter::Linear) {
    2282           0 :     StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::LINEAR);
    2283             :   }
    2284             :   else {
    2285           0 :     StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::CIRCULAR);
    2286             :   }
    2287             :   
    2288           0 :   initializeToSky(theImage,weight,*vb);
    2289             :   //This call is a NOP for all weighting schemes except for cube-briggs-perchanweightdensity
    2290           0 :   initBriggsWeightor(vi);
    2291             :   // Loop over the visibilities, putting VisBuffers
    2292           0 :   for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
    2293           0 :     for (vi.origin(); vi.more(); vi.next()) {
    2294             :       
    2295           0 :       switch(type) {
    2296           0 :       case FTMachine::RESIDUAL:
    2297           0 :         vb->setVisCube(vb->visCubeCorrected()-vb->visCubeModel());
    2298           0 :         put(*vb, -1, false);
    2299           0 :         break;
    2300           0 :       case FTMachine::MODEL:
    2301           0 :         vb->setVisCube(vb->visCubeModel());
    2302           0 :         put(*vb, -1, false);
    2303           0 :         break;
    2304           0 :       case FTMachine::CORRECTED:
    2305           0 :         vb->setVisCube(vb->visCubeCorrected());
    2306           0 :         put(*vb, -1, false);
    2307           0 :         break;
    2308           0 :       case FTMachine::PSF:
    2309           0 :         vb->setVisCube(Complex(1.0,0.0));
    2310           0 :         put(*vb, -1, true);
    2311           0 :         break;
    2312           0 :       case FTMachine::OBSERVED:
    2313             :       default:
    2314           0 :         put(*vb, -1, false);
    2315           0 :         break;
    2316             :       }
    2317             :     }
    2318             :   }
    2319           0 :   finalizeToSky();
    2320             :   // Normalize by dividing out weights, etc.
    2321           0 :   getImage(weight, true);
    2322           0 : }
    2323             : 
    2324           0 : Bool MosaicFT::getXYPos(const vi::VisBuffer2& vb, Int row) {
    2325             :   
    2326           0 :   MSColumns mscol(vb.ms());
    2327           0 :   const MSPointingColumns& act_mspc=mscol.pointing();
    2328           0 :   Int pointIndex=getIndex(act_mspc, vb.time()(row), vb.timeInterval()(row));
    2329           0 :   if((pointIndex<0)||pointIndex>=Int(act_mspc.time().nrow())) {
    2330             :     //    ostringstream o;
    2331             :     //    o << "Failed to find pointing information for time " <<
    2332             :     //      MVTime(vb.time()(row)/86400.0);
    2333             :     //    logIO_p << LogIO::DEBUGGING << String(o) << LogIO::POST;
    2334             :     //    logIO_p << String(o) << LogIO::POST;
    2335             :     
    2336           0 :     worldPosMeas = mscol.field().phaseDirMeas(vb.fieldId()(0));
    2337             :   }
    2338             :   else {
    2339             :    
    2340           0 :       worldPosMeas=act_mspc.directionMeas(pointIndex);
    2341             :       // Make a machine to convert from the worldPosMeas to the output
    2342             :       // Direction Measure type for the relevant frame
    2343             :  
    2344             : 
    2345             :  
    2346             :   }
    2347             : 
    2348           0 :   if(!pointingToImage) {
    2349             :     // Set the frame - choose the first antenna. For e.g. VLBI, we
    2350             :     // will need to reset the frame per antenna
    2351           0 :     FTMachine::mLocation_p=mscol.antenna().positionMeas()(0);
    2352           0 :     mFrame_p=MeasFrame(MEpoch(Quantity(vb.time()(row), "s")),
    2353           0 :                        FTMachine::mLocation_p);
    2354           0 :     MDirection::Ref outRef(directionCoord.directionType(), mFrame_p);
    2355           0 :     pointingToImage = new MDirection::Convert(worldPosMeas, outRef);
    2356             :   
    2357           0 :     if(!pointingToImage) {
    2358           0 :       logIO_p << "Cannot make direction conversion machine" << LogIO::EXCEPTION;
    2359             :     }
    2360           0 :   }
    2361             :   else {
    2362           0 :     mFrame_p.resetEpoch(MEpoch(Quantity(vb.time()(row), "s")));
    2363             :   }
    2364             :   
    2365           0 :   worldPosMeas=(*pointingToImage)(worldPosMeas);
    2366             :  
    2367           0 :   Bool result=directionCoord.toPixel(xyPos, worldPosMeas);
    2368           0 :   if(!result) {
    2369           0 :     logIO_p << "Failed to find pixel location for " 
    2370           0 :             << worldPosMeas.getAngle().getValue() << LogIO::EXCEPTION;
    2371           0 :     return false;
    2372             :   }
    2373           0 :   return result;
    2374             :   
    2375           0 : }
    2376             : // Get the index into the pointing table for this time. Note that the 
    2377             : // in the pointing table, TIME specifies the beginning of the spanned
    2378             : // time range, whereas for the main table, TIME is the centroid.
    2379             : // Note that the behavior for multiple matches is not specified! i.e.
    2380             : // if there are multiple matches, the index returned depends on the
    2381             : // history of previous matches. It is deterministic but not obvious.
    2382             : // One could cure this by searching but it would be considerably
    2383             : // costlier.
    2384           0 : Int MosaicFT::getIndex(const MSPointingColumns& mspc, const Double& time,
    2385             :                        const Double& /*interval*/) {
    2386           0 :   Int start=lastIndex_p;
    2387             :   // Search forwards
    2388           0 :   Int nrows=mspc.time().nrow();
    2389           0 :   if(nrows<1) {
    2390             :     //    logIO_p << "No rows in POINTING table - cannot proceed" << LogIO::EXCEPTION;
    2391           0 :     return -1;
    2392             :   }
    2393           0 :   for (Int i=start;i<nrows;i++) {
    2394           0 :     Double midpoint = mspc.time()(i); // time in POINTING table is midpoint
    2395             :     // If the interval in the pointing table is negative, use the last
    2396             :     // entry. Note that this may be invalid (-1) but in that case 
    2397             :     // the calling routine will generate an error
    2398           0 :     if(mspc.interval()(i)<0.0) {
    2399           0 :       return lastIndex_p;
    2400             :     }
    2401             :     // Pointing table interval is specified so we have to do a match
    2402             :     else {
    2403             :       // Is the midpoint of this pointing table entry within the specified
    2404             :       // tolerance of the main table entry?
    2405           0 :       if(abs(midpoint-time) < (mspc.interval()(i)/2.0)) {
    2406           0 :         lastIndex_p=i;
    2407           0 :         return i;
    2408             :       }
    2409             :     }
    2410             :   }
    2411             :   // Search backwards
    2412           0 :   for (Int i=start;i>=0;i--) {
    2413           0 :     Double midpoint = mspc.time()(i); // time in POINTING table is midpoint
    2414           0 :     if(mspc.interval()(i)<0.0) {
    2415           0 :       return lastIndex_p;
    2416             :     }
    2417             :     // Pointing table interval is specified so we have to do a match
    2418             :     else {
    2419             :       // Is the midpoint of this pointing table entry within the specified
    2420             :       // tolerance of the main table entry?
    2421           0 :       if(abs(midpoint-time) < (mspc.interval()(i)/2.0)) {
    2422           0 :         lastIndex_p=i;
    2423           0 :         return i;
    2424             :       }
    2425             :     }
    2426             :   }
    2427             :   // No match!
    2428           0 :   return -1;
    2429             : }
    2430             : 
    2431             : 
    2432             : 
    2433             : 
    2434           0 : void MosaicFT::addBeamCoverage(ImageInterface<Complex>& pbImage){
    2435             : 
    2436           0 :   CoordinateSystem cs(pbImage.coordinates());
    2437             :   //  IPosition blc(4,0,0,0,0);
    2438             :   //  IPosition trc(pbImage.shape());
    2439             :   //  trc(0)=trc(0)-1;
    2440             :   //  trc(1)=trc(1)-1;
    2441             :   // trc(2)=0;
    2442             :   //  trc(3)=0;
    2443           0 :   WCBox *wbox= new WCBox(LCBox(pbImage.shape()), cs);
    2444           0 :   SubImage<Float> toAddTo(*skyCoverage_p, ImageRegion(wbox), true);
    2445           0 :   TempImage<Float> beamStokes(pbImage.shape(), cs);
    2446           0 :   StokesImageUtil::To(beamStokes, pbImage);
    2447             :   //  toAddTo.copyData((LatticeExpr<Float>)(toAddTo + beamStokes ));
    2448           0 :   skyCoverage_p->copyData((LatticeExpr<Float>)(*skyCoverage_p + beamStokes ));
    2449             : 
    2450             : 
    2451           0 : }
    2452             : 
    2453             : 
    2454             : 
    2455             : 
    2456             : 
    2457           0 : String MosaicFT::name() const {
    2458           0 :   return machineName_p;
    2459             : }
    2460             : 
    2461             : } // REFIM ends
    2462             : } //# NAMESPACE CASA - END
    2463             : 
    2464             : 

Generated by: LCOV version 1.16