LCOV - code coverage report
Current view: top level - synthesis/TransformMachines - SimplePBConvFunc.cc (source / functions) Hit Total Coverage
Test: casacpp_coverage.info Lines: 0 508 0.0 %
Date: 2024-10-12 00:35:29 Functions: 0 20 0.0 %

          Line data    Source code
       1             : //# SimplePBConvFunc.cc: implementation of SimplePBConvFunc
       2             : //# Copyright (C) 2007
       3             : //# Associated Universities, Inc. Washington DC, USA.
       4             : //#
       5             : //# This library is free software; you can redistribute it and/or modify it
       6             : //# under the terms of the GNU Library General Public License as published by
       7             : //# the Free Software Foundation; either version 2 of the License, or (at your
       8             : //# option) any later version.
       9             : //#
      10             : //# This library is distributed in the hope that it will be useful, but WITHOUT
      11             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      12             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
      13             : //# License for more details.
      14             : //#
      15             : //# You should have received a copy of the GNU Library General Public License
      16             : //# along with this library; if not, write to the Free Software Foundation,
      17             : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
      18             : //#
      19             : //# Correspondence concerning AIPS++ should be adressed as follows:
      20             : //#        Internet email: casa-feedback@nrao.edu.
      21             : //#        Postal address: AIPS++ Project Office
      22             : //#                        National Radio Astronomy Observatory
      23             : //#                        520 Edgemont Road
      24             : //#                        Charlottesville, VA 22903-2475 USA
      25             : //#
      26             : //#
      27             : //# $Id$
      28             : #include <casacore/casa/Arrays/ArrayMath.h>
      29             : #include <casacore/casa/Arrays/Array.h>
      30             : #include <casacore/casa/Arrays/MaskedArray.h>
      31             : #include <casacore/casa/Arrays/Vector.h>
      32             : #include <casacore/casa/Arrays/Slice.h>
      33             : #include <casacore/casa/Arrays/Slicer.h>
      34             : #include <casacore/casa/Arrays/Matrix.h>
      35             : #include <casacore/casa/Arrays/Cube.h>
      36             : #include <casacore/casa/OS/Timer.h>
      37             : #include <casacore/casa/Utilities/Assert.h>
      38             : #include <casacore/casa/Quanta/MVTime.h>
      39             : #include <casacore/casa/Quanta/MVAngle.h>
      40             : #include <casacore/measures/Measures/MeasTable.h>
      41             : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
      42             : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
      43             : 
      44             : #include <casacore/images/Images/ImageInterface.h>
      45             : #include <casacore/images/Images/PagedImage.h>
      46             : #include <casacore/images/Images/TempImage.h>
      47             : #include <casacore/images/Images/SubImage.h>
      48             : #include <casacore/casa/Logging/LogIO.h>
      49             : #include <casacore/casa/Logging/LogSink.h>
      50             : #include <casacore/casa/Logging/LogMessage.h>
      51             : 
      52             : #include <casacore/ms/MeasurementSets/MSColumns.h>
      53             : #include <casacore/lattices/Lattices/ArrayLattice.h>
      54             : #include <casacore/lattices/Lattices/SubLattice.h>
      55             : #include <casacore/lattices/LRegions/LCBox.h>
      56             : #include <casacore/lattices/LEL/LatticeExpr.h>
      57             : #include <casacore/lattices/Lattices/LatticeCache.h>
      58             : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
      59             : 
      60             : #include <casacore/scimath/Mathematics/ConvolveGridder.h>
      61             : #include <casacore/scimath/Mathematics/MathFunc.h>
      62             : #include <msvis/MSVis/VisBuffer.h>
      63             : #include <msvis/MSVis/VisibilityIterator.h>
      64             : 
      65             : #include <synthesis/TransformMachines/SimplePBConvFunc.h>
      66             : #include <synthesis/TransformMachines/SkyJones.h>
      67             : #include <synthesis/TransformMachines/PBMath1DNumeric.h>
      68             : #include <casacore/scimath/Mathematics/FFTPack.h>
      69             : #include <casacore/scimath/Mathematics/FFTW.h>
      70             : #include <casacore/scimath/Mathematics/FFTServer.h>
      71             : #include <casacore/casa/Utilities/CompositeNumber.h>
      72             : #include <math.h>
      73             : #ifdef _OPENMP
      74             : #include <omp.h>
      75             : #endif
      76             : using namespace casacore;
      77             : namespace casa { //# NAMESPACE CASA - BEGIN
      78             : 
      79           0 : SimplePBConvFunc::SimplePBConvFunc(): nchan_p(-1),
      80           0 :         npol_p(-1), pointToPix_p(), directionIndex_p(-1), thePix_p(0),
      81           0 :         filledFluxScale_p(false),doneMainConv_p(0),
      82             :                                       
      83           0 :         calcFluxScale_p(true), actualConvIndex_p(-1), convSize_p(0), convSupport_p(0), pointingPix_p()  {
      84             :     //
      85             : 
      86           0 :     pbClass_p=PBMathInterface::COMMONPB;
      87           0 :     ft_p=FFT2D(true);
      88           0 : }
      89             : 
      90           0 :   SimplePBConvFunc::SimplePBConvFunc(const PBMathInterface::PBClass typeToUse): 
      91           0 :     nchan_p(-1),npol_p(-1),pointToPix_p(),
      92           0 :     directionIndex_p(-1), thePix_p(0), filledFluxScale_p(false),doneMainConv_p(0), 
      93           0 :      calcFluxScale_p(true), actualConvIndex_p(-1), convSize_p(0), convSupport_p(0), pointingPix_p() {
      94             :     //
      95           0 :     pbClass_p=typeToUse;
      96           0 :     ft_p=FFT2D(true);
      97           0 :   }
      98           0 :   SimplePBConvFunc::SimplePBConvFunc(const RecordInterface& rec, const Bool calcfluxneeded)
      99           0 :   : nchan_p(-1),npol_p(-1),pointToPix_p(), directionIndex_p(-1), thePix_p(0), filledFluxScale_p(false),
     100           0 :     doneMainConv_p(0), 
     101           0 :     calcFluxScale_p(calcfluxneeded), actualConvIndex_p(-1), convSize_p(0), convSupport_p(0), pointingPix_p()
     102             :   {
     103           0 :     String err;
     104           0 :     fromRecord(err, rec, calcfluxneeded);
     105           0 :     ft_p=FFT2D(true);
     106           0 :   }
     107           0 :   SimplePBConvFunc::~SimplePBConvFunc(){
     108             :     //
     109             :    
     110             : 
     111           0 :   }
     112             : 
     113           0 :   void SimplePBConvFunc::storeImageParams(const ImageInterface<Complex>& iimage,
     114             :                                           const VisBuffer& vb){
     115             :     //image signature changed...rather simplistic for now
     116           0 :     if((iimage.shape().product() != nx_p*ny_p*nchan_p*npol_p) || nchan_p < 1){
     117           0 :       csys_p=iimage.coordinates();
     118           0 :       Int coordIndex=csys_p.findCoordinate(Coordinate::DIRECTION);
     119           0 :       AlwaysAssert(coordIndex>=0, AipsError);
     120           0 :       directionIndex_p=coordIndex;
     121           0 :       dc_p=csys_p.directionCoordinate(directionIndex_p);
     122           0 :       ObsInfo imInfo=csys_p.obsInfo();
     123           0 :       String tel= imInfo.telescope();
     124           0 :       MPosition pos;
     125           0 :       if (vb.msColumns().observation().nrow() > 0) {
     126           0 :         tel = vb.msColumns().observation().telescopeName()(vb.msColumns().observationId()(0));
     127             :       }
     128           0 :       if (tel.length() == 0 || !tel.contains("VLA") ||
     129           0 :           !MeasTable::Observatory(pos,tel)) {
     130             :         // unknown observatory, use first antenna
     131           0 :           Int ant1=vb.antenna1()(0);
     132           0 :           pos=vb.msColumns().antenna().positionMeas()(ant1);
     133             :       }
     134             :       //cout << "TELESCOPE " << tel << endl;
     135             :       //Store this to build epochs via the time access of visbuffer later
     136           0 :       timeMType_p=MEpoch::castType(vb.msColumns().timeMeas()(0).getRef().getType());
     137           0 :       timeUnit_p=Unit(vb.msColumns().timeMeas().measDesc().getUnits()(0).getName());
     138             :       // timeUnit_p=Unit("s");
     139             :       //cout << "UNIT " << timeUnit_p.getValue() << " name " << timeUnit_p.getName()  << endl;
     140           0 :       pointFrame_p=MeasFrame(imInfo.obsDate(), pos);
     141           0 :       MDirection::Ref elRef(dc_p.directionType(), pointFrame_p);
     142             :       //For now we set the conversion from this direction 
     143           0 :       pointToPix_p=MDirection::Convert( MDirection(), elRef);
     144           0 :       nx_p=iimage.shape()(coordIndex);
     145           0 :       ny_p=iimage.shape()(coordIndex+1);
     146           0 :       pointingPix_p.resize(nx_p, ny_p);
     147           0 :       pointingPix_p.set(false);
     148           0 :       coordIndex=csys_p.findCoordinate(Coordinate::SPECTRAL);
     149           0 :       Int pixAxis=csys_p.pixelAxes(coordIndex)[0];
     150           0 :       nchan_p=iimage.shape()(pixAxis);
     151           0 :       coordIndex=csys_p.findCoordinate(Coordinate::STOKES);
     152           0 :       pixAxis=csys_p.pixelAxes(coordIndex)[0];
     153           0 :       npol_p=iimage.shape()(pixAxis);
     154           0 :       if(calcFluxScale_p){
     155           0 :           if(fluxScale_p.shape().nelements()==0){
     156           0 :                   fluxScale_p=TempImage<Float>(IPosition(4,nx_p,ny_p,npol_p,nchan_p), csys_p);
     157           0 :                   fluxScale_p.set(0.0);
     158             :           }
     159           0 :           filledFluxScale_p=false;
     160             :       }
     161             :       
     162           0 :     }
     163             : 
     164           0 :   }
     165             : 
     166           0 :   void SimplePBConvFunc::toPix(const VisBuffer& vb){
     167           0 :     thePix_p.resize(2);
     168             : 
     169           0 :     const MDirection& p1=pointingDirAnt1(vb);
     170           0 :     if(dc_p.directionType() !=  MDirection::castType(p1.getRef().getType())){
     171             :       //pointToPix_p.setModel(theDir);
     172             : 
     173           0 :         String tel= csys_p.obsInfo().telescope();
     174           0 :         if(!tel.contains("VLA")) {
     175             :                 //use first antenna as direction1_p is used to calculate pointing
     176             :                 // as only VLA uses observatory pos for calculations
     177           0 :                   Int ant1=vb.antenna1()(0);
     178           0 :                   MPosition pos=vb.msColumns().antenna().positionMeas()(ant1);
     179           0 :                   pointFrame_p.resetPosition(pos);
     180           0 :         }
     181           0 :       MEpoch timenow(Quantity(vb.time()(0), timeUnit_p), timeMType_p);
     182             :       //cerr << "Ref " << p1.toString() <<  " ep " << timenow.getRefString() << " time " << MVTime(timenow.getValue().getTime()).string(MVTime::YMD) << endl; 
     183           0 :       pointFrame_p.resetEpoch(timenow);
     184             :       ///////////////////////////
     185             :       //MDirection some=pointToPix_p(vb.direction1()(0));
     186             :       //MVAngle mvRa=some.getAngle().getValue()(0);
     187             :       //MVAngle mvDec=some.getAngle().getValue()(1);
     188             :       
     189             :       //cout  << mvRa(0.0).string(MVAngle::TIME,8) << "   ";
     190             :       // cout << mvDec.string(MVAngle::DIG2,8) << "   ";
     191             :       //cout << MDirection::showType(some.getRefPtr()->getType()) << endl;
     192             : 
     193             :       //////////////////////////
     194             :       //pointToPix holds pointFrame_p by reference...
     195             :       //thus good to go for conversion
     196           0 :       direction1_p=pointToPix_p(p1);
     197             :       //cerr << "converted " << direction1_p.toString() << endl;
     198             :       //direction2_p=pointToPix_p(vb.direction2()(0));
     199           0 :       direction2_p=direction1_p;
     200           0 :       dc_p.toPixel(thePix_p, direction1_p);
     201             : 
     202           0 :     }
     203             :     else{
     204           0 :       direction1_p=p1;
     205             :       //direction2_p=vb.direction2()(0);
     206             :       //For now 
     207           0 :       direction2_p=direction1_p;
     208           0 :       dc_p.toPixel(thePix_p, direction1_p);
     209             :     }
     210           0 :   }
     211             : 
     212           0 :   void SimplePBConvFunc::setWeightImage(CountedPtr<TempImage<Float> >& wgtimage){
     213           0 :     convWeightImage_p=wgtimage;
     214           0 :     calcFluxScale_p=true;
     215             : 
     216           0 :   }
     217             :  
     218           0 :   void SimplePBConvFunc::reset(){
     219           0 :     doneMainConv_p.resize();
     220           0 :     convFunctions_p.resize(0, true);
     221           0 :     convWeights_p.resize(0, true);
     222           0 :     convSizes_p.resize(0, true);
     223           0 :     convSupportBlock_p.resize(0, true);
     224           0 :     convFunctionMap_p.clear();
     225           0 :     vbConvIndex_p.clear();
     226           0 :   }
     227             : 
     228             : 
     229             : 
     230           0 :   Int SimplePBConvFunc::convIndex(const VisBuffer& vb){
     231           0 :           String elkey=String::toString(vb.msId())+String("_")+String::toString(vb.spectralWindow());
     232           0 :           if(vbConvIndex_p.count(elkey) > 0){
     233           0 :                   return vbConvIndex_p[elkey];
     234             :           }
     235           0 :           Int val=vbConvIndex_p.size();
     236           0 :           vbConvIndex_p[elkey]=val;
     237           0 :           return val;
     238           0 :   }
     239             : 
     240           0 :   const MDirection& SimplePBConvFunc::pointingDirAnt1(const VisBuffer& vb){
     241           0 :     std::ostringstream oss;
     242           0 :     oss << vb.msId() << "_" << vb.antenna1()(0) << "_";
     243           0 :     oss.precision(13);
     244           0 :     oss << vb.time()(0);
     245           0 :     String elkey=oss.str();
     246             :     //  String elkey=String::toString(vb.msId())+String("_")+String::toString(vb.antenna1()(0))+String("_")
     247             :     //                                                                    +String::toString(vb.time()(0));
     248             : 
     249             :     //cerr << "key " << elkey << " count " << ant1PointVal_p.count(elkey)  << " size " << ant1PointVal_p.size() << "  " << ant1PointingCache_p.nelements() << endl;
     250           0 :     if(ant1PointVal_p.count(elkey) > 0){
     251           0 :       return ant1PointingCache_p[ant1PointVal_p[elkey]];
     252             : 
     253             :     }
     254           0 :     Int val=ant1PointingCache_p.nelements();
     255           0 :     ant1PointingCache_p.resize(val+1, true);
     256           0 :     ant1PointingCache_p[val]=vb.firstDirection1();
     257             :     //cerr << "POINTINGDIRANT1 " << vb.firstDirection1() << "   " << vb.direction1()(0) << endl;
     258             :     //ant1PointingCache_p[val]=vbUtil_p.getPointingDir(vb, vb.antenna1()(0), 0);
     259           0 :     ant1PointVal_p[elkey]=val;
     260           0 :     return ant1PointingCache_p[val];
     261             : 
     262           0 :   }
     263           0 : void SimplePBConvFunc::findConvFunction(const ImageInterface<Complex>& iimage, 
     264             :                                         const VisBuffer& vb,
     265             :                                         const Int& convSampling,
     266             :                                         const Vector<Double>& visFreq, 
     267             :                                           Array<Complex>& convFunc, 
     268             :                                           Array<Complex>& weightConvFunc, 
     269             :                                           Vector<Int>& convsize,
     270             :                                           Vector<Int>& convSupport,
     271             :                                           Vector<Int>& convFuncPolMap,
     272             :                                           Vector<Int>& convFuncChanMap,
     273             :                                           Vector<Int>& convFuncRowMap,
     274             :                                           const Bool /*getConjFunc*/
     275             :                                           ){
     276             : 
     277             : 
     278             : 
     279           0 :     LogIO os;
     280           0 :      os << LogOrigin("SimplePBConv", "findConvFunction")  << LogIO::NORMAL;
     281             :      /////////////////////////
     282             :      os<< LogIO::DEBUG1 
     283           0 :        << "msID " << vb.msId()  <<  " ANT1 id" << vb.antenna1()(0) 
     284           0 :        << " direction " << vb.firstDirection1().toString() << " ANT2 id" 
     285           0 :        << vb.antenna2()(0) << " direction " << vb.direction2()(0).toString() 
     286           0 :        << LogIO::POST ; 
     287             :     //////////////////////
     288           0 :   Int convSamp=2*convSampling;
     289           0 :   storeImageParams(iimage, vb);
     290           0 :   convFuncChanMap.resize(vb.nChannel());
     291           0 :   Vector<Double> beamFreqs;
     292           0 :   findUsefulChannels(convFuncChanMap, beamFreqs, vb, visFreq);
     293             :   //cerr << "CHANMAP " << convFuncChanMap << " beamFreqs " << beamFreqs << endl;
     294           0 :   Int nBeamChans=beamFreqs.nelements();
     295             :   //indgen(convFuncChanMap);
     296           0 :   convFuncPolMap.resize(vb.nCorr());
     297           0 :   convFuncPolMap.set(0);
     298             :   //Only one plane in this version
     299           0 :   convFuncRowMap.resize();
     300           0 :   convFuncRowMap=Vector<Int>(vb.nRow(),0);
     301             :   //break reference
     302           0 :   convFunc.resize();
     303           0 :   weightConvFunc.resize();
     304             :  
     305             :   
     306             :   // Get the coordinate system
     307           0 :   CoordinateSystem coords(iimage.coordinates());
     308             :   
     309             :   
     310           0 :   actualConvIndex_p=convIndex(vb);
     311             :   //cerr << "In findConv " << actualConvIndex_p << endl;
     312             :   // Make a two dimensional image to calculate the
     313             :   // primary beam. We want this on a fine grid in the
     314             :   // UV plane 
     315           0 :   Int directionIndex=directionIndex_p;
     316           0 :     AlwaysAssert(directionIndex>=0, AipsError);
     317             :     
     318             :     // Set up the convolution function.
     319           0 :     Int nx=nx_p;
     320           0 :     Int ny=ny_p;
     321             :     //    convSize_p=max(nx,ny)*convSampling;
     322             :     //cerr << "size " << nx << "  " << ny << endl;
     323             :     //3 times the support size
     324           0 :     if(doneMainConv_p.shape()[0] < (actualConvIndex_p+1)){
     325             :       // cerr << "resizing DONEMAIN " <<   doneMainConv_p.shape()[0] << endl;
     326           0 :       doneMainConv_p.resize(actualConvIndex_p+1, true);
     327           0 :       doneMainConv_p[actualConvIndex_p]=false;
     328             :     }
     329             : 
     330           0 :     if(!(doneMainConv_p[actualConvIndex_p])){
     331             : 
     332             :       //convSize_p=4*(sj_p->support(vb, coords));
     333           0 :       convSize_p=Int(max(nx_p, ny_p)/2)*2*convSamp;
     334             :       //cerr << "convSize_p " << convSize_p << " support " << sj_p->support(vb, coords) << endl;
     335             :       // Make this a nice composite number, to speed up FFTs
     336             :       //cerr << "convSize_p 0 " <<  convSize_p << " convSamp " << convSamp<< endl;
     337           0 :       CompositeNumber cn(uInt(convSize_p*2.0));  
     338             :      
     339           0 :       convSize_p  = cn.nextLargerEven(Int(convSize_p));
     340             :       //cerr << "convSize : " << convSize_p << endl;
     341             : 
     342           0 :     }
     343             :     
     344             :    
     345           0 :     toPix(vb);
     346             :     //Timer tim;
     347             :     //tim.mark();
     348           0 :     addPBToFlux(vb);
     349             :     //tim.show("After addPBToFlux");
     350           0 :     DirectionCoordinate dc=dc_p;
     351             : 
     352             :     //where in the image in pixels is this pointing
     353           0 :     Vector<Double> pixFieldDir(2);
     354           0 :     pixFieldDir=thePix_p;
     355             : 
     356             :     //cerr << "pix of pointing " << pixFieldDir << endl;
     357           0 :     MDirection fieldDir=direction1_p;
     358             :     //shift from center
     359           0 :     pixFieldDir(0) = pixFieldDir(0) - Double(nx / 2);
     360           0 :     pixFieldDir(1) = pixFieldDir(1) - Double(ny / 2);
     361             : 
     362             :     //phase gradient per pixel to apply
     363           0 :     pixFieldDir(0) = -pixFieldDir(0)*2.0*C::pi/Double(nx)/Double(convSampling);
     364           0 :     pixFieldDir(1) = -pixFieldDir(1)*2.0*C::pi/Double(ny)/Double(convSampling);
     365             : 
     366             :     //cerr << "DonemainConv " << doneMainConv_p[actualConvIndex_p] << endl;
     367           0 :     if(!doneMainConv_p[actualConvIndex_p]){
     368           0 :       Vector<Double> sampling;
     369           0 :       sampling = dc.increment();
     370           0 :       sampling*=Double(convSamp);
     371           0 :       sampling(0)*=Double(nx)/Double(convSize_p);
     372           0 :       sampling(1)*=Double(ny)/Double(convSize_p);
     373           0 :       dc.setIncrement(sampling);
     374             :       
     375             :       
     376           0 :       Vector<Double> unitVec(2);
     377           0 :       unitVec=convSize_p/2;
     378           0 :       dc.setReferencePixel(unitVec);
     379             :       
     380             :       
     381             :       //make sure we are using the same units
     382           0 :       fieldDir.set(dc.worldAxisUnits()(0));
     383             :       //cerr << "SPB direction " << fieldDir << endl;
     384             : 
     385           0 :       dc.setReferenceValue(fieldDir.getAngle().getValue());
     386             :       
     387           0 :       coords.replaceCoordinate(dc, directionIndex);
     388           0 :       Int spind=coords.findCoordinate(Coordinate::SPECTRAL);
     389           0 :       SpectralCoordinate spCoord=coords.spectralCoordinate(spind);
     390           0 :       spCoord.setReferencePixel(Vector<Double>(1,0.0));
     391           0 :       spCoord.setReferenceValue(Vector<Double>(1, beamFreqs(0)));
     392           0 :       if(beamFreqs.nelements() >1)
     393           0 :         spCoord.setIncrement(Vector<Double>(1, beamFreqs(1)-beamFreqs(0)));
     394           0 :       coords.replaceCoordinate(spCoord, spind);
     395             : 
     396             : 
     397           0 :       CoordinateSystem coordLastPlane= coords;
     398           0 :       spCoord.setReferenceValue(Vector<Double>(1, beamFreqs(nBeamChans-1)));
     399           0 :       coordLastPlane.replaceCoordinate(spCoord, spind);
     400             :       //cerr << "BEAM freqs " << beamFreqs << endl;
     401             : 
     402             :       //  coords.list(logIO(), MDoppler::RADIO, IPosition(), IPosition());
     403             :       
     404           0 :       Int tempConvSize=((convSize_p/4/(convSamp/convSampling))/2)*2;
     405           0 :       IPosition pbShape(4, tempConvSize, tempConvSize, 1, nBeamChans);
     406             :       
     407           0 :       Long memtot=HostInfo::memoryFree();
     408           0 :       Double memtobeused= Double(memtot)*1024.0;
     409             :       //cerr << "Mem to be used " << memtobeused << " arraysize " << Double(convSize_p*convSize_p)*8.0 << endl;;
     410             :       //check for 32 bit OS and limit it to 2Gbyte
     411             :       if( sizeof(void*) == 4){
     412             :           if(memtot > 2000000)
     413             :                   memtot=2000000;
     414             :       }
     415           0 :       if(memtot <= 2000000)
     416           0 :           memtobeused=0.0;
     417             :       //cerr << "mem to be used " << memtobeused << endl;
     418             :       //tim.mark();
     419           0 :       IPosition start(4, 0, 0, 0, 0);
     420             :       //IPosition pbSlice(4, convSize_p, convSize_p, 1, 1);
     421             :       //cerr << "pbshape " << pbShape << endl;
     422           0 :       TempImage<Complex> twoDPB(TiledShape(pbShape, IPosition(4, pbShape(0), pbShape(1), 1, 1)), coords, memtobeused/10.0);
     423             : 
     424             :       //tim.show("after making one image");
     425           0 :       convFunc_p.resize(tempConvSize, tempConvSize);
     426           0 :       convFunc_p=0.0;
     427             :       
     428             :       
     429             : 
     430             :       // Accumulate terms 
     431             :       //Matrix<Complex> screen(convSize_p, convSize_p);
     432             :       //screen=1.0;
     433             :       // Either the SkyJones
     434             :       //tim.mark();
     435             :       //twoDPB.set(Complex(1.0,0.0));
     436             :       //for (Int k=0; k < nBeamChans; ++k){
     437             :       //blcin[3]=k;
     438             :       //trcin[3]=k;
     439             :       //Slicer slin(blcin, trcin, Slicer::endIsLast);
     440             :       //SubImage<Complex> subim(twoDPB, slin, true);
     441           0 :       TempImage<Complex> subim(IPosition(4, convSize_p, convSize_p, 1, 1), coordLastPlane, memtobeused/2.2);
     442           0 :       subim.set(Complex(1.0,0.0));
     443             :      
     444             : 
     445             :       /////Let's make a spheroid taper...
     446             :       /* Vector<Float> valsph(300);
     447             :       valsph.set(1.0);
     448             :       Sph_Conv<Float> d(1.0, 1.0);
     449             :        for(Int k=200; k <300; ++k){
     450             :         valsph(k)=d.value((Float)(k-200)/100.0);
     451             :        }
     452             :        Quantity fulrad(27,"arcmin");
     453             :        Quantity lefreq(1.4, "GHz");
     454             :        PBMath1DNumeric taper(valsph,fulrad,lefreq);
     455             :        taper.applyPB(subim, subim, fieldDir); 
     456             :       */
     457             :       ///////////////
     458             :        /*{
     459             :          PagedImage<Complex> thisScreen(subim.shape(), subim.coordinates(), "ELTAPER");
     460             :         //LatticeExpr<Float> le(abs(twoDPB2));
     461             :          thisScreen.copyData(subim);
     462             :          }*/
     463             :       //twoDPB.putSlice(screen, start);
     464             :       /////////////
     465             :       //Double wtime0=omp_get_wtime();
     466             :       //////////
     467           0 :       sj_p->apply(subim, subim, vb, 0);
     468             :       ///////////////
     469             :       //Double wtime1=omp_get_wtime();
     470             :       /////////////////
     471             :       /*{
     472             :          PagedImage<Complex> thisScreen(subim.shape(), subim.coordinates(), "BEAM-ELTAPER");
     473             :         //LatticeExpr<Float> le(abs(twoDPB2));
     474             :          thisScreen.copyData(subim);
     475             :          }*/
     476             :       //LatticeFFT::cfft2d(subim);
     477             :       //////
     478             :    
     479             : 
     480           0 :       ft_p.c2cFFT(subim);
     481             :      
     482             :       //   cerr << "make pb " << wtime1-wtime0 << " fft " << omp_get_wtime()-wtime1 << endl;
     483             :         //  }
     484             :       //tim.show("after an apply" );
     485             :       //tim.mark();
     486           0 :       TempImage<Float> screen2(TiledShape(IPosition(4, convSize_p, convSize_p, 1, 1)), coordLastPlane, memtobeused/2.2);
     487           0 :       screen2.set(1.0);
     488           0 :       TempImage<Complex> subout(TiledShape(IPosition(4, convSize_p, convSize_p, 1, 1)), coordLastPlane, memtobeused/2.2);
     489             :       
     490             :       //////Making a reference on half of the lattice as on the Mac rcfft is failing for some 
     491             :       //////reason
     492           0 :       SubImage<Complex> halfsubout(subout, Slicer(IPosition(4, 0, 0, 0, 0), 
     493           0 :                                                   IPosition(4, convSize_p/2, convSize_p-1, 0, 0), 
     494           0 :                                                   Slicer::endIsLast), true);
     495             :       /////////////////////
     496             :       //wtime0=omp_get_wtime();
     497             :       ///////////////////
     498           0 :       sj_p->applySquare(screen2, screen2, vb, 0); 
     499             :       ///////////////
     500             :       //wtime1=omp_get_wtime();
     501             :       /////////////////
     502             :       /////////
     503             :       /*  arr.resize();
     504             :       Array<Float> arrf;
     505             :       screen2.get(arrf, false);
     506             :       isRef=halfsubout.get (arr,false);
     507             :       Bool delf;
     508             :       Float* scrf=arrf.getStorage(delf);
     509             :       scr= arr.getStorage(del);
     510             :       
     511             :       ft_p.r2cFFT(scr, scrf, Long(convSize_p), Long(convSize_p));
     512             :       cerr << "isRef " << isRef << " iscopy " << del << " delf " << delf << endl; 
     513             :       if(del)
     514             :         arr.putStorage(scr, del);
     515             :       if(delf)
     516             :         arrf.putStorage(scrf, delf);
     517             :       if(!isRef)
     518             :         halfsubout.put(arr.reform(IPosition(4, convSize_p/2+1, convSize_p, 1, 1)));
     519             :       */
     520             :       ////////////////////
     521           0 :       ft_p.r2cFFT(halfsubout, screen2);
     522             :       //LatticeFFT::rcfft(halfsubout, screen2, true, false);
     523             :       //cerr << "make pb2 " << wtime1-wtime0 << " fft " << omp_get_wtime()-wtime1 << endl;
     524             :       //Real FFT fills only first half of the array
     525             :       //making it look like a Complex to Complex FFT
     526           0 :       IPosition iblc(4, 0, 3*subout.shape()(1)/8, 0, 0);
     527           0 :       IPosition itrc(4, 0, 5*subout.shape()(1)/8, 0, 0);
     528           0 :       for(Int x=subout.shape()(0)/2; x <(5*subout.shape()(0)/8); ++x){
     529             :         
     530           0 :         iblc[0]=x-subout.shape()(0)/2;
     531           0 :         itrc[0]=x-subout.shape()(0)/2;
     532           0 :         Slicer isl(iblc, itrc, Slicer::endIsLast);
     533           0 :         iblc[0]=x;
     534           0 :         subout.putSlice(subout.getSlice(isl), iblc);
     535           0 :       }
     536           0 :       for(Int x=subout.shape()(0)/2+1; x <(5*subout.shape()(0)/8); ++x){
     537             :         
     538           0 :         iblc[0]=x;
     539           0 :         itrc[0]=x;
     540           0 :         Slicer isl(iblc, itrc, Slicer::endIsLast);
     541           0 :         iblc[0]=subout.shape()(0)-x;
     542           0 :         subout.putSlice(subout.getSlice(isl), iblc);
     543           0 :         if(x==(subout.shape()(0)-1)){
     544           0 :           iblc[0]=0;
     545           0 :           subout.putSlice(subout.getSlice(isl), iblc);
     546             :         }
     547           0 :       }
     548             :       //End of FFT's
     549             :       //tim.show("After apply2 ");
     550             :        /////////////////////
     551             :       //wtime0=omp_get_wtime();
     552             :       ///////////////////
     553           0 :       TempImage<Complex> twoDPB2(TiledShape(pbShape, IPosition(4, pbShape(0), pbShape(1), 1, 1)), coords, memtobeused/10.0);
     554             :       
     555           0 :       IPosition blcout(4, 0, 0, 0, nBeamChans-1);
     556           0 :       IPosition trcout(4, pbShape(0)-1, pbShape(1)-1, 0,nBeamChans-1);
     557           0 :       Slicer outsl(blcout, trcout, Slicer::endIsLast);
     558             :       
     559           0 :       IPosition blcin(4, convSize_p/2-pbShape(0)/2, convSize_p/2-pbShape(1)/2, 0, 0);
     560           0 :       IPosition trcin(4, convSize_p/2+pbShape(0)/2-1, convSize_p/2+pbShape(1)/2-1, 0, 0);
     561           0 :       Slicer insl(blcin, trcin, Slicer::endIsLast);
     562             :       {
     563           0 :         SubImage<Complex> subtwoDPB(twoDPB, outsl, true);
     564           0 :         SubImage<Complex> intwoDPB(subim, insl, false);
     565           0 :         subtwoDPB.copyData(intwoDPB);
     566           0 :       }
     567             :       {
     568           0 :         SubImage<Complex> subtwoDPB2(twoDPB2, outsl, true);
     569           0 :         SubImage<Complex> intwoDPB2(subout, insl, false);
     570           0 :         subtwoDPB2.copyData(intwoDPB2);
     571           0 :       }
     572             :       
     573           0 :       if(nBeamChans > 0){
     574           0 :         blcin=IPosition(4,0,0,0, nBeamChans-1);
     575           0 :         trcin=IPosition(4, pbShape(0)-1, pbShape(1)-1, 0, nBeamChans-1);
     576           0 :         Slicer slin(blcin, trcin, Slicer::endIsLast);
     577           0 :         SubImage<Complex> origPB(twoDPB, slin, false);
     578           0 :         IPosition elshape= origPB.shape();
     579           0 :         Matrix<Complex> i1=origPB.get(true);
     580           0 :         SubImage<Complex> origPB2(twoDPB2, slin, false);
     581           0 :         Matrix<Complex> i2=origPB2.get(true);
     582           0 :         Int cenX=i1.shape()(0)/2;
     583           0 :         Int cenY=i1.shape()(1)/2;
     584             :         
     585             :            
     586           0 :         for (Int kk=0; kk < nBeamChans; ++kk){
     587           0 :           Double fratio=beamFreqs(kk)/beamFreqs(nBeamChans-1);
     588             :           //cerr << "fratio " << fratio << endl;
     589           0 :           Float convRatio=convSamp/convSampling;
     590           0 :           blcin[3]=kk;
     591           0 :           trcin[3]=kk;
     592             :           //Slicer slout(blcin, trcin, Slicer::endIsLast);
     593           0 :           Matrix<Complex> o1(i1.shape(), Complex(0.0));
     594           0 :           Matrix<Complex> o2(i2.shape(), Complex(0.0));
     595           0 :           for (Int yy=0;  yy < i1.shape()(1); ++yy){
     596             :             //Int nyy= (Double(yy-cenY)*fratio) + cenY; 
     597           0 :             Double nyy= (Double((yy-cenY)*convRatio)/fratio) + cenY;
     598           0 :             Double cyy=ceil(nyy);
     599           0 :             Double fyy= floor(nyy);
     600           0 :             Int iy=nyy > fyy+0.5 ? Int(cyy) : Int(fyy); 
     601           0 :             if(cyy <2*cenY && fyy >=0.0)
     602           0 :               for(Int xx=0; xx < i1.shape()(0); ++ xx){
     603             :                 //Int nxx= Int(Double(xx-cenX)*fratio) + cenX; 
     604           0 :                 Double nxx= Int(Double((xx-cenX)*convRatio)/fratio) + cenX;
     605           0 :                 Double cxx=ceil(nxx);
     606           0 :                 Double fxx= floor(nxx);
     607           0 :                 Int ix=nxx > fxx+0.5 ? Int(cxx) : Int(fxx) ;
     608           0 :                 if(cxx < 2*cenX && fxx >=0.0 ){
     609             :                   //Double dist=sqrt((nxx-cxx)*(nxx-cxx)+(nyy-cyy)*(nyy-cyy))/sqrt(2.0);
     610             :                   //o1(xx, yy)=float(1-dist)*i1(fxx, fyy)+ dist*i1(cxx,cyy);
     611           0 :                   o1(xx, yy)=i1( ix, iy);
     612             :                   //o2(xx, yy)=i2(nxx, nyy);
     613             :                   //o2(xx, yy)=float(1-dist)*i2(fxx, fyy)+ dist*i2(cxx,cyy);
     614           0 :                   o2(xx, yy)=i2(ix, iy);
     615             :                 }
     616             :               }
     617             :           }
     618           0 :           twoDPB.putSlice(o1.reform(elshape), blcin);
     619           0 :           twoDPB2.putSlice(o2.reform(elshape), blcin);
     620           0 :         }
     621             :         
     622           0 :       }
     623             :       //cerr << "make multfreq beam " << omp_get_wtime()-wtime0 << endl;
     624             :       /*
     625             :       {
     626             :         TempImage<Float> screen2(TiledShape(pbShape, IPosition(4, pbShape(0), pbShape(1), 1, 1)), coords, memtobeused);
     627             :           //    Matrix<Float> screenoo(convSize_p, convSize_p);
     628             :           //screenoo.set(1.0);
     629             :           //screen2.putSlice(screenoo,start);
     630             :           //screen2.set(1.0);
     631             :           for (Int k=0; k < nBeamChans; ++k){
     632             :             blcin[3]=k;
     633             :             trcin[3]=k;
     634             :             Slicer slin(blcin, trcin, Slicer::endIsLast);
     635             :             SubImage<Float> subim(screen2, slin, true);
     636             :             SubImage<Complex> subout(twoDPB2, slin, true);
     637             :             subim.set(1.0);
     638             :             //twoDPB.putSlice(screen, start);
     639             :             sj_p->applySquare(subim, subim, vb, 0); 
     640             :             //// LatticeExpr<Complex> le(subim);
     641             :             //// subout.copyData(le);
     642             :             ///// LatticeFFT::cfft2d(subout);
     643             :            
     644             :             LatticeFFT::rcfft(subout, subim, true, false);
     645             :             IPosition iblc(4, 0, 3*subout.shape()(1)/8, 0, 0);
     646             :             IPosition itrc(4, 0, 5*subout.shape()(1)/8, 0, 0);
     647             :             for(Int x=subout.shape()(0)/2; x <(5*subout.shape()(0)/8); ++x){
     648             :               
     649             :               iblc[0]=x-subout.shape()(0)/2;
     650             :               itrc[0]=x-subout.shape()(0)/2;
     651             :               Slicer isl(iblc, itrc, Slicer::endIsLast);
     652             :               iblc[0]=x;
     653             :               subout.putSlice(subout.getSlice(isl), iblc);
     654             :             }
     655             :             for(Int x=subout.shape()(0)/2+1; x <(5*subout.shape()(0)/8); ++x){
     656             :               
     657             :               iblc[0]=x;
     658             :               itrc[0]=x;
     659             :               Slicer isl(iblc, itrc, Slicer::endIsLast);
     660             :               iblc[0]=subout.shape()(0)-x;
     661             :               subout.putSlice(subout.getSlice(isl), iblc);
     662             :               if(x==(subout.shape()(0)-1)){
     663             :                 iblc[0]=0;
     664             :                 subout.putSlice(subout.getSlice(isl), iblc);
     665             :               }
     666             :             }
     667             :             
     668             :           }
     669             :       
     670             :           //sj_p->applySquare(screen2, screen2, vb, 0);
     671             :           //LatticeExpr<Complex> le(screen2);
     672             :           //twoDPB2.copyData(le);
     673             :       }
     674             :       
     675             :       */ 
     676             :       
     677             :       /*
     678             :       if(0) {
     679             :         CoordinateSystem ftCoords(coords);
     680             :         //directionIndex=ftCoords.findCoordinate(Coordinate::DIRECTION);
     681             :         //AlwaysAssert(directionIndex>=0, AipsError);
     682             :         dc=coords.directionCoordinate(directionIndex);
     683             :         Vector<Bool> axes(2); axes(0)=true;axes(1)=true;
     684             :         Vector<Int> shape(2); shape(0)=convSize_p;shape(1)=convSize_p;
     685             :         Coordinate* ftdc=dc.makeFourierCoordinate(axes,shape);
     686             :         //ftCoords.replaceCoordinate(*ftdc, directionIndex);
     687             :         delete ftdc; ftdc=0;
     688             :         ostringstream os1;
     689             :         os1 << "Screen_" << vb.fieldId() ;
     690             :         PagedImage<Complex> thisScreen(twoDPB2.shape(), ftCoords, String(os1));
     691             :         //LatticeExpr<Float> le(abs(twoDPB2));
     692             :         thisScreen.copyData(twoDPB2);
     693             :         LatticeFFT::cfft2d(thisScreen, false);
     694             :         PagedImage<Complex> thisScreen0(twoDPB.shape(), ftCoords, String("PB_")+String(os1));
     695             :         thisScreen0.copyData(twoDPB);
     696             :         LatticeFFT::cfft2d(thisScreen0, false);
     697             :       }
     698             :       */
     699             :       /* 
     700             :       // Now FFT and get the result back
     701             :       //LatticeFFT::cfft2d(twoDPB);
     702             :       //LatticeFFT::cfft2d(twoDPB2);
     703             :       
     704             :       // Write out FT of screen as an image
     705             :       if(0) {
     706             :         CoordinateSystem ftCoords(coords);
     707             :         directionIndex=ftCoords.findCoordinate(Coordinate::DIRECTION);
     708             :         AlwaysAssert(directionIndex>=0, AipsError);
     709             :         dc=coords.directionCoordinate(directionIndex);
     710             :         Vector<Bool> axes(2); axes(0)=true;axes(1)=true;
     711             :         Vector<Int> shape(2); shape(0)=convSize_p;shape(1)=convSize_p;
     712             :         Coordinate* ftdc=dc.makeFourierCoordinate(axes,shape);
     713             :         ftCoords.replaceCoordinate(*ftdc, directionIndex);
     714             :         delete ftdc; ftdc=0;
     715             :         ostringstream os1;
     716             :         os1 << "FTScreen_" << vb.fieldId() ;
     717             :         PagedImage<Float> thisScreen(pbShape, ftCoords, String(os1));
     718             :         LatticeExpr<Float> le(abs(twoDPB2));
     719             :         thisScreen.copyData(le);
     720             :       }
     721             :       */
     722             :       //cerr << "twoDPB shape " << twoDPB.shape() << " slice shape " << IPosition(4, tempConvSize, tempConvSize, 1, 1) << endl;
     723           0 :       convFunc_p=twoDPB.getSlice(IPosition(4,0,0,0,0), IPosition(4, tempConvSize, tempConvSize, 1, 1), true);
     724           0 :       weightConvFunc_p.resize();
     725           0 :       weightConvFunc_p=twoDPB2.getSlice(IPosition(4,0,0,0,nBeamChans-1), IPosition(4, tempConvSize, tempConvSize, 1, 1), true);
     726             :       //convFunc/=max(abs(convFunc));
     727           0 :       Float maxAbsConvFunc=max(amplitude(weightConvFunc_p));
     728             :       
     729           0 :       Float minAbsConvFunc=min(amplitude(weightConvFunc_p));
     730             :       //cerr << "min max " << minAbsConvFunc << "  " <<  maxAbsConvFunc << endl;
     731           0 :       convSupport_p=-1;
     732           0 :       Bool found=false;
     733             :       //Bool found2=true;
     734             :       //Int trial2=0;
     735           0 :       Int trial=0;
     736           0 :       for (trial=tempConvSize/2-2;trial>0;trial--) {
     737             :         //Searching down a diagonal
     738           0 :         if(abs(weightConvFunc_p(tempConvSize/2-trial, tempConvSize/2-trial)) >  (1.0e-3*maxAbsConvFunc)) {
     739           0 :           found=true;
     740           0 :           trial=Int(sqrt(2.0*Float(trial*trial)));
     741           0 :           break;
     742             :         }
     743             :       }
     744           0 :       if(!found){
     745           0 :         if((maxAbsConvFunc-minAbsConvFunc) > (1.0e-3*maxAbsConvFunc)) 
     746           0 :           found=true;
     747             :         // if it drops by more than 2 magnitudes per pixel
     748           0 :         trial=( tempConvSize > (10*convSampling)) ? 5*convSampling : (tempConvSize/2 - 4*convSampling);
     749             :       }
     750             : 
     751           0 :       if(trial < 5*convSampling) 
     752           0 :         trial=( tempConvSize > (10*convSampling)) ? 5*convSampling : (tempConvSize/2 - 4*convSampling);
     753             :       
     754           0 :       if(found) {
     755           0 :         convSupport_p=Int(0.5+Float(trial)/Float(convSampling))+1;
     756             :       }
     757             :       else {
     758             :         os << "Convolution function is misbehaved - support seems to be zero\n"
     759             :            << "Reasons can be: \n(1)The image definition not covering one or more of the pointings selected\n"
     760             :            << "(2) No unflagged data in a given pointing\n"
     761             :            << "(3) The entries in the POINTING subtable do not match the field being imaged."
     762             :            << "Please check, and try again with an empty POINTING subtable.)\n"
     763           0 :            << LogIO::EXCEPTION;
     764             :       }
     765             : 
     766             :       // Normalize such that plane 0 sums to 1 (when jumping in
     767             :       // steps of convSampling)
     768             :       
     769           0 :       Double pbSum=0.0;
     770             :       //Double pbSum2=0.0;
     771             :    
     772             : 
     773           0 :       for (Int iy=-convSupport_p;iy<=convSupport_p;iy++) {
     774           0 :         for (Int ix=-convSupport_p;ix<=convSupport_p;ix++) {
     775           0 :           Complex val=convFunc_p(ix*convSampling+tempConvSize/2,
     776           0 :                                  iy*convSampling+tempConvSize/2);
     777           0 :           pbSum+=real(val);
     778             :       
     779             :           //pbSum+=sqrt(real(val)*real(val)+ imag(val)*imag(val));
     780             :         }
     781             :       }
     782             :       
     783             :  
     784             :       //pbSum=sum(amplitude(convFunc_p))/Double(convSampling)/Double(convSampling);
     785             : 
     786           0 :       if(pbSum>0.0) {
     787           0 :         convFunc_p*=Complex(1.0/pbSum,0.0);
     788             :       }
     789             :       else {
     790             :         os << "Convolution function integral is not positive"
     791           0 :            << LogIO::EXCEPTION;
     792             :       }
     793             :       
     794             :       //##########################################
     795             :       os << "Convolution support = " << convSupport_p
     796             :          << " pixels in Fourier plane"
     797           0 :          << LogIO::POST;
     798             :       
     799           0 :       convSupportBlock_p.resize(actualConvIndex_p+1);
     800           0 :       convSizes_p.resize(actualConvIndex_p+1);
     801             :       //Only one beam for now...but later this should be able to
     802             :       // take all the beams for the different antennas.
     803           0 :       convSupportBlock_p[actualConvIndex_p]=new Vector<Int>(1);
     804           0 :       convSizes_p[actualConvIndex_p]=new Vector<Int> (1);
     805           0 :       (*(convSupportBlock_p[actualConvIndex_p]))[0]=convSupport_p;
     806           0 :       convFunctions_p.resize(actualConvIndex_p+1);
     807           0 :       convWeights_p.resize(actualConvIndex_p+1);
     808           0 :       convFunctions_p[actualConvIndex_p]= new Array<Complex>();
     809           0 :       convWeights_p[actualConvIndex_p]= new Array<Complex>();
     810           0 :       Int newConvSize=2*(convSupport_p+2)*convSampling;
     811             :       //NEED to chop this right ...and in the centre
     812           0 :       if(newConvSize >=tempConvSize)
     813           0 :           newConvSize=tempConvSize;
     814             : 
     815           0 :       IPosition blc(4, (tempConvSize/2)-(newConvSize/2),
     816           0 :                   (tempConvSize/2)-(newConvSize/2), 0, 0);
     817           0 :       IPosition trc(4, (tempConvSize/2)+(newConvSize/2-1),
     818           0 :                       (tempConvSize/2)+(newConvSize/2-1), 0, nBeamChans-1);
     819           0 :       convFunctions_p[actualConvIndex_p]->resize(IPosition(5, newConvSize, newConvSize, 1, nBeamChans,1));
     820             :       //cerr << "convFunc shape " << (convFunctions_p[actualConvIndex_p])->shape() << 
     821             :       //"  " << " twoDPB shape " <<twoDPB.get(false)(blc,trc).shape() << endl;
     822           0 :       convFunctions_p[actualConvIndex_p]->copyMatchingPart(twoDPB.get(false)(blc,trc));//*Complex(1.0/pbSum,0.0));
     823           0 :       convSize_p=newConvSize;
     824           0 :       convWeights_p[actualConvIndex_p]->resize(IPosition(5, newConvSize, newConvSize, 1, nBeamChans,1));
     825           0 :       convWeights_p[actualConvIndex_p]->copyMatchingPart(twoDPB2.get(false)(blc,trc));//*Complex(1.0/pbSum2,0.0));
     826           0 :       blc.resize(5, false);
     827           0 :       trc.resize(5,false);
     828           0 :       blc=IPosition(5, 0, 0, 0, 0,0);
     829           0 :       trc=IPosition(5, newConvSize-1, newConvSize-1, 0, 0,0);
     830           0 :       for (Int bc=0; bc< nBeamChans; ++bc){
     831           0 :         blc[3]=bc;
     832           0 :         trc[3]=bc;
     833           0 :         pbSum=real(sum((*convFunctions_p[actualConvIndex_p])(blc,trc)))/Double(convSampling)/Double(convSampling);
     834           0 :         (*convFunctions_p[actualConvIndex_p])(blc,trc)=         (*convFunctions_p[actualConvIndex_p])(blc,trc)/pbSum;
     835           0 :         (*convWeights_p[actualConvIndex_p])(blc,trc)=   (*convWeights_p[actualConvIndex_p])(blc,trc)/pbSum;
     836             :       }
     837             : 
     838             : 
     839             : 
     840           0 :       convFunc_p.resize();//break any reference
     841           0 :       (*convSizes_p[actualConvIndex_p])[0]=convSize_p;
     842           0 :       doneMainConv_p[actualConvIndex_p]=true;
     843             :       
     844           0 :     }
     845             :     else{
     846           0 :       convSize_p=(*convSizes_p[actualConvIndex_p])[0];
     847             : 
     848             :     }
     849             : 
     850             :     //Apply the shift phase gradient
     851           0 :     convFunc.resize();
     852           0 :     weightConvFunc.resize();
     853           0 :     convFunc.assign(*(convFunctions_p[actualConvIndex_p]));
     854           0 :     weightConvFunc.assign(*(convWeights_p[actualConvIndex_p]));
     855             :     Bool copyconv, copywgt;
     856           0 :     Complex *cv=convFunc.getStorage(copyconv);
     857           0 :     Complex *wcv=weightConvFunc.getStorage(copywgt);
     858             :     //cerr << "Field " << vb.fieldId() << " spw " << vb.spectralWindow() << " phase grad: " << pixFieldDir << endl;
     859             :    
     860           0 :     for (Int nc=0; nc < nBeamChans; ++nc){ 
     861           0 :         Int planeoffset=nc*convSize_p*convSize_p;
     862           0 :         for (Int iy=0;iy<convSize_p;iy++) {
     863             :                 Double cy, sy;
     864             :                 Int offset;
     865             :                
     866           0 :                 SINCOS(Double(iy-convSize_p/2)*pixFieldDir(1), sy, cy);
     867           0 :                 Complex phy(cy,sy) ;
     868           0 :                 offset = iy*convSize_p+planeoffset;
     869           0 :                 for (Int ix=0;ix<convSize_p;ix++) {
     870             :                         Double cx, sx;
     871           0 :                         SINCOS(Double(ix-convSize_p/2)*pixFieldDir(0), sx, cx);
     872           0 :                         Complex phx(cx,sx) ;
     873           0 :                         cv[ix+offset]= cv[ix+offset]*phx*phy;
     874           0 :                         wcv[ix+offset]= wcv[ix+offset]*phx*phy;
     875             : 
     876             :                 }
     877             :         }
     878             :     }
     879           0 :     convFunc.putStorage(cv, copyconv);
     880           0 :     weightConvFunc.putStorage(wcv, copywgt);
     881           0 :     convsize.resize();
     882           0 :     convsize=*(convSizes_p[actualConvIndex_p]);
     883           0 :     convSupport.resize();
     884           0 :     convSupport=(*(convSupportBlock_p[actualConvIndex_p]));
     885             :     
     886             :     
     887           0 :   }
     888             : 
     889           0 :   void SimplePBConvFunc::setSkyJones(SkyJones* sj){
     890           0 :     sj_p=sj;
     891           0 :   }
     892             : 
     893           0 :   void SimplePBConvFunc::findUsefulChannels(Vector<Int>& chanMap, Vector<Double>& chanFreqs,  const VisBuffer& vb, const Vector<Double>& freq){
     894           0 :     chanMap.resize(freq.nelements());
     895           0 :     Vector<Double> localfreq=vb.frequency();
     896           0 :     Double minfreq=min(freq);
     897             :     
     898           0 :     Double origwidth=freq.nelements()==1 ? 1e12 : (max(freq)-min(freq))/(freq.nelements()-1);
     899             :     ///Fractional bandwidth which will trigger mutiple PB in one spw
     900           0 :     Double tol=(max(freq))*0.5/100;
     901           0 :     Int nchan=Int(lround((max(freq)-min(freq))/tol));
     902             :    
     903             :     //cerr  << "TOLERA " << tol << " nchan " << nchan << " vb.nchan " << vb.nChannel() << endl;
     904             :     //Number of beams that matters are the ones in the data
     905           0 :     if(nchan > vb.nChannel())
     906           0 :       nchan=vb.nChannel();
     907             : 
     908           0 :     if(tol < origwidth) tol=origwidth;
     909           0 :     chanFreqs.resize();
     910           0 :     if(nchan >= (Int)(freq.nelements()-1)) { indgen(chanMap); chanFreqs=freq; return;}
     911           0 :     if((nchan==0) || (freq.nelements()==1)) { chanFreqs=Vector<Double>(1, freq[0]);chanMap.set(0); return;}
     912             : 
     913             :     //readjust the tolerance...
     914           0 :     tol=(max(freq)-min(freq)+origwidth)/Double(nchan);
     915           0 :     chanFreqs.resize(nchan);
     916           0 :     for (Int k=0; k < nchan; ++k)
     917           0 :       chanFreqs[k]=minfreq-origwidth+tol/2.0+tol*Double(k);
     918           0 :     Int activechan=0;
     919           0 :     chanMap.set(-1);
     920           0 :     for (uInt k=0; k < chanMap.nelements(); ++k){
     921             :      
     922           0 :       while((activechan< nchan) && Float(fabs(freq[k]-chanFreqs[activechan])) > Float(tol/2.0)){
     923             :         //              cerr << "k " << k << " atcivechan " << activechan << " comparison " 
     924             :         //     << freq[k] << "    " << chanFreqs[activechan]  << endl;        
     925           0 :         ++activechan;
     926             :       }
     927           0 :       if(activechan != nchan)
     928           0 :         chanMap[k]=activechan;
     929             :       //////////////////
     930             :       //if(chanMap[k] < 0)
     931             :       //cerr << "freq diffs " << freq[k]-chanFreqs << "  TOL " << tol/2.0 << endl;
     932             : 
     933             :       ///////////////////////////
     934           0 :       activechan=0;
     935             :     }
     936             : 
     937           0 :     return;
     938           0 :   }
     939             : 
     940             : 
     941           0 :   Bool SimplePBConvFunc::checkPBOfField(const VisBuffer& vb){
     942             :     //Int fieldid=vb.fieldId();
     943           0 :     String msid=vb.msName(true);
     944             :     /*
     945             :      if(convFunctionMap_p.ndefined() > 0){
     946             :       if (((fluxScale_p.shape()[3] != nchan_p) || (fluxScale_p.shape()[2] != npol_p)) && calcFluxScale_p){
     947             :         convFunctionMap_p.clear();
     948             :       }
     949             :     }
     950             :     // if you rename the ms might be a problem
     951             :     String mapid=msid+String("_")+String::toString(fieldid);
     952             :     if(convFunctionMap_p.ndefined() == 0){
     953             :       convFunctionMap_p.define(mapid, 0);    
     954             :       actualConvIndex_p=0;
     955             :       if(calcFluxScale_p){
     956             :         // 0ne channel only is needed to keep track of pb coverage
     957             :         if(fluxScale_p.shape().nelements()==0){
     958             :           fluxScale_p=TempImage<Float>(IPosition(4,nx_p,ny_p,npol_p,1), csys_p);
     959             :           fluxScale_p.set(0.0);
     960             :         }
     961             :         filledFluxScale_p=false;
     962             :       }
     963             :       return false;
     964             :     }
     965             :     
     966             :     if(!convFunctionMap_p.isDefined(mapid)){
     967             :       actualConvIndex_p=convFunctionMap_p.ndefined();
     968             :       convFunctionMap_p.define(mapid, actualConvIndex_p);
     969             :       return false;
     970             :     }
     971             :     else{
     972             :       actualConvIndex_p=convFunctionMap_p(mapid);
     973             :       convFunc_p.resize(); // break any reference
     974             :       weightConvFunc_p.resize(); 
     975             :       //Here we will need to use the right xyPlane for different PA range.
     976             :       //convFunc_p.reference(convFunctions_p[actualConvIndex_p]->xyPlane(0));
     977             :       //weightConvFunc_p.reference(convWeights_p[actualConvIndex_p]->xyPlane(0));
     978             :       //Again this for one time of antenna only later should be fixed for all 
     979             :       // antennas independently
     980             :       convSupport_p=(*convSupportBlock_p[actualConvIndex_p])[0];
     981             :       convSize_p=(*convSizes_p[actualConvIndex_p])[0];
     982             : 
     983             :   }
     984             : */
     985             :  
     986           0 :  return true;
     987             : 
     988             : 
     989             : 
     990           0 :   }
     991           0 :   ImageInterface<Float>&  SimplePBConvFunc::getFluxScaleImage(){
     992             : 
     993           0 :     if(!calcFluxScale_p)
     994           0 :       throw(AipsError("Programmer error: Cannot get flux scale"));
     995           0 :     if(!filledFluxScale_p){
     996           0 :       IPosition blc=fluxScale_p.shape();
     997           0 :       IPosition trc=fluxScale_p.shape();
     998           0 :       blc(0)=0; blc(1)=0; trc(0)=nx_p-1; trc(1)=ny_p-1;
     999             :       
    1000           0 :       for (Int j=0; j < fluxScale_p.shape()(2); ++j){
    1001           0 :         for (Int k=0; k < fluxScale_p.shape()(3) ; ++k){
    1002             :           
    1003           0 :           blc(2)=j; trc(2)=j;
    1004           0 :           blc(3)=k; trc(3)=k;
    1005           0 :           Slicer sl(blc, trc, Slicer::endIsLast);
    1006           0 :           SubImage<Float> fscalesub(fluxScale_p, sl, true);
    1007             :           Float planeMax;
    1008           0 :           LatticeExprNode LEN = max( fscalesub );
    1009           0 :           planeMax =  LEN.getFloat();
    1010           0 :           if(planeMax !=0){
    1011           0 :             fscalesub.copyData( (LatticeExpr<Float>) (fscalesub/planeMax));
    1012             :             
    1013             :           }
    1014           0 :         }
    1015             :       }
    1016             :       
    1017           0 :       filledFluxScale_p=true;
    1018           0 :     }
    1019             :       
    1020             : 
    1021           0 :     return fluxScale_p;
    1022             :     
    1023             :   }
    1024             : 
    1025             : 
    1026           0 :   Bool SimplePBConvFunc::toRecord(RecordInterface& rec){
    1027           0 :     Int numConv=convFunctions_p.nelements();
    1028             :     // not saving the protected variables as they are generated by
    1029             :     // the first  call to storeImageParams 
    1030             :     try{
    1031           0 :       rec.define("name", "SimplePBConvFunc");
    1032           0 :       rec.define("numconv", numConv);
    1033             :       //cerr << "num of conv " << numConv << "  " << convFunctionMap_p.ndefined() << "  " <<convFunctions_p.nelements() << endl;
    1034           0 :       std::map<String, Int>::iterator it=vbConvIndex_p.begin();
    1035           0 :       for (Int k=0; k < numConv; ++k){
    1036           0 :         rec.define("convfunctions"+String::toString(k), *(convFunctions_p[k]));
    1037           0 :         rec.define("convweights"+String::toString(k), *(convWeights_p[k]));
    1038           0 :         rec.define("convsizes"+String::toString(k), *(convSizes_p[k]));
    1039           0 :         rec.define("convsupportblock"+String::toString(k), *(convSupportBlock_p[k]));
    1040             :         //cerr << "k " << k << " key " << convFunctionMap_p.getKey(k) << " val " << convFunctionMap_p.getVal(k) << endl;
    1041           0 :         rec.define(String("key")+String::toString(k), it->first);
    1042           0 :         rec.define(String("val")+String::toString(k), it->second);
    1043           0 :         it++;
    1044             :       }
    1045           0 :       rec.define("pbclass", Int(pbClass_p));
    1046           0 :       rec.define("actualconvindex",  actualConvIndex_p);
    1047           0 :       rec.define("donemainconv", doneMainConv_p);
    1048             :       //The following is not needed ..can be regenerated
    1049             :       //rec.define("pointingpix", pointingPix_p);
    1050             :     }
    1051           0 :     catch(AipsError &x) {
    1052           0 :       return false;
    1053           0 :     }
    1054           0 :     return true;
    1055             :   }
    1056             : 
    1057           0 :   Bool SimplePBConvFunc::fromRecord(String& err, const RecordInterface& rec, Bool calcFluxneeded){
    1058           0 :      Int numConv=0;
    1059             :      //make sure storeImageParams is triggered
    1060           0 :      nchan_p=0;
    1061             :      
    1062             :      try{
    1063           0 :        if(!rec.isDefined("name") || rec.asString("name") != "SimplePBConvFunc"){
    1064           0 :          throw(AipsError("Wrong record to recover SimplePBConvFunc from"));
    1065             :         }
    1066           0 :        rec.get("numconv", numConv);
    1067           0 :        convFunctions_p.resize(numConv, true, false);
    1068           0 :        convSupportBlock_p.resize(numConv, true, false);
    1069           0 :        convWeights_p.resize(numConv, true, false);
    1070           0 :        convSizes_p.resize(numConv, true, false);
    1071           0 :        convFunctionMap_p.clear( );
    1072           0 :        vbConvIndex_p.erase(vbConvIndex_p.begin(), vbConvIndex_p.end());
    1073           0 :        for (Int k=0; k < numConv; ++k){
    1074           0 :          convFunctions_p[k]=new Array<Complex>();
    1075           0 :          convWeights_p[k]=new Array<Complex>();
    1076           0 :          convSizes_p[k]=new Vector<Int>();
    1077           0 :          convSupportBlock_p[k]=new Vector<Int>();
    1078           0 :          rec.get("convfunctions"+String::toString(k), *(convFunctions_p[k]));
    1079           0 :          rec.get("convsupportblock"+String::toString(k), *(convSupportBlock_p[k]));
    1080           0 :          rec.get("convweights"+String::toString(k), *(convWeights_p[k]));
    1081           0 :          rec.get("convsizes"+String::toString(k), *(convSizes_p[k]));
    1082           0 :          String key;
    1083             :          Int val;
    1084           0 :          rec.get(String("key")+String::toString(k), key);
    1085           0 :          rec.get(String("val")+String::toString(k), val);
    1086           0 :          vbConvIndex_p[key]=val;
    1087           0 :          ant1PointVal_p.clear();
    1088           0 :          ant1PointingCache_p.resize();
    1089             :          //convFunctionMap_p.define(key,val);
    1090           0 :        }
    1091           0 :        pbClass_p=static_cast<PBMathInterface::PBClass>(rec.asInt("pbclass"));
    1092           0 :        rec.get("actualconvindex",  actualConvIndex_p);
    1093           0 :        pointingPix_p.resize();
    1094             :        //rec.get("pointingpix", pointingPix_p);
    1095           0 :        calcFluxScale_p=calcFluxneeded;
    1096             : 
    1097             :      }
    1098           0 :      catch(AipsError & x) {
    1099           0 :        err=x.getMesg();
    1100           0 :        return false;
    1101           0 :      }
    1102           0 :      return true;
    1103             :      
    1104             :   }
    1105           0 :   void SimplePBConvFunc::addPBToFlux(const VisBuffer& vb){
    1106             :     
    1107           0 :     if(calcFluxScale_p){
    1108           0 :       if(fluxScale_p.shape().nelements()==0)
    1109             :         { //cerr << "nx_p " << nx_p << endl;
    1110           0 :         calcFluxScale_p=False;
    1111           0 :         return;
    1112             :       }
    1113           0 :       Vector<Int> pixdepoint(2, -100000);
    1114           0 :       convertArray(pixdepoint, thePix_p);
    1115           0 :       if((pixdepoint(0) >=0) && (pixdepoint(0) < pointingPix_p.shape()[0]) && (pixdepoint(1) >=0) && (pixdepoint(1) < pointingPix_p.shape()[1])  && !pointingPix_p(pixdepoint(0), pixdepoint(1))){
    1116           0 :          TempImage<Float> thispb(fluxScale_p.shape(), fluxScale_p.coordinates());
    1117           0 :          thispb.set(1.0);
    1118           0 :          sj_p->applySquare(thispb, thispb, vb, 0);
    1119           0 :          LatticeExpr<Float> le(fluxScale_p+thispb);
    1120           0 :          fluxScale_p.copyData(le);
    1121           0 :          pointingPix_p(pixdepoint(0), pixdepoint(1))=true;
    1122             :          //LatticeExprNode LEN = max(fluxScale_p);
    1123             :          //Float maxsca=LEN.getFloat();
    1124             :          //Tempporary fix when cubesky is chunking...do not add on 
    1125             :          //already defined position
    1126             :          //if(maxsca > 1.98){
    1127             :          //  cerr << "avoiding subtract " << endl;
    1128             :         //fluxScale_p.copyData(LatticeExpr<Float>(fluxScale_p-thispb));
    1129             : 
    1130             :          //}      
    1131             :       /*
    1132             :       if(0) {
    1133             :         ostringstream os1;
    1134             :         os1 << "SINGLE_" << vb.fieldId() ;
    1135             :         PagedImage<Float> thisScreen(fluxScale_p.shape(), fluxScale_p.coordinates(), String(os1));
    1136             :         thisScreen.copyData(thispb);
    1137             :         ostringstream os2;
    1138             :         os2 << "ALL_" << vb.fieldId() ;
    1139             :         PagedImage<Float> thisScreen2(fluxScale_p.shape(), fluxScale_p.coordinates(), String(os2));
    1140             :         thisScreen2.copyData(fluxScale_p);
    1141             :       }
    1142             :       */
    1143           0 :        }
    1144           0 :     }
    1145             : 
    1146             :   }
    1147             : 
    1148           0 :   void SimplePBConvFunc::sliceFluxScale(Int npol) {
    1149           0 :      IPosition fshp=fluxScale_p.shape();
    1150           0 :      if (fshp(2)>npol){
    1151           0 :        npol_p=npol;
    1152             :        // use first npol planes...
    1153           0 :        IPosition blc(4,0,0,0,0);
    1154           0 :        IPosition trc(4,fluxScale_p.shape()(0)-1, fluxScale_p.shape()(1)-1,npol-1,fluxScale_p.shape()(3)-1);
    1155           0 :        Slicer sl=Slicer(blc, trc, Slicer::endIsLast);
    1156             :        //writeable if possible
    1157           0 :        SubImage<Float> fluxScaleSub = SubImage<Float> (fluxScale_p, sl, true);
    1158           0 :        fluxScale_p = TempImage<Float>(fluxScaleSub.shape(),fluxScaleSub.coordinates());
    1159           0 :        LatticeExpr<Float> le(fluxScaleSub);
    1160           0 :        fluxScale_p.copyData(le);
    1161           0 :      }
    1162           0 :   }
    1163             : 
    1164             : } //# NAMESPACE CASA - END
    1165             : 
    1166             : 
    1167             : 
    1168             : 
    1169             : 
    1170             : 
    1171             : 

Generated by: LCOV version 1.16