Line data Source code
1 : // -*- C++ -*-
2 : //# nPBWProjectFT.cc: Implementation of nPBWProjectFT class
3 : //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
4 : //# Associated Universities, Inc. Washington DC, USA.
5 : //#
6 : //# This library is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU Library General Public License as published by
8 : //# the Free Software Foundation; either version 2 of the License, or (at your
9 : //# option) any later version.
10 : //#
11 : //# This library is distributed in the hope that it will be useful, but WITHOUT
12 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
14 : //# License for more details.
15 : //#
16 : //# You should have received a copy of the GNU Library General Public License
17 : //# along with this library; if not, write to the Free Software Foundation,
18 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning AIPS++ should be addressed as follows:
21 : //# Internet email: casa-feedback@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : //#
27 : //# $Id$
28 :
29 : #include <msvis/MSVis/VisibilityIterator.h>
30 : #include <casacore/casa/Quanta/UnitMap.h>
31 : #include <casacore/casa/Quanta/MVTime.h>
32 : #include <casacore/casa/Quanta/UnitVal.h>
33 : #include <casacore/measures/Measures/Stokes.h>
34 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
35 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
36 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
37 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
38 : #include <casacore/coordinates/Coordinates/Projection.h>
39 : #include <casacore/ms/MeasurementSets/MSColumns.h>
40 : #include <casacore/ms/MeasurementSets/MSRange.h>
41 : #include <casacore/casa/BasicSL/Constants.h>
42 : #include <casacore/scimath/Mathematics/FFTServer.h>
43 : #include <synthesis/MeasurementComponents/nPBWProjectFT.h>
44 : #include <casacore/scimath/Mathematics/RigidVector.h>
45 : #include <msvis/MSVis/StokesVector.h>
46 : #include <synthesis/TransformMachines/StokesImageUtil.h>
47 : #include <msvis/MSVis/VisBuffer.h>
48 : #include <msvis/MSVis/VisSet.h>
49 : #include <casacore/images/Images/ImageInterface.h>
50 : #include <casacore/images/Images/ImageRegrid.h>
51 : #include <casacore/images/Images/PagedImage.h>
52 : #include <casacore/casa/Containers/Block.h>
53 : #include <casacore/casa/Containers/Record.h>
54 : #include <casacore/casa/Arrays/ArrayLogical.h>
55 : #include <casacore/casa/Arrays/ArrayMath.h>
56 : #include <casacore/casa/Arrays/Array.h>
57 : #include <casacore/casa/Arrays/MaskedArray.h>
58 : #include <casacore/casa/Arrays/Vector.h>
59 : #include <casacore/casa/Arrays/Slice.h>
60 : #include <casacore/casa/Arrays/Matrix.h>
61 : #include <casacore/casa/Arrays/Cube.h>
62 : #include <casacore/casa/Arrays/MatrixIter.h>
63 : #include <casacore/casa/BasicSL/String.h>
64 : #include <casacore/casa/Utilities/Assert.h>
65 : #include <casacore/casa/Exceptions/Error.h>
66 : #include <casacore/lattices/Lattices/ArrayLattice.h>
67 : #include <casacore/lattices/Lattices/SubLattice.h>
68 : #include <casacore/lattices/LRegions/LCBox.h>
69 : #include <casacore/lattices/LEL/LatticeExpr.h>
70 : #include <casacore/lattices/Lattices/LatticeCache.h>
71 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
72 : #include <casacore/lattices/Lattices/LatticeIterator.h>
73 : #include <casacore/lattices/Lattices/LatticeStepper.h>
74 : #include <casacore/casa/Utilities/CompositeNumber.h>
75 : #include <casacore/casa/OS/Timer.h>
76 : #include <casacore/casa/OS/HostInfo.h>
77 : #include <sstream>
78 :
79 : #include <synthesis/TransformMachines/VLACalcIlluminationConvFunc.h>
80 : #include <synthesis/TransformMachines/IlluminationConvFunc.h>
81 : #include <synthesis/MeasurementComponents/ExpCache.h>
82 : #include <synthesis/MeasurementComponents/CExp.h>
83 : #include <synthesis/TransformMachines/Utils.h>
84 : #include <synthesis/TransformMachines/SynthesisError.h>
85 : #include <casacore/measures/Measures/MEpoch.h>
86 : #include <casacore/measures/Measures/MeasTable.h>
87 : #include <casacore/scimath/Mathematics/MathFunc.h>
88 :
89 : #include <casacore/casa/System/ProgressMeter.h>
90 :
91 : #define CONVSIZE (1024*2)
92 : #define CONVWTSIZEFACTOR 1.0
93 : #define OVERSAMPLING 20
94 : #define THRESHOLD 1E-4
95 : #define USETABLES 0 // If equal to 1, use tabulated exp() and
96 : // complex exp() functions.
97 : #define MAXPOINTINGERROR 250.0 // Max. pointing error in arcsec used to
98 : // determine the resolution of the
99 : // tabulated exp() function.
100 : using namespace casacore;
101 : namespace casa { //# NAMESPACE CASA - BEGIN
102 :
103 : #define NEED_UNDERSCORES
104 : extern "C"
105 : {
106 : //
107 : // The Gridding Convolution Function (GCF) used by the underlying
108 : // gridder written in FORTRAN.
109 : //
110 : // The arguments must all be pointers and the value of the GCF at
111 : // the given (u,v) point is returned in the weight variable. Making
112 : // this a function which returns a complex value (namely the weight)
113 : // has problems when called in FORTRAN - I (SB) don't understand
114 : // why.
115 : //
116 : #if defined(NEED_UNDERSCORES)
117 : #define nwcppeij nwcppeij_
118 : #endif
119 : //
120 : //---------------------------------------------------------------
121 : //
122 : IlluminationConvFunc nwEij;
123 0 : void nwcppeij(Double *griduvw, Double *area,
124 : Double *raoff1, Double *decoff1,
125 : Double *raoff2, Double *decoff2,
126 : Int *doGrad,
127 : Complex *weight,
128 : Complex *dweight1,
129 : Complex *dweight2,
130 : Double *currentCFPA)
131 : {
132 0 : Complex w,d1,d2;
133 0 : nwEij.getValue(griduvw, raoff1, raoff2, decoff1, decoff2,
134 : area,doGrad,w,d1,d2,*currentCFPA);
135 0 : *weight = w;
136 0 : *dweight1 = d1;
137 0 : *dweight2 = d2;
138 0 : }
139 : }
140 : //
141 : //---------------------------------------------------------------
142 : //
143 : #define FUNC(a) (((a)))
144 0 : nPBWProjectFT::nPBWProjectFT(Int nWPlanes, Long icachesize,
145 : String& cfCacheDirName,
146 : Bool applyPointingOffset,
147 : Bool doPBCorr,
148 : Int itilesize,
149 : Float /*paSteps*/,
150 : Float pbLimit,
151 0 : Bool usezero)
152 0 : : FTMachine(), padding_p(1.0), nWPlanes_p(nWPlanes),
153 0 : imageCache(0), cachesize(icachesize), tilesize(itilesize),
154 0 : gridder(0), isTiled(false), arrayLattice( ), lattice( ),
155 0 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
156 0 : mspc(0), msac(0), pointingToImage(0), usezero_p(usezero),
157 0 : doPBCorrection(doPBCorr),
158 0 : Second("s"),Radian("rad"),Day("d"), noOfPASteps(0),
159 0 : pbNormalized(false),resetPBs(true), rotateAperture_p(true),
160 0 : cfCache(), paChangeDetector(), cfStokes(),Area(),
161 0 : avgPBSaved(false),avgPBReady(false)
162 : {
163 0 : epJ=NULL;
164 0 : convSize=0;
165 0 : tangentSpecified_p=false;
166 0 : lastIndex_p=0;
167 0 : paChangeDetector.reset();
168 0 : pbLimit_p=pbLimit;
169 : //
170 : // Get various parameters from the visibilities.
171 : //
172 0 : bandID_p=-1;
173 0 : if (applyPointingOffset) doPointing=1; else doPointing=0;
174 :
175 0 : convFuncCacheReady=false;
176 0 : PAIndex = -1;
177 0 : maxConvSupport=-1;
178 : //
179 : // Set up the Conv. Func. disk cache manager object.
180 : //
181 0 : cfCache.setCacheDir(cfCacheDirName.data());
182 0 : cfCache.initCache();
183 0 : convSampling=OVERSAMPLING;
184 0 : convSize=CONVSIZE;
185 0 : Long hostRAM = (HostInfo::memoryTotal(true)*1024); // In bytes
186 0 : hostRAM = hostRAM/(sizeof(Float)*2); // In complex pixels
187 0 : if (cachesize > hostRAM) cachesize=hostRAM;
188 0 : }
189 : //
190 : //---------------------------------------------------------------
191 : //
192 0 : nPBWProjectFT::nPBWProjectFT(const RecordInterface& stateRec)
193 0 : : FTMachine(),Second("s"),Radian("rad"),Day("d")
194 : {
195 : //
196 : // Construct from the input state record
197 : //
198 0 : String error;
199 :
200 0 : if (!fromRecord(error, stateRec)) {
201 0 : throw (AipsError("Failed to create nPBWProjectFT: " + error));
202 : };
203 0 : bandID_p = -1;
204 0 : PAIndex = -1;
205 0 : maxConvSupport=-1;
206 0 : convSampling=OVERSAMPLING;
207 0 : convSize=CONVSIZE;
208 0 : }
209 : //
210 : //----------------------------------------------------------------------
211 : //
212 0 : nPBWProjectFT::nPBWProjectFT(const nPBWProjectFT& other):FTMachine()
213 : {
214 0 : operator=(other);
215 0 : }
216 : //
217 : //---------------------------------------------------------------
218 : //
219 0 : nPBWProjectFT& nPBWProjectFT::operator=(const nPBWProjectFT& other)
220 : {
221 0 : if(this!=&other)
222 : {
223 : //Do the base parameters
224 0 : FTMachine::operator=(other);
225 :
226 :
227 0 : padding_p=other.padding_p;
228 :
229 0 : nWPlanes_p=other.nWPlanes_p;
230 0 : imageCache=other.imageCache;
231 0 : cachesize=other.cachesize;
232 0 : tilesize=other.tilesize;
233 0 : cfRefFreq_p = other.cfRefFreq_p;
234 0 : if(other.gridder==0) gridder=0;
235 : else
236 : {
237 0 : uvScale.resize();
238 0 : uvOffset.resize();
239 0 : uvScale=other.uvScale;
240 0 : uvOffset=other.uvOffset;
241 0 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
242 0 : uvScale, uvOffset,
243 0 : "SF");
244 : }
245 :
246 0 : isTiled=other.isTiled;
247 0 : lattice=0;
248 0 : arrayLattice=0;
249 :
250 0 : maxAbsData=other.maxAbsData;
251 0 : centerLoc=other.centerLoc;
252 0 : offsetLoc=other.offsetLoc;
253 0 : pointingToImage=other.pointingToImage;
254 0 : usezero_p=other.usezero_p;
255 :
256 :
257 0 : padding_p=other.padding_p;
258 0 : nWPlanes_p=other.nWPlanes_p;
259 0 : imageCache=other.imageCache;
260 0 : cachesize=other.cachesize;
261 0 : tilesize=other.tilesize;
262 0 : isTiled=other.isTiled;
263 0 : maxAbsData=other.maxAbsData;
264 0 : centerLoc=other.centerLoc;
265 0 : offsetLoc=other.offsetLoc;
266 0 : mspc=other.mspc;
267 0 : msac=other.msac;
268 0 : pointingToImage=other.pointingToImage;
269 0 : usezero_p=other.usezero_p;
270 0 : doPBCorrection = other.doPBCorrection;
271 0 : maxConvSupport= other.maxConvSupport;
272 :
273 0 : epJ=other.epJ;
274 0 : convSize=other.convSize;
275 0 : lastIndex_p=other.lastIndex_p;
276 0 : paChangeDetector=other.paChangeDetector;
277 0 : pbLimit_p=other.pbLimit_p;
278 : //
279 : // Get various parameters from the visibilities.
280 : //
281 0 : bandID_p = other.bandID_p;
282 0 : doPointing=other.doPointing;
283 :
284 0 : convFuncCacheReady=other.convFuncCacheReady;
285 0 : PAIndex = other.PAIndex;
286 0 : maxConvSupport=other.maxConvSupport;
287 : //
288 : // Set up the Conv. Func. disk cache manager object.
289 : //
290 0 : cfCache=other.cfCache;
291 0 : convSampling=other.convSampling;
292 0 : convSize=other.convSize;
293 0 : cachesize=other.cachesize;
294 :
295 0 : resetPBs=other.resetPBs;
296 0 : pbNormalized=other.pbNormalized;
297 0 : currentCFPA=other.currentCFPA;
298 0 : lastPAUsedForWtImg = other.lastPAUsedForWtImg;
299 0 : cfStokes=other.cfStokes;
300 0 : Area=other.Area;
301 0 : avgPB = other.avgPB;
302 0 : avgPBReady = other.avgPBReady;
303 : };
304 0 : return *this;
305 : };
306 : //
307 : //---------------------------------------------------------------------
308 : //
309 0 : int nPBWProjectFT::getVisParams(const VisBuffer& vb)
310 : {
311 : Double Freq;
312 0 : Vector<String> telescopeNames=vb.msColumns().observation().telescopeName().getColumn();
313 0 : for(uInt nt=0;nt<telescopeNames.nelements();nt++)
314 : {
315 0 : if ((telescopeNames(nt) != "VLA") && (telescopeNames(nt) != "EVLA"))
316 : {
317 0 : String mesg="pbwproject algorithm can handle only (E)VLA antennas for now.\n";
318 0 : mesg += "Erroneous telescope name = " + telescopeNames(nt) + ".";
319 0 : SynthesisError err(mesg);
320 0 : throw(err);
321 0 : }
322 0 : if (telescopeNames(nt) != telescopeNames(0))
323 : {
324 0 : String mesg="pbwproject algorithm does not (yet) handle inhomogeneous arrays!\n";
325 0 : mesg += "Not yet a \"priority\"!!";
326 0 : SynthesisError err(mesg);
327 0 : throw(err);
328 0 : }
329 : }
330 : // MSSpWindowColumns mssp(vb.msColumns().spectralWindow());
331 0 : Freq = vb.msColumns().spectralWindow().refFrequency()(0);
332 0 : Diameter_p=0;
333 0 : Nant_p = vb.msColumns().antenna().nrow();
334 0 : for (Int i=0; i < Nant_p; i++)
335 0 : if (!vb.msColumns().antenna().flagRow()(i))
336 : {
337 0 : Diameter_p = vb.msColumns().antenna().dishDiameter()(i);
338 0 : break;
339 : }
340 0 : if (Diameter_p == 0)
341 : {
342 0 : logIO() << LogOrigin("nPBWProjectFT", "getVisParams")
343 : << "No valid or finite sized antenna found in the antenna table. "
344 : << "Assuming diameter = 25m."
345 : << LogIO::WARN
346 0 : << LogIO::POST;
347 0 : Diameter_p=25.0;
348 : }
349 :
350 0 : Double Lambda=C::c/Freq;
351 0 : HPBW = Lambda/(Diameter_p*sqrt(log(2.0)));
352 0 : sigma = 1.0/(HPBW*HPBW);
353 0 : nwEij.setSigma(sigma);
354 0 : Int bandID = BeamCalc::Instance()->getBandID(Freq,telescopeNames(0),"");
355 : // Int bandID = getVLABandID(Freq,telescopeNames(0));
356 0 : return bandID;
357 0 : }
358 : //
359 : //----------------------------------------------------------------------
360 : //
361 0 : void nPBWProjectFT::init()
362 : {
363 0 : nx = image->shape()(0);
364 0 : ny = image->shape()(1);
365 0 : npol = image->shape()(2);
366 0 : nchan = image->shape()(3);
367 :
368 0 : if(image->shape().product()>cachesize)
369 0 : isTiled=true;
370 : else
371 0 : isTiled=false;
372 :
373 0 : sumWeight.resize(npol, nchan);
374 :
375 0 : wConvSize=max(1, nWPlanes_p);
376 0 : if (!convFuncCacheReady)
377 : {
378 0 : if (convSupport.shape()(0) != wConvSize)
379 : {
380 0 : convSupport.resize(wConvSize,1,1,true);
381 0 : convSupport=0;
382 : }
383 : }
384 :
385 0 : uvScale.resize(3);
386 0 : uvScale=0.0;
387 0 : uvScale(0)=Float(nx)*image->coordinates().increment()(0);
388 0 : uvScale(1)=Float(ny)*image->coordinates().increment()(1);
389 0 : uvScale(2)=Float(wConvSize)*abs(image->coordinates().increment()(0));
390 :
391 0 : uvOffset.resize(3);
392 0 : uvOffset(0)=nx/2;
393 0 : uvOffset(1)=ny/2;
394 0 : uvOffset(2)=0;
395 :
396 0 : if(gridder) delete gridder; gridder=0;
397 0 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
398 0 : uvScale, uvOffset,
399 0 : "SF");
400 :
401 : // Set up image cache needed for gridding.
402 0 : if(imageCache) delete imageCache; imageCache=0;
403 :
404 : // The tile size should be large enough that the
405 : // extended convolution function can fit easily
406 0 : if(isTiled)
407 : {
408 0 : Float tileOverlap=0.5;
409 0 : tilesize=min(256,tilesize);
410 0 : IPosition tileShape=IPosition(4,tilesize,tilesize,npol,nchan);
411 0 : Vector<Float> tileOverlapVec(4);
412 0 : tileOverlapVec=0.0;
413 0 : tileOverlapVec(0)=tileOverlap;
414 0 : tileOverlapVec(1)=tileOverlap;
415 : if (sizeof(long) < 4) // 32-bit machine
416 : {
417 : Int tmpCacheVal=static_cast<Int>(cachesize);
418 : imageCache=new LatticeCache <Complex> (*image, tmpCacheVal, tileShape,
419 : tileOverlapVec,
420 : (tileOverlap>0.0));
421 : }
422 : else // 64-bit machine
423 : {
424 0 : Long tmpCacheVal=cachesize;
425 0 : imageCache=new LatticeCache <Complex> (*image, tmpCacheVal, tileShape,
426 : tileOverlapVec,
427 0 : (tileOverlap>0.0));
428 : }
429 0 : }
430 :
431 : #if(USETABLES)
432 : Double StepSize;
433 : Int N=500000;
434 : StepSize = abs((((2*nx)/uvScale(0))/(sigma) +
435 : MAXPOINTINGERROR*1.745329E-02*(sigma)/3600.0))/N;
436 : if (!nwEij.isReady())
437 : {
438 : logIO() << LogOrigin("nPBWProjectFT","init")
439 : << "Making lookup table for exp function with a resolution of "
440 : << StepSize << " radians. "
441 : << "Memory used: " << sizeof(Float)*N/(1024.0*1024.0)<< " MB."
442 : << LogIO::NORMAL
443 : <<LogIO::POST;
444 :
445 : nwEij.setSigma(sigma);
446 : nwEij.initExpTable(N,StepSize);
447 : // ExpTab.build(N,StepSize);
448 :
449 : logIO() << LogOrigin("nPBWProjectFT","init")
450 : << "Making lookup table for complex exp function with a resolution of "
451 : << 2*M_PI/N << " radians. "
452 : << "Memory used: " << 2*sizeof(Float)*N/(1024.0*1024.0) << " MB."
453 : << LogIO::NORMAL
454 : << LogIO::POST;
455 : nwEij.initCExpTable(N);
456 : // CExpTab.build(N);
457 : }
458 : #endif
459 : // vpSJ->reset();
460 0 : paChangeDetector.reset();
461 0 : makingPSF = false;
462 0 : }
463 : //
464 : //---------------------------------------------------------------
465 : //
466 : // This is nasty, we should use CountedPointers here.
467 0 : nPBWProjectFT::~nPBWProjectFT()
468 : {
469 0 : if(imageCache) delete imageCache; imageCache=0;
470 0 : if(gridder) delete gridder; gridder=0;
471 :
472 0 : Int NCF=convFuncCache.nelements();
473 0 : for(Int i=0;i<NCF;i++) delete convFuncCache[i];
474 0 : NCF=convWeightsCache.nelements();
475 0 : for(Int i=0;i<NCF;i++) delete convWeightsCache[i];
476 0 : }
477 : //
478 : //---------------------------------------------------------------
479 : //
480 0 : Int nPBWProjectFT::makePBPolnCoords(CoordinateSystem& squintCoord,const VisBuffer&vb)
481 : {
482 : //
483 : // Make an image with circular polarization axis.
484 : //
485 0 : Int NPol=0,M,N=0;
486 0 : M=polMap.nelements();
487 0 : for(Int i=0;i<M;i++) if (polMap(i) > -1) NPol++;
488 0 : Vector<Int> poln(NPol);
489 :
490 : Int index;
491 0 : Vector<Int> inStokes;
492 0 : index = squintCoord.findCoordinate(Coordinate::STOKES);
493 0 : inStokes = squintCoord.stokesCoordinate(index).stokes();
494 0 : N = 0;
495 : try
496 : {
497 0 : for(Int i=0;i<M;i++) if (polMap(i) > -1) {poln(N) = vb.corrType()(i);N++;}
498 0 : StokesCoordinate polnCoord(poln);
499 0 : Int StokesIndex = squintCoord.findCoordinate(Coordinate::STOKES);
500 0 : squintCoord.replaceCoordinate(polnCoord,StokesIndex);
501 0 : cfStokes = poln;
502 0 : }
503 0 : catch(AipsError& x)
504 : {
505 0 : throw(SynthesisFTMachineError("Likely cause: Discrepancy between the poln. "
506 0 : "axis of the data and the image specifications."));
507 0 : }
508 :
509 0 : return NPol;
510 0 : }
511 : //
512 : //---------------------------------------------------------------
513 : //
514 0 : MDirection::Convert nPBWProjectFT::makeCoordinateMachine(const VisBuffer& vb,
515 : const MDirection::Types& From,
516 : const MDirection::Types& To,
517 : MEpoch& last)
518 : {
519 0 : Double time = getCurrentTimeStamp(vb);
520 :
521 0 : MEpoch epoch(Quantity(time,Second),MEpoch::TAI);
522 : // epoch = MEpoch(Quantity(time,Second),MEpoch::TAI);
523 : //
524 : // ...now make an object to hold the observatory position info...
525 : //
526 0 : MPosition pos;
527 0 : String ObsName=vb.msColumns().observation().telescopeName()(vb.arrayId());
528 :
529 0 : if (!MeasTable::Observatory(pos,ObsName))
530 0 : throw(AipsError("Observatory position for "+ ObsName + " not found"));
531 : //
532 : // ...now make a Frame object out of the observatory position and
533 : // time objects...
534 : //
535 0 : MeasFrame frame(epoch,pos);
536 : //
537 : // ...finally make the convert machine.
538 : //
539 0 : MDirection::Convert mac(MDirection::Ref(From,frame),
540 0 : MDirection::Ref(To,frame));
541 :
542 0 : MEpoch::Convert toLAST = MEpoch::Convert(MEpoch::Ref(MEpoch::TAI,frame),
543 0 : MEpoch::Ref(MEpoch::LAST,frame));
544 0 : last = toLAST(epoch);
545 :
546 0 : return mac;
547 0 : }
548 : //
549 : //---------------------------------------------------------------
550 : //
551 0 : int nPBWProjectFT::findPointingOffsets(const VisBuffer& vb,
552 : Array<Float> &l_off,
553 : Array<Float> &m_off,
554 : Bool Evaluate)
555 : {
556 :
557 : // throw(AipsError("PBWProject::findPointingOffsets temporarily disabled. (gmoellen 06Nov10)"));
558 :
559 0 : Int NAnt = 0;
560 0 : MEpoch LAST;
561 0 : Double thisTime = getCurrentTimeStamp(vb);
562 : // Array<Float> pointingOffsets = epJ->nearest(thisTime);
563 0 : if (epJ==NULL) return 0;
564 0 : Array<Float> pointingOffsets; epJ->nearest(thisTime,pointingOffsets);
565 0 : NAnt=pointingOffsets.shape()(2);
566 0 : l_off.resize(IPosition(3,1,1,NAnt)); // Poln x NChan x NAnt
567 0 : m_off.resize(IPosition(3,1,1,NAnt)); // Poln x NChan x NAnt
568 : // Can't figure out how to do the damn slicing of [Pol,NChan,NAnt,1] array
569 : // into [Pol,NChan,NAnt] array
570 : //
571 : // l_off = pointingOffsets(Slicer(IPosition(4,0,0,0,0),IPosition(4,1,1,NAnt+1,0)));
572 : // m_off = pointingOffsets(Slicer(IPosition(4,1,0,0,0),IPosition(4,1,1,NAnt+1,0)));
573 0 : IPosition tndx(3,0,0,0), sndx(4,0,0,0,0);
574 0 : for(tndx(2)=0;tndx(2)<NAnt; tndx(2)++,sndx(2)++)
575 : // for(Int j=0;j<NAnt;j++)
576 : {
577 : // l_off(IPosition(3,0,0,j)) = pointingOffsets(IPosition(4,0,0,j,0));
578 : // m_off(IPosition(3,0,0,j)) = pointingOffsets(IPosition(4,1,0,j,0));
579 :
580 0 : sndx(0)=0; l_off(tndx) = pointingOffsets(sndx);
581 0 : sndx(0)=2; m_off(tndx) = pointingOffsets(sndx);
582 : }
583 0 : return NAnt;
584 : if (!Evaluate) return NAnt;
585 :
586 : // cout << "AzEl Offsets: " << pointingOffsets << endl;
587 : //
588 : // Make a Coordinate Conversion Machine to go from (Az,El) to
589 : // (HA,Dec).
590 : //
591 : MDirection::Convert toAzEl = makeCoordinateMachine(vb,MDirection::HADEC,
592 : MDirection::AZEL,
593 : LAST);
594 : MDirection::Convert toHADec = makeCoordinateMachine(vb,MDirection::AZEL,
595 : MDirection::HADEC,
596 : LAST);
597 : //
598 : // ...and now hope that it all works and works correctly!!!
599 : //
600 : Quantity dAz(0,Radian),dEl(0,Radian);
601 : //
602 : // An array of shape [2,1,1]!
603 : //
604 : Array<Double> phaseDir = vb.msColumns().field().phaseDir().getColumn();
605 : Double RA0 = phaseDir(IPosition(3,0,0,0));
606 : Double Dec0 = phaseDir(IPosition(3,1,0,0));
607 : //
608 : // Compute reference (HA,Dec)
609 : //
610 : Double LST = LAST.get(Day).getValue();
611 : Double SDec0 = sin(Dec0), CDec0=cos(Dec0);
612 : LST -= floor(LST); // Extract the fractional day
613 : LST *= 2*C::pi;// Convert to Raidan
614 :
615 : Double HA0;
616 : HA0 = LST - RA0;
617 : Quantity QHA0(HA0,Radian), QDEC0(Dec0,Radian);
618 : //
619 : // Convert reference (HA,Dec) to reference (Az,El)
620 : //
621 : MDirection PhaseCenter(QHA0, QDEC0,MDirection::Ref(MDirection::HADEC));
622 : MDirection AzEl0 = toAzEl(PhaseCenter);
623 :
624 : MDirection tmpHADec = toHADec(AzEl0);
625 :
626 : Double Az0_Rad = AzEl0.getAngle(Radian).getValue()(0);
627 : Double El0_Rad = AzEl0.getAngle(Radian).getValue()(1);
628 :
629 : //
630 : // Convert the antenna pointing offsets from (Az,El)-->(RA,Dec)-->(l,m)
631 : //
632 :
633 : for(IPosition n(3,0,0,0);n(2)<=NAnt;n(2)++)
634 : {
635 : //
636 : // From (Az,El) -> (HA,Dec)
637 : //
638 : // Add (Az,El) offsets to the reference (Az,El)
639 : //
640 : dAz.setValue(l_off(n)+Az0_Rad); dEl.setValue(m_off(n)+El0_Rad);
641 : // dAz.setValue(0.0+Az0_Rad); dEl.setValue(0.0+El0_Rad);
642 : MDirection AzEl(dAz,dEl,MDirection::Ref(MDirection::AZEL));
643 : //
644 : // Convert offsetted (Az,El) to (HA,Dec) and then to (RA,Dec)
645 : //
646 : MDirection HADec = toHADec(AzEl);
647 : Double HA,Dec,RA, dRA;
648 : HA = HADec.getAngle(Radian).getValue()(0);
649 : Dec = HADec.getAngle(Radian).getValue()(1);
650 : RA = LST - HA;
651 : dRA = RA - RA0;
652 : //
653 : // Convert offsetted (RA,Dec) -> (l,m)
654 : //
655 : l_off(n) = sin(dRA)*cos(Dec);
656 : m_off(n) = sin(Dec)*CDec0-cos(Dec)*SDec0*cos(dRA);
657 : // cout << "FindOff: " << n(2) << " " << l_offsets(n) << " " << m_offsets(n) << endl;
658 :
659 : // cout << l_off(n) << " " << m_off(n) << " "
660 : // << " " << HA << " " << Dec
661 : // << " " << LST << " " << RA0 << " " << Dec0
662 : // << " " << RA << " " << Dec
663 : // << endl;
664 : }
665 :
666 : return NAnt+1;
667 0 : }
668 : //
669 : //---------------------------------------------------------------
670 : //
671 0 : int nPBWProjectFT::findPointingOffsets(const VisBuffer& vb,
672 : Cube<Float>& pointingOffsets,
673 : Array<Float> &l_off,
674 : Array<Float> &m_off,
675 : Bool Evaluate)
676 : {
677 0 : Int NAnt = 0;
678 : // TBD: adapt the following to VisCal mechanism:
679 0 : MEpoch LAST;
680 :
681 0 : NAnt=pointingOffsets.shape()(2);
682 0 : l_off.resize(IPosition(3,2,1,NAnt));
683 0 : m_off.resize(IPosition(3,2,1,NAnt));
684 0 : IPosition ndx(3,0,0,0),ndx1(3,0,0,0);
685 0 : for(ndx(2)=0;ndx(2)<NAnt;ndx(2)++)
686 : {
687 0 : ndx1=ndx;
688 0 : ndx(0)=0;ndx1(0)=0; l_off(ndx) = pointingOffsets(ndx1);//Axis_0,Pol_0,Ant_i
689 0 : ndx(0)=1;ndx1(0)=1; l_off(ndx) = pointingOffsets(ndx1);//Axis_0,Pol_1,Ant_i
690 0 : ndx(0)=0;ndx1(0)=2; m_off(ndx) = pointingOffsets(ndx1);//Axis_1,Pol_0,Ant_i
691 0 : ndx(0)=1;ndx1(0)=3; m_off(ndx) = pointingOffsets(ndx1);//Axis_1,Pol_1,Ant_i
692 : }
693 :
694 : // l_off = pointingOffsets(IPosition(3,0,0,0),IPosition(3,0,0,NAnt));
695 : // m_off = pointingOffsets(IPosition(3,1,0,0),IPosition(3,1,0,NAnt));
696 : /*
697 : IPosition shp(pointingOffsets.shape());
698 : IPosition shp1(l_off.shape()),shp2(m_off.shape());
699 : for(Int ii=0;ii<NAnt;ii++)
700 : {
701 : IPosition ndx(3,0,0,0);
702 : ndx(2)=ii;
703 : cout << "Pointing Offsets: " << ii << " "
704 : << l_off(ndx)*57.295*60.0 << " "
705 : << m_off(ndx)*57.295*60.0 << endl;
706 : }
707 : */
708 0 : return NAnt;
709 : if (!Evaluate) return NAnt;
710 :
711 : //
712 : // Make a Coordinate Conversion Machine to go from (Az,El) to
713 : // (HA,Dec).
714 : //
715 : MDirection::Convert toAzEl = makeCoordinateMachine(vb,MDirection::HADEC,
716 : MDirection::AZEL,
717 : LAST);
718 : MDirection::Convert toHADec = makeCoordinateMachine(vb,MDirection::AZEL,
719 : MDirection::HADEC,
720 : LAST);
721 : //
722 : // ...and now hope that it all works and works correctly!!!
723 : //
724 : Quantity dAz(0,Radian),dEl(0,Radian);
725 : //
726 : // An array of shape [2,1,1]!
727 : //
728 : Array<Double> phaseDir = vb.msColumns().field().phaseDir().getColumn();
729 : Double RA0 = phaseDir(IPosition(3,0,0,0));
730 : Double Dec0 = phaseDir(IPosition(3,1,0,0));
731 : //
732 : // Compute reference (HA,Dec)
733 : //
734 : Double LST = LAST.get(Day).getValue();
735 : Double SDec0 = sin(Dec0), CDec0=cos(Dec0);
736 : LST -= floor(LST); // Extract the fractional day
737 : LST *= 2*C::pi;// Convert to Raidan
738 :
739 : Double HA0;
740 : HA0 = LST - RA0;
741 : Quantity QHA0(HA0,Radian), QDEC0(Dec0,Radian);
742 : //
743 : // Convert reference (HA,Dec) to reference (Az,El)
744 : //
745 : MDirection PhaseCenter(QHA0, QDEC0,MDirection::Ref(MDirection::HADEC));
746 : MDirection AzEl0 = toAzEl(PhaseCenter);
747 :
748 : MDirection tmpHADec = toHADec(AzEl0);
749 :
750 : Double Az0_Rad = AzEl0.getAngle(Radian).getValue()(0);
751 : Double El0_Rad = AzEl0.getAngle(Radian).getValue()(1);
752 :
753 : //
754 : // Convert the antenna pointing offsets from (Az,El)-->(RA,Dec)-->(l,m)
755 : //
756 :
757 : for(IPosition n(3,0,0,0);n(2)<=NAnt;n(2)++)
758 : {
759 : //
760 : // From (Az,El) -> (HA,Dec)
761 : //
762 : // Add (Az,El) offsets to the reference (Az,El)
763 : //
764 : dAz.setValue(l_off(n)+Az0_Rad); dEl.setValue(m_off(n)+El0_Rad);
765 : // dAz.setValue(0.0+Az0_Rad); dEl.setValue(0.0+El0_Rad);
766 : MDirection AzEl(dAz,dEl,MDirection::Ref(MDirection::AZEL));
767 : //
768 : // Convert offsetted (Az,El) to (HA,Dec) and then to (RA,Dec)
769 : //
770 : MDirection HADec = toHADec(AzEl);
771 : Double HA,Dec,RA, dRA;
772 : HA = HADec.getAngle(Radian).getValue()(0);
773 : Dec = HADec.getAngle(Radian).getValue()(1);
774 : RA = LST - HA;
775 : dRA = RA - RA0;
776 : //
777 : // Convert offsetted (RA,Dec) -> (l,m)
778 : //
779 : l_off(n) = sin(dRA)*cos(Dec);
780 : m_off(n) = sin(Dec)*CDec0-cos(Dec)*SDec0*cos(dRA);
781 : // cout << "FindOff: " << n(2) << " " << l_offsets(n) << " " << m_offsets(n) << endl;
782 :
783 : // cout << l_off(n) << " " << m_off(n) << " "
784 : // << " " << HA << " " << Dec
785 : // << " " << LST << " " << RA0 << " " << Dec0
786 : // << " " << RA << " " << Dec
787 : // << endl;
788 : }
789 :
790 : return NAnt+1;
791 0 : }
792 : //
793 : //---------------------------------------------------------------
794 : //
795 0 : void nPBWProjectFT::makeSensitivityImage(Lattice<Complex>& wtImage,
796 : ImageInterface<Float>& sensitivityImage,
797 : const Matrix<Float>& sumWt,
798 : const Bool& doFFTNorm)
799 : {
800 0 : Bool doSumWtNorm=true;
801 0 : if (sumWt.shape().nelements()==0) doSumWtNorm=false;
802 :
803 0 : if ((sumWt.shape().nelements() < 2) ||
804 0 : (sumWt.shape()(0) != wtImage.shape()(2)) ||
805 0 : (sumWt.shape()(1) != wtImage.shape()(3)))
806 0 : throw(AipsError("makeSensitivityImage(): "
807 0 : "Sum of weights per poln and chan required"));
808 0 : Float sumWtVal=1.0;
809 :
810 0 : LatticeFFT::cfft2d(wtImage,false);
811 0 : Int sizeX=wtImage.shape()(0), sizeY=wtImage.shape()(1);
812 0 : sensitivityImage.resize(wtImage.shape());
813 0 : Array<Float> senBuf;
814 0 : sensitivityImage.get(senBuf,false);
815 0 : ArrayLattice<Float> senLat(senBuf, true);
816 :
817 : //
818 : // Copy one 2D plane at a time, normalizing by the sum of weights
819 : // and possibly 2D FFT.
820 : //
821 : // Set up Lattice iteratos on wtImage and sensitivityImage
822 : //
823 0 : IPosition axisPath(4, 0, 1, 2, 3);
824 0 : IPosition cursorShape(4, sizeX, sizeY, 1, 1);
825 0 : LatticeStepper wtImStepper(wtImage.shape(), cursorShape, axisPath);
826 0 : LatticeIterator<Complex> wtImIter(wtImage, wtImStepper);
827 0 : LatticeStepper senImStepper(senLat.shape(), cursorShape, axisPath);
828 0 : LatticeIterator<Float> senImIter(senLat, senImStepper);
829 : //
830 : // Iterate over channel and polarization axis
831 : //
832 0 : if (!doFFTNorm) sizeX=sizeY=1;
833 0 : for(wtImIter.reset(),senImIter.reset(); !wtImIter.atEnd(); wtImIter++,senImIter++)
834 : {
835 0 : Int pol=wtImIter.position()(2), chan=wtImIter.position()(3);
836 0 : if (doSumWtNorm) sumWtVal=sumWt(pol,chan);
837 0 : senImIter.rwCursor() = (real(wtImIter.rwCursor())
838 0 : *Float(sizeX)*Float(sizeY)/sumWtVal);
839 : }
840 : //
841 : // The following code is averaging RR and LL planes and writing
842 : // the result to back to both planes. This needs to be
843 : // generalized for full-pol case.
844 : //
845 0 : IPosition start0(4,0,0,0,0), start1(4,0,0,1,0), length(4,sizeX,sizeY,1,1);
846 0 : Slicer slicePol0(start0,length), slicePol1(start1,length);
847 0 : Array<Float> polPlane0, polPlane1;
848 0 : senLat.getSlice(polPlane0,slicePol0);
849 0 : senLat.getSlice(polPlane1,slicePol1);
850 0 : polPlane0=(polPlane0+polPlane1)/2.0;
851 0 : polPlane1=polPlane0;
852 : // senLat.putSlice(polPlane0,IPosition(4,0,0,0,0));
853 : // senLat.putSlice(polPlane1,IPosition(4,0,0,1,0));
854 : // cerr << "Pol0: " << polPlane0.shape() << " " << max(polPlane0) << endl;
855 : // cerr << "Pol1: " << polPlane1.shape() << " " << max(polPlane1) << endl;
856 0 : }
857 : //
858 : //---------------------------------------------------------------
859 : //
860 0 : void nPBWProjectFT::normalizeAvgPB()
861 : {
862 0 : if (!pbNormalized)
863 : {
864 0 : pbPeaks.resize(avgPB.shape()(2),true);
865 0 : if (makingPSF) pbPeaks = 1.0;
866 0 : else pbPeaks /= (Float)noOfPASteps;
867 0 : pbPeaks = 1.0;
868 0 : logIO() << LogOrigin("nPBWProjectFT", "normalizeAvgPB")
869 : << "Normalizing the average PBs to " << 1.0
870 : << LogIO::NORMAL
871 0 : << LogIO::POST;
872 :
873 0 : IPosition avgPBShape(avgPB.shape()),ndx(4,0,0,0,0);
874 0 : Vector<Float> peak(avgPBShape(2));
875 :
876 : Bool isRefF;
877 0 : Array<Float> avgPBBuf;
878 0 : isRefF=avgPB.get(avgPBBuf);
879 :
880 0 : Float pbMax = max(avgPBBuf);
881 :
882 0 : ndx=0;
883 0 : for(ndx(1)=0;ndx(1)<avgPBShape(1);ndx(1)++)
884 0 : for(ndx(0)=0;ndx(0)<avgPBShape(0);ndx(0)++)
885 : {
886 0 : IPosition plane1(ndx);
887 0 : plane1=ndx;
888 0 : plane1(2)=1; // The other poln. plane
889 0 : avgPBBuf(ndx) = (avgPBBuf(ndx) + avgPBBuf(plane1))/2.0;
890 0 : }
891 0 : for(ndx(1)=0;ndx(1)<avgPBShape(1);ndx(1)++)
892 0 : for(ndx(0)=0;ndx(0)<avgPBShape(0);ndx(0)++)
893 : {
894 0 : IPosition plane1(ndx);
895 0 : plane1=ndx;
896 0 : plane1(2)=1; // The other poln. plane
897 0 : avgPBBuf(plane1) = avgPBBuf(ndx);
898 0 : }
899 0 : if (fabs(pbMax-1.0) > 1E-3)
900 : {
901 : // avgPBBuf = avgPBBuf/noOfPASteps;
902 0 : for(ndx(3)=0;ndx(3)<avgPBShape(3);ndx(3)++)
903 0 : for(ndx(2)=0;ndx(2)<avgPBShape(2);ndx(2)++)
904 : {
905 0 : peak(ndx(2)) = 0;
906 0 : for(ndx(1)=0;ndx(1)<avgPBShape(1);ndx(1)++)
907 0 : for(ndx(0)=0;ndx(0)<avgPBShape(0);ndx(0)++)
908 0 : if (abs(avgPBBuf(ndx)) > peak(ndx(2)))
909 0 : peak(ndx(2)) = avgPBBuf(ndx);
910 :
911 0 : for(ndx(1)=0;ndx(1)<avgPBShape(1);ndx(1)++)
912 0 : for(ndx(0)=0;ndx(0)<avgPBShape(0);ndx(0)++)
913 0 : avgPBBuf(ndx) *= (pbPeaks(ndx(2))/peak(ndx(2)));
914 : }
915 0 : if (isRefF) avgPB.put(avgPBBuf);
916 : }
917 0 : }
918 0 : pbNormalized = true;
919 0 : }
920 : //
921 : //---------------------------------------------------------------
922 : //
923 0 : Bool nPBWProjectFT::makeAveragePB0(const VisBuffer& vb,
924 : const ImageInterface<Complex>& image,
925 : Int& /*polInUse*/,
926 : //TempImage<Float>& thesquintPB,
927 : TempImage<Float>& theavgPB)
928 : {
929 0 : TempImage<Float> localPB;
930 :
931 0 : logIO() << LogOrigin("nPBWProjecFT","makeAveragePB")
932 0 : << LogIO::NORMAL;
933 :
934 0 : localPB.resize(image.shape()); localPB.setCoordinateInfo(image.coordinates());
935 0 : localPB.setMaximumCacheSize(cachesize);
936 : // cerr << "Max. cache size = " << localPB.maximumCacheSize() << " " << cachesize << endl;
937 : //
938 : // If this is the first time, resize the average PB
939 : //
940 0 : if (resetPBs)
941 : {
942 : logIO() << "Initializing the average PBs"
943 : << LogIO::NORMAL
944 0 : << LogIO::POST;
945 0 : theavgPB.resize(localPB.shape());
946 0 : theavgPB.setCoordinateInfo(localPB.coordinates());
947 0 : theavgPB.set(0.0);
948 0 : noOfPASteps = 0;
949 0 : pbPeaks.resize(theavgPB.shape()(2));
950 0 : pbPeaks.set(0.0);
951 0 : resetPBs=false;
952 : }
953 : //
954 : // Make the Stokes PB
955 : //
956 0 : localPB.set(1.0);
957 :
958 : // The following replaces the simple vlaPB.applyPB() call.
959 : // The functional interface of VLACalcIllumination was
960 : // modified to apply one polarization PB at a time. As a
961 : // result, the co-ord. sys. of the image going to applyPB()
962 : // should have only one Poln. plane.
963 : //
964 : // Before this change, the localPB image was directly sent to
965 : // applyPB(). Now, in the code below, we go over the poln. planes
966 : // of this image, make a temp. image with one poln. planes of
967 : // localPB at a time, applyPB() on the temp. image, and copy the
968 : // result in the appropriate planes of the localPB image. The rest
969 : // of the code therefore does not see a difference.
970 : {
971 0 : VLACalcIlluminationConvFunc vlaPB;
972 0 : if (bandID_p == -1) bandID_p=getVisParams(vb);
973 0 : CoordinateSystem coords=localPB.coordinates();
974 : //----------------------------------------------------------------------
975 0 : IPosition PolnPlane(4,0,0,0,0);
976 0 : CoordinateSystem singlePolCoords=coords;
977 0 : Int index, outNdx, whichPolPlane=0;
978 0 : Vector<Int> inStokes, outStokes(1);
979 0 : index = coords.findCoordinate(Coordinate::STOKES);
980 0 : inStokes = coords.stokesCoordinate(index).stokes();
981 0 : for (uInt i=0; i<inStokes.nelements(); i++)
982 : {
983 0 : outStokes(0)=inStokes(i);
984 : // Make a temp. image with a single Stokes along Pol. axis.
985 0 : StokesCoordinate polnCoord(outStokes);
986 0 : outNdx = singlePolCoords.findCoordinate(Coordinate::STOKES);
987 0 : singlePolCoords.replaceCoordinate(polnCoord, outNdx);
988 0 : IPosition singlePolShape = localPB.shape();
989 0 : singlePolShape(2)=1;
990 0 : Bool doSquint=false;
991 : {
992 0 : TempImage<Float> singlePolImg(singlePolShape, singlePolCoords);
993 : // Copy screen to the single pol. image. Apply PB to it. Copy
994 : // the single pol. image plane to the twoDPB image.
995 0 : singlePolImg.set(1.0);
996 0 : Double pa=getPA(vb);
997 0 : Vector<Double> chanFreq = vb.frequency();
998 : //freqHi = max(chanFreq);
999 0 : Double Freq = max(chanFreq);
1000 :
1001 0 : vlaPB.applyPB(singlePolImg, pa, bandID_p, doSquint, Freq);
1002 0 : PolnPlane(2)=whichPolPlane; localPB.putSlice(singlePolImg.get(), PolnPlane);
1003 0 : }
1004 0 : whichPolPlane++;
1005 0 : }
1006 : //----------------------------------------------------------------------
1007 :
1008 : // vlaPB.applyPB(localPB, vb, bandID_p);
1009 0 : }
1010 :
1011 0 : IPosition twoDPBShape(localPB.shape());
1012 0 : TempImage<Complex> localTwoDPB(twoDPBShape,localPB.coordinates());
1013 0 : localTwoDPB.setMaximumCacheSize(cachesize);
1014 0 : Float peak=0;
1015 : Int NAnt;
1016 0 : noOfPASteps++;
1017 0 : NAnt=1;
1018 :
1019 : // logIO() << " Shape of localPB Cube : " << twoDPBShape << LogIO::POST;
1020 : // logIO() << " Shape of avgPB Cube : " << theavgPB.shape() << LogIO::POST;
1021 :
1022 0 : for(Int ant=0;ant<NAnt;ant++)
1023 : { //Ant loop
1024 : {
1025 0 : IPosition ndx(4,0,0,0,0);
1026 0 : for(ndx(0)=0; ndx(0)<twoDPBShape(0); ndx(0)++)
1027 0 : for(ndx(1)=0; ndx(1)<twoDPBShape(1); ndx(1)++)
1028 0 : for(ndx(2)=0; ndx(2)<twoDPBShape(2); ndx(2)++)
1029 0 : for(ndx(3)=0; ndx(3)<twoDPBShape(3); ndx(3)++)
1030 0 : localTwoDPB.putAt(Complex((localPB(ndx)),0.0),ndx);
1031 0 : }
1032 : //
1033 : // If antenna pointing errors are not applied, no shifting
1034 : // (which can be expensive) is required.
1035 : //
1036 : //
1037 : // Accumulate the shifted PBs
1038 : //
1039 : {
1040 : Bool isRefF;
1041 0 : Array<Float> fbuf;
1042 0 : Array<Complex> cbuf;
1043 0 : isRefF=theavgPB.get(fbuf);
1044 : //isRefC=
1045 0 : localTwoDPB.get(cbuf);
1046 :
1047 0 : IPosition fs(fbuf.shape());
1048 : {
1049 0 : IPosition ndx(4,0,0,0,0),avgNDX(4,0,0,0,0);
1050 0 : for(ndx(3)=0,avgNDX(3)=0;ndx(3)<fs(3);ndx(3)++,avgNDX(3)++)
1051 : {
1052 0 : for(ndx(2)=0,avgNDX(2)=0;ndx(2)<twoDPBShape(2);ndx(2)++,avgNDX(2)++)
1053 : {
1054 0 : for(ndx(0)=0,avgNDX(0)=0;ndx(0)<fs(0);ndx(0)++,avgNDX(0)++)
1055 0 : for(ndx(1)=0,avgNDX(1)=0;ndx(1)<fs(1);ndx(1)++,avgNDX(1)++)
1056 : {
1057 : Float val;
1058 0 : val = real(cbuf(ndx));
1059 0 : fbuf(avgNDX) += val;
1060 0 : if (fbuf(avgNDX) > peak) peak=fbuf(avgNDX);
1061 : }
1062 : }
1063 : }
1064 0 : }
1065 0 : if (!isRefF) theavgPB.put(fbuf);
1066 0 : pbPeaks += peak;
1067 0 : }
1068 : }
1069 0 : theavgPB.setCoordinateInfo(localPB.coordinates());
1070 0 : return true; // i.e., an average PB was made and is in the mem. cache
1071 0 : }
1072 : //
1073 : //---------------------------------------------------------------
1074 : //
1075 : //
1076 : //---------------------------------------------------------------
1077 : //
1078 : // Locate a convlution function in either mem. or disk cache.
1079 : // Return 1 if found in the disk cache.
1080 : // 2 if found in the mem. cache.
1081 : // <0 if not found in either cache. In this case, absolute of
1082 : // the return value corresponds to the index in the list of
1083 : // conv. funcs. where this conv. func. should be filled
1084 : //
1085 0 : Int nPBWProjectFT::locateConvFunction(Int Nw, Int /*polInUse*/,
1086 : const VisBuffer& vb, Float &pa)
1087 : {
1088 : Int i; Bool found;
1089 : // Commented out to prevent compiler warning
1090 : // Float dPA = paChangeDetector.getParAngleTolerance().getValue("rad");
1091 0 : found = cfCache.searchConvFunction(vb,paChangeDetector,i,pa);
1092 0 : if (found)
1093 : {
1094 0 : Vector<Float> sampling;
1095 0 : PAIndex=i;
1096 : // CoordinateSystem csys;
1097 0 : if (cfCache.loadConvFunction(i,Nw,convFuncCache,convSupport,sampling,
1098 0 : cfRefFreq_p,convFuncCS_p))
1099 : {
1100 0 : convSampling = (Int)sampling[0];
1101 0 : convFunc.reference(*convFuncCache[PAIndex]);
1102 0 : cfCache.loadConvFunction(i,Nw,convWeightsCache,convSupport,sampling,cfRefFreq_p,
1103 0 : convFuncCS_p, "/CFWT");
1104 0 : convWeights.reference(*convWeightsCache[PAIndex]);
1105 0 : if (PAIndex < (Int)convFuncCache.nelements())
1106 : logIO() << "Loaded from disk cache: Conv. func. # "
1107 0 : << PAIndex << LogIO::POST;
1108 0 : return 1;
1109 : }
1110 0 : convFunc.reference(*convFuncCache[PAIndex]);
1111 0 : convWeights.reference(*convWeightsCache[PAIndex]);
1112 0 : return 2;
1113 0 : }
1114 0 : return i;
1115 : }
1116 0 : void nPBWProjectFT::makeCFPolMap(const VisBuffer& vb, Vector<Int>& polM)
1117 : {
1118 0 : Vector<Int> msStokes = vb.corrType();
1119 0 : Int nPol = msStokes.nelements();
1120 0 : polM.resize(polMap.shape());
1121 0 : polM = -1;
1122 :
1123 0 : for(Int i=0;i<nPol;i++)
1124 0 : for(uInt j=0;j<cfStokes.nelements();j++)
1125 0 : if (cfStokes(j) == msStokes(i))
1126 0 : {polM(i) = j;break;}
1127 0 : }
1128 : //
1129 : //---------------------------------------------------------------
1130 : //
1131 : // Given a polMap (mapping of which Visibility polarization is
1132 : // gridded onto which grid plane), make a map of the conjugate
1133 : // planes of the grid E.g, for Stokes-I and -V imaging, the two
1134 : // planes of the uv-grid are [LL,RR]. For input VisBuffer
1135 : // visibilites in order [RR,RL,LR,LL], polMap = [1,-1,-1,0]. The
1136 : // conjugate map will be [0,-1,-1,1].
1137 : //
1138 0 : void nPBWProjectFT::makeConjPolMap(const VisBuffer& vb,
1139 : const Vector<Int> cfPolMap,
1140 : Vector<Int>& conjPolMap)
1141 : {
1142 : //
1143 : // All the Natak (Drama) below with slicers etc. is to extract the
1144 : // Poln. info. for the first IF only (not much "information
1145 : // hiding" for the code to slice arrays in a general fashion).
1146 : //
1147 : // Extract the shape of the array to be sliced.
1148 : //
1149 0 : Array<Int> stokesForAllIFs = vb.msColumns().polarization().corrType().getColumn();
1150 0 : IPosition stokesShape(stokesForAllIFs.shape());
1151 0 : IPosition firstIFStart(stokesShape),firstIFLength(stokesShape);
1152 : //
1153 : // Set up the start and length IPositions to extract only the
1154 : // first column of the array. The following is required since the
1155 : // array could have only one column as well.
1156 : //
1157 0 : firstIFStart(0)=0;firstIFLength(0)=stokesShape(0);
1158 0 : for(uInt i=1;i<stokesShape.nelements();i++) {firstIFStart(i)=0;firstIFLength(i)=1;}
1159 : //
1160 : // Construct the slicer and produce the slice. .nonDegenerate
1161 : // required to ensure the result of slice is a pure vector.
1162 : //
1163 0 : Vector<Int> visStokes = stokesForAllIFs(Slicer(firstIFStart,firstIFLength)).nonDegenerate();
1164 :
1165 0 : conjPolMap = cfPolMap;
1166 :
1167 0 : Int i,j,N = cfPolMap.nelements();
1168 0 : for(i=0;i<N;i++)
1169 0 : if (cfPolMap[i] > -1)
1170 : {
1171 0 : if (visStokes[i] == Stokes::RR)
1172 : {
1173 0 : conjPolMap[i]=-1;
1174 0 : for(j=0;j<N;j++) if (visStokes[j] == Stokes::LL) break;
1175 0 : conjPolMap[i]=cfPolMap[j];
1176 : }
1177 0 : else if (visStokes[i] == Stokes::LL)
1178 : {
1179 0 : conjPolMap[i]=-1;
1180 0 : for(j=0;j<N;j++) if (visStokes[j] == Stokes::RR) break;
1181 0 : conjPolMap[i]=cfPolMap[j];
1182 : }
1183 0 : else if (visStokes[i] == Stokes::LR)
1184 : {
1185 0 : conjPolMap[i]=-1;
1186 0 : for(j=0;j<N;j++) if (visStokes[j] == Stokes::RL) break;
1187 0 : conjPolMap[i]=cfPolMap[j];
1188 : }
1189 0 : else if (visStokes[i] == Stokes::RL)
1190 : {
1191 0 : conjPolMap[i]=-1;
1192 0 : for(j=0;j<N;j++) if (visStokes[j] == Stokes::LR) break;
1193 0 : conjPolMap[i]=cfPolMap[j];
1194 : }
1195 : }
1196 0 : }
1197 : //
1198 : //---------------------------------------------------------------
1199 : //
1200 0 : void nPBWProjectFT::findConvFunction(const ImageInterface<Complex>& image,
1201 : const VisBuffer& vb)
1202 : {
1203 0 : if (!paChangeDetector.changed(vb,0)) return;
1204 :
1205 0 : logIO() << LogOrigin("nPBWProjectFT", "findConvFunction") << LogIO::NORMAL;
1206 :
1207 0 : ok();
1208 :
1209 :
1210 0 : CoordinateSystem coords(image.coordinates());
1211 :
1212 : //
1213 : // Make a two dimensional image to calculate auto-correlation of
1214 : // the ideal illumination pattern. We want this on a fine grid in
1215 : // the UV plane
1216 : //
1217 0 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
1218 0 : AlwaysAssert(directionIndex>=0, AipsError);
1219 0 : DirectionCoordinate dc=coords.directionCoordinate(directionIndex);
1220 0 : directionCoord=coords.directionCoordinate(directionIndex);
1221 0 : Vector<Double> sampling;
1222 0 : sampling = dc.increment();
1223 0 : sampling*=Double(convSampling);
1224 0 : sampling*=Double(nx)/Double(convSize);
1225 0 : dc.setIncrement(sampling);
1226 :
1227 :
1228 0 : Vector<Double> unitVec(2);
1229 0 : unitVec=convSize/2;
1230 0 : dc.setReferencePixel(unitVec);
1231 :
1232 : // Set the reference value to that of the image
1233 0 : coords.replaceCoordinate(dc, directionIndex);
1234 :
1235 : //
1236 : // Make an image with circular polarization axis. Return the
1237 : // no. of vis. poln. planes that will be used in making the user
1238 : // defined Stokes image.
1239 : //
1240 0 : polInUse=makePBPolnCoords(coords,vb);
1241 :
1242 : Float pa;
1243 0 : Int cfSource=locateConvFunction(wConvSize, polInUse, vb, pa);
1244 : // cerr << "CFS: " << cfSource << " " << PAIndex << " "
1245 : // << convFuncCache.nelements() << " "
1246 : // << convWeightsCache.nelements() << " "
1247 : // << endl;
1248 0 : lastPAUsedForWtImg = currentCFPA = pa;
1249 :
1250 0 : Bool pbMade=false;
1251 0 : if (cfSource==1) // CF found and loaded from the disk cache
1252 : {
1253 : // cout << "### New CFPA = " << currentCFPA << endl;
1254 0 : polInUse = convFunc.shape()(3);
1255 0 : wConvSize = convFunc.shape()(2);
1256 : try
1257 : {
1258 0 : cfCache.loadAvgPB(avgPB);
1259 0 : avgPBReady=true;
1260 : }
1261 0 : catch (AipsError& err)
1262 : {
1263 : logIO() << "Average PB does not exist in the cache. A fresh one will be made."
1264 0 : << LogIO::NORMAL << LogIO::POST;
1265 0 : pbMade=makeAveragePB0(vb, image, polInUse, avgPB);
1266 0 : pbNormalized=false; normalizeAvgPB(); pbNormalized=true;
1267 0 : }
1268 : }
1269 0 : else if (cfSource==2) // CF found in the mem. cache
1270 : {
1271 : }
1272 : else // CF not found in either cache
1273 : {
1274 : //
1275 : // Make the CF, update the average PB and update the CF and
1276 : // the avgPB disk cache
1277 : //
1278 0 : PAIndex = abs(cfSource);
1279 : //
1280 : // Load the average PB from the disk since it's going to be
1281 : // updated in memory and on the disk. Without loading it from
1282 : // the disk (from a potentially more complete existing cache),
1283 : // the average PB can get inconsistant with the rest of the
1284 : // cache.
1285 : //
1286 : // logIO() << LogOrigin("nPBWProjectFT::findConvFunction()","")
1287 : // << "Making the convolution function for PA=" << pa << "deg."
1288 : // << LogIO::NORMAL
1289 : // << LogIO::POST;
1290 0 : makeConvFunction(image,vb,pa);
1291 : try
1292 : {
1293 0 : cfCache.loadAvgPB(avgPB);
1294 0 : resetPBs = false;
1295 0 : avgPBReady=true;
1296 : }
1297 0 : catch(SynthesisFTMachineError &err)
1298 : {
1299 0 : logIO() << LogOrigin("nPBWProjectFT::findConvFunction()","")
1300 : << "Average PB does not exist in the cache. A fresh one will be made."
1301 : << LogIO::NORMAL
1302 0 : << LogIO::POST;
1303 0 : pbMade=makeAveragePB0(vb, image, polInUse,avgPB);
1304 0 : }
1305 :
1306 : // makeAveragePB(vb, image, polInUse,avgPB);
1307 0 : pbNormalized=false;
1308 0 : normalizeAvgPB();
1309 0 : pbNormalized=true;
1310 0 : Int index=coords.findCoordinate(Coordinate::SPECTRAL);
1311 0 : SpectralCoordinate spCS = coords.spectralCoordinate(index);
1312 0 : Vector<Double> refValue; refValue.resize(1);refValue(0)=cfRefFreq_p;
1313 0 : spCS.setReferenceValue(refValue);
1314 0 : coords.replaceCoordinate(spCS,index);
1315 0 : cfCache.cacheConvFunction(PAIndex, pa, convFunc, coords, convFuncCS_p,
1316 0 : convSize, convSupport,convSampling);
1317 0 : Cube<Int> convWtSize=convSupport*CONVWTSIZEFACTOR;
1318 0 : cfCache.cacheConvFunction(PAIndex, pa, convWeights, coords, convFuncCS_p,
1319 0 : convSize, convWtSize,convSampling,"WT",false);
1320 0 : cfCache.finalize(); // Write the aux info file
1321 0 : if (pbMade) cfCache.finalize(avgPB); // Save the AVG PB and write the aux info.
1322 0 : }
1323 :
1324 0 : verifyShapes(avgPB.shape(), image.shape());
1325 :
1326 0 : Int lastPASlot = PAIndex;
1327 :
1328 0 : if (paChangeDetector.changed(vb,0)) paChangeDetector.update(vb,0);
1329 : //
1330 : // If mem. cache not yet ready and the latest CF was loaded from
1331 : // the disk cache, compute and give some user useful info.
1332 : //
1333 0 : if ((!convFuncCacheReady) && (cfSource != 2))
1334 : {
1335 : //
1336 : // Compute the aggregate memory used by the cached convolution
1337 : // functions.
1338 : //
1339 0 : Int maxMemoryMB=HostInfo::memoryTotal(true)/1024;
1340 0 : Float memoryKB=0;
1341 0 : String unit(" KB");
1342 0 : for(uInt iPA=0;iPA<convFuncCache.nelements();iPA++)
1343 : {
1344 0 : Int volume=1;
1345 0 : if (convFuncCache[iPA] != NULL)
1346 : {
1347 0 : IPosition shape=(*convFuncCache[iPA]).shape();
1348 0 : volume=1;
1349 0 : for(uInt i=0;i<shape.nelements();i++)
1350 0 : volume*=shape(i);
1351 0 : memoryKB += volume;
1352 0 : }
1353 : }
1354 :
1355 :
1356 0 : memoryKB = Int(memoryKB*sizeof(Complex)*2/1024.0+0.5);
1357 0 : if (memoryKB > 1024) {memoryKB /=1024; unit=" MB";}
1358 :
1359 : logIO() << "Memory used in gridding functions = "
1360 0 : << (Int)(memoryKB+0.5) << unit << " out of a maximum of "
1361 0 : << maxMemoryMB << " MB" << LogIO::POST;
1362 :
1363 : //
1364 : // Show the list of support sizes along the w-axis for the current PA.
1365 : //
1366 0 : IPosition sliceStart(3,0,0,lastPASlot),
1367 0 : sliceLength(3,wConvSize,1,1);
1368 : logIO() << "Convolution support [CF#= " << lastPASlot
1369 0 : << ", PA=" << pa*180.0/M_PI << "deg"
1370 : <<"] = "
1371 0 : << convSupport(Slicer(sliceStart,sliceLength)).nonDegenerate()
1372 : << " pixels in Fourier plane"
1373 0 : << LogIO::POST;
1374 0 : }
1375 :
1376 0 : IPosition shp(convFunc.shape());
1377 0 : IPosition ndx(shp);
1378 0 : ndx =0;
1379 0 : Area.resize(Area.nelements()+1,true);
1380 0 : Area(lastPASlot)=0;
1381 :
1382 0 : for(ndx(0)=0;ndx(0)<shp(0);ndx(0)++)
1383 0 : for(ndx(1)=0;ndx(1)<shp(1);ndx(1)++)
1384 0 : Area(lastPASlot)+=convFunc(ndx);
1385 0 : }
1386 : //
1387 : //---------------------------------------------------------------
1388 : //
1389 0 : void nPBWProjectFT::makeConvFunction(const ImageInterface<Complex>& image,
1390 : const VisBuffer& vb,Float pa)
1391 : {
1392 0 : if (bandID_p == -1) bandID_p=getVisParams(vb);
1393 0 : Int NNN=0;
1394 0 : logIO() << LogOrigin("nPBWProjectFT", "makeConvFunction")
1395 : << "Making a new convolution function for PA="
1396 0 : << pa*(180/C::pi) << "deg"
1397 : << LogIO::NORMAL
1398 0 : << LogIO::POST;
1399 :
1400 0 : if(wConvSize>NNN) {
1401 0 : logIO() << "Using " << wConvSize << " planes for W-projection" << LogIO::POST;
1402 : Double maxUVW;
1403 0 : maxUVW=0.25/abs(image.coordinates().increment()(0));
1404 : logIO() << "Estimating maximum possible W = " << maxUVW
1405 0 : << " (wavelengths)" << LogIO::POST;
1406 :
1407 0 : Double invLambdaC=vb.frequency()(0)/C::c;
1408 : // logIO() << "Typical wavelength = " << 1.0/invLambdaC
1409 : // << " (m)" << LogIO::POST;
1410 0 : Double invMinL = vb.frequency()((vb.frequency().nelements())-1)/C::c;
1411 : logIO() << "wavelength range = " << 1.0/invLambdaC << " (m) to "
1412 0 : << 1.0/invMinL << " (m)" << LogIO::POST;
1413 0 : if (wConvSize > 1)
1414 : {
1415 0 : uvScale(2)=Float((wConvSize-1)*(wConvSize-1))/maxUVW;
1416 0 : logIO() << "Scaling in W (at maximum W) = " << 1.0/uvScale(2)
1417 0 : << " wavelengths per pixel" << LogIO::POST;
1418 : }
1419 : }
1420 : //
1421 : // Get the coordinate system
1422 : //
1423 0 : CoordinateSystem coords(image.coordinates());
1424 :
1425 : //
1426 : // Set up the convolution function.
1427 : //
1428 0 : if(wConvSize>NNN)
1429 : {
1430 0 : if(wConvSize>256)
1431 : {
1432 0 : convSampling=4;
1433 0 : convSize=min(nx,512);
1434 : }
1435 : else
1436 : {
1437 0 : convSampling=4;
1438 0 : convSize=min(nx,2048);
1439 : }
1440 : }
1441 : else
1442 : {
1443 0 : convSampling=4;
1444 0 : convSize=nx;
1445 : }
1446 0 : convSampling=OVERSAMPLING;
1447 0 : convSize=CONVSIZE;
1448 : //
1449 : // Make a two dimensional image to calculate auto-correlation of
1450 : // the ideal illumination pattern. We want this on a fine grid in
1451 : // the UV plane
1452 : //
1453 0 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
1454 0 : AlwaysAssert(directionIndex>=0, AipsError);
1455 0 : DirectionCoordinate dc=coords.directionCoordinate(directionIndex);
1456 0 : directionCoord=coords.directionCoordinate(directionIndex);
1457 0 : Vector<Double> sampling;
1458 0 : sampling = dc.increment();
1459 : // sampling /= Double(2.0);
1460 :
1461 0 : sampling*=Double(convSampling);
1462 0 : sampling*=Double(nx)/Double(convSize);
1463 :
1464 : // cerr << "Sampling on the sky = " << dc.increment() << " " << nx << "x" << ny << endl;
1465 : // cerr << "Sampling on the PB = " << sampling << " " << convSize << "x" << convSize << endl;
1466 0 : dc.setIncrement(sampling);
1467 :
1468 0 : Vector<Double> unitVec(2);
1469 0 : unitVec=convSize/2;
1470 0 : dc.setReferencePixel(unitVec);
1471 :
1472 : // Set the reference value to that of the image
1473 0 : coords.replaceCoordinate(dc, directionIndex);
1474 :
1475 : //
1476 : // Make an image with circular polarization axis. Return the
1477 : // no. of vis. poln. planes that will be used in making the user
1478 : // defined Stokes image.
1479 : //
1480 0 : polInUse=makePBPolnCoords(coords,vb);
1481 0 : convSupport.resize(wConvSize,polInUse,PAIndex+1,true);
1482 0 : Int N=convFuncCache.nelements();
1483 0 : convFuncCache.resize(PAIndex+1,true);
1484 0 : for(Int i=N;i<PAIndex;i++) convFuncCache[i]=NULL;
1485 0 : convWeightsCache.resize(PAIndex+1,true);
1486 0 : for(Int i=N;i<PAIndex;i++) convWeightsCache[i]=NULL;
1487 : //------------------------------------------------------------------
1488 : //
1489 : // Make the sky Stokes PB. This will be used in the gridding
1490 : // correction
1491 : //
1492 : //------------------------------------------------------------------
1493 0 : IPosition pbShape(4, convSize, convSize, polInUse, 1);
1494 0 : TempImage<Complex> twoDPB(pbShape, coords);
1495 :
1496 0 : IPosition pbSqShp(pbShape);
1497 : // pbSqShp[0] *= 2; pbSqShp[1] *= 2;
1498 :
1499 0 : unitVec=pbSqShp[0]/2;
1500 0 : dc.setReferencePixel(unitVec);
1501 : // sampling *= Double(2.0);
1502 : // dc.setIncrement(sampling);
1503 0 : coords.replaceCoordinate(dc, directionIndex);
1504 :
1505 0 : TempImage<Complex> twoDPBSq(pbSqShp,coords);
1506 0 : twoDPB.setMaximumCacheSize(cachesize);
1507 0 : twoDPB.set(Complex(1.0,0.0));
1508 0 : twoDPBSq.setMaximumCacheSize(cachesize);
1509 0 : twoDPBSq.set(Complex(1.0,0.0));
1510 : //
1511 : // Accumulate the various terms that constitute the gridding
1512 : // convolution function.
1513 : //
1514 : //------------------------------------------------------------------
1515 :
1516 0 : Int inner=convSize/convSampling;
1517 : // inner = convSize/2;
1518 :
1519 0 : Vector<Double> cfUVScale(3,0),cfUVOffset(3,0);
1520 :
1521 0 : cfUVScale(0)=Float(twoDPB.shape()(0))*sampling(0);
1522 0 : cfUVScale(1)=Float(twoDPB.shape()(1))*sampling(1);
1523 0 : cfUVOffset(0)=Float(twoDPB.shape()(0))/2;
1524 0 : cfUVOffset(1)=Float(twoDPB.shape()(1))/2;
1525 : // cerr << uvScale << " " << cfUVScale << endl;
1526 : // cerr << uvOffset << " " << cfUVOffset << endl;
1527 : ConvolveGridder<Double, Complex>
1528 : // ggridder(IPosition(2, inner, inner), cfUVScale, cfUVOffset, "SF");
1529 0 : ggridder(IPosition(2, inner, inner), uvScale, uvOffset, "SF");
1530 :
1531 : // convFuncCache[PAIndex] = new Array<Complex>(IPosition(4,convSize/2,convSize/2,
1532 : // wConvSize,polInUse));
1533 : // convWeightsCache[PAIndex] = new Array<Complex>(IPosition(4,convSize/2,convSize/2,
1534 : // wConvSize,polInUse));
1535 0 : convFuncCache[PAIndex] = new Array<Complex>(IPosition(4,convSize,convSize,
1536 0 : wConvSize,polInUse));
1537 0 : convWeightsCache[PAIndex] = new Array<Complex>(IPosition(4,convSize,convSize,
1538 0 : wConvSize,polInUse));
1539 0 : convFunc.reference(*convFuncCache[PAIndex]);
1540 0 : convWeights.reference(*convWeightsCache[PAIndex]);
1541 0 : convFunc=0.0;
1542 0 : convWeights=0.0;
1543 :
1544 0 : IPosition start(4, 0, 0, 0, 0);
1545 0 : IPosition pbSlice(4, convSize, convSize, 1, 1);
1546 :
1547 0 : Matrix<Complex> screen(convSize, convSize);
1548 0 : if (paChangeDetector.changed(vb,0)) paChangeDetector.update(vb,0);
1549 0 : VLACalcIlluminationConvFunc vlaPB;
1550 0 : vlaPB.setMaximumCacheSize(cachesize);
1551 :
1552 0 : for (Int iw=0;iw<wConvSize;iw++)
1553 : {
1554 0 : screen = 1.0;
1555 :
1556 : /*
1557 : screen=0.0;
1558 : // First the spheroidal function
1559 : //
1560 : // inner=convSize/2;
1561 : // screen = 0.0;
1562 : Vector<Complex> correction(inner);
1563 : for (Int iy=-inner/2;iy<inner/2;iy++)
1564 : {
1565 : ggridder.correctX1D(correction, iy+inner/2);
1566 : for (Int ix=-inner/2;ix<inner/2;ix++)
1567 : screen(ix+convSize/2,iy+convSize/2)=correction(ix+inner/2);
1568 : // if (iy==0)
1569 : // for (Int ii=0;ii<inner;ii++)
1570 : // cout << ii << " " << correction(ii) << endl;
1571 : }
1572 : */
1573 : //
1574 : // Now the w term
1575 : //
1576 0 : if(wConvSize>1)
1577 : {
1578 0 : logIO() << LogOrigin("nPBWProjectFT", "")
1579 0 : << "Computing WPlane " << iw << LogIO::POST;
1580 :
1581 0 : Double twoPiW=2.0*C::pi*Double(iw*iw)/uvScale(2);
1582 :
1583 0 : for (Int iy=-inner/2;iy<inner/2;iy++)
1584 : {
1585 0 : Double m=sampling(1)*Double(iy);
1586 0 : Double msq=m*m;
1587 0 : for (Int ix=-inner/2;ix<inner/2;ix++)
1588 : {
1589 0 : Double l=sampling(0)*Double(ix);
1590 0 : Double rsq=l*l+msq;
1591 0 : if(rsq<1.0)
1592 : {
1593 0 : Double phase=twoPiW*(sqrt(1.0-rsq)-1.0);
1594 0 : screen(ix+convSize/2,iy+convSize/2)*=Complex(cos(phase),sin(phase));
1595 : }
1596 : }
1597 : }
1598 : }
1599 :
1600 : //
1601 : // Fill the complex image with the w-term...
1602 : //
1603 0 : IPosition PolnPlane(4,0,0,0,0);
1604 0 : IPosition ndx(4,0,0,0,0);
1605 :
1606 0 : for(Int i=0;i<polInUse;i++)
1607 : {
1608 0 : PolnPlane(2)=i;
1609 0 : twoDPB.putSlice(screen, PolnPlane);
1610 0 : twoDPBSq.putSlice(screen, PolnPlane);
1611 : }
1612 : // {
1613 : // Vector<IPosition> posMax(twoDPB.shape()(2));
1614 : // posMax(0)(0)=pbShape(0)/2;
1615 : // posMax(0)(1)=pbShape(1)/2;
1616 : // posMax(1)(0)=pbShape(0)/2;
1617 : // posMax(1)(1)=pbShape(1)/2;
1618 : // getVisParams(vb);
1619 : // applyAntiAliasingOp(twoDPB,posMax,0);
1620 : // }
1621 : //
1622 : // Apply the PB...
1623 : //
1624 :
1625 : // The following replaces the simple vlaPB.applyPB() call.
1626 : // The functional interface of VLACalcIllumination was
1627 : // modified to apply one polarization PB at a time. As a
1628 : // result, the co-ord. sys. of the image going to applyPB()
1629 : // should have only one Poln. plane.
1630 : //
1631 : // Before this change, the twoDPB and twoDPBSq images were
1632 : // directly sent to applyPB(). Now, in the code below, we go
1633 : // over the poln. planes of these images, make a temp. image
1634 : // with the poln. planes of these images, applyPB() on the
1635 : // temp. image, and copy the result in the appropriate planes
1636 : // of the twoD* images. The rest of the code therefore does
1637 : // not see a difference.
1638 : {
1639 0 : CoordinateSystem singlePolCoords=coords;
1640 0 : Int index, outNdx, whichPolPlane=0;
1641 0 : Vector<Int> inStokes, outStokes(1);
1642 0 : index = coords.findCoordinate(Coordinate::STOKES);
1643 0 : inStokes = coords.stokesCoordinate(index).stokes();
1644 0 : for (uInt i=0; i<inStokes.nelements(); i++)
1645 : {
1646 0 : outStokes(0)=inStokes(i);
1647 : // Make a temp. image with a single Stokes along Pol. axis.
1648 0 : StokesCoordinate polnCoord(outStokes);
1649 0 : outNdx = singlePolCoords.findCoordinate(Coordinate::STOKES);
1650 0 : singlePolCoords.replaceCoordinate(polnCoord, outNdx);
1651 0 : IPosition singlePolShape = pbSqShp;
1652 0 : singlePolShape(2)=1;
1653 0 : Bool doSquint=true;
1654 0 : Double pa=getPA(vb);
1655 0 : Vector<Double> chanFreq = vb.frequency();
1656 : //freqHi = max(chanFreq);
1657 0 : Double Freq = max(chanFreq);
1658 : {
1659 0 : TempImage<Complex> singlePolImg(singlePolShape, singlePolCoords);
1660 : // Copy screen to the single pol. image. Apply PB to it. Copy
1661 : // the single pol. image plane to the twoDPB image.
1662 0 : doSquint=true;
1663 0 : PolnPlane(2)=0; singlePolImg.putSlice(screen, PolnPlane);
1664 0 : vlaPB.applyPB(singlePolImg, pa, doSquint,bandID_p, 0, Freq);
1665 0 : PolnPlane(2)=whichPolPlane; twoDPB.putSlice(singlePolImg.get(), PolnPlane);
1666 0 : }
1667 : {
1668 0 : TempImage<Complex> singlePolImg(singlePolShape, singlePolCoords);
1669 : // Copy screen to the single pol. image. Apply PB to it. Copy
1670 : // the single pol. image plane to the twoDPBSq image.
1671 0 : doSquint = false;
1672 0 : PolnPlane(2)=0; singlePolImg.putSlice(screen, PolnPlane);
1673 : //vlaPB.applyPBSq(twoDPBSq, vb, bandID_p, doSquint);
1674 0 : vlaPB.applyPB(singlePolImg, pa, doSquint, bandID_p, 0, Freq);
1675 0 : PolnPlane(2)=whichPolPlane; twoDPBSq.putSlice(singlePolImg.get(), PolnPlane);
1676 0 : }
1677 :
1678 0 : whichPolPlane++;
1679 0 : }
1680 0 : }
1681 : /*
1682 : // twoDPB.put(abs(twoDPB.get()));
1683 : // twoDPBSq.put(abs(twoDPBSq.get()));
1684 : */
1685 :
1686 : // {
1687 : // String name("twoDPB.before.im");
1688 : // storeImg(name,twoDPB);
1689 : // }
1690 : // {
1691 : // //
1692 : // // Apply (multiply) by the Spheroidal functions
1693 : // //
1694 : // Vector<Float> maxVal(twoDPB.shape()(2));
1695 : // Vector<IPosition> posMax(twoDPB.shape()(2));
1696 :
1697 : // SynthesisUtils::findLatticeMax(twoDPB,maxVal,posMax);
1698 : // posMax(0)(0)+=1;
1699 : // posMax(0)(1)+=1;
1700 : // posMax(1)(0)+=1;
1701 : // posMax(1)(1)+=1;
1702 : // // applyAntiAliasingOp(twoDPB,posMax,1);
1703 :
1704 : // SynthesisUtils::findLatticeMax(twoDPBSq,maxVal,posMax);
1705 : // posMax(0)(0)+=1;
1706 : // posMax(0)(1)+=1;
1707 : // posMax(1)(0)+=1;
1708 : // posMax(1)(1)+=1;
1709 : // // applyAntiAliasingOp(twoDPBSq,posMax,1,true);
1710 : // }
1711 :
1712 0 : Complex cpeak=max(twoDPB.get());
1713 0 : twoDPB.put(twoDPB.get()/cpeak);
1714 0 : cpeak=max(twoDPBSq.get());
1715 0 : twoDPBSq.put(twoDPBSq.get()/cpeak);
1716 : // twoDPBSq.set(1.0);
1717 : // {
1718 : // String name("twoDPB.im");
1719 : // storeImg(name,twoDPB);
1720 : // }
1721 :
1722 0 : CoordinateSystem cs=twoDPB.coordinates();
1723 0 : Int index= twoDPB.coordinates().findCoordinate(Coordinate::SPECTRAL);
1724 0 : SpectralCoordinate SpCS = twoDPB.coordinates().spectralCoordinate(index);
1725 :
1726 0 : cfRefFreq_p=SpCS.referenceValue()(0);
1727 : //
1728 : // Now FFT and get the result back
1729 : //
1730 : // {
1731 : // String name("twoDPB.im");
1732 : // storeImg(name,twoDPB);
1733 : // }
1734 0 : LatticeFFT::cfft2d(twoDPB);
1735 0 : LatticeFFT::cfft2d(twoDPBSq);
1736 : // {
1737 : // String name("twoDPBFT.im");
1738 : // storeImg(name,twoDPB);
1739 : // }
1740 : //
1741 : // Fill the convolution function planes with the result.
1742 : //
1743 : {
1744 : // IPosition start(4, convSize/4, convSize/4, 0, 0),
1745 : // pbSlice(4, convSize/2-1, convSize/2-1, polInUse, 1);
1746 : // IPosition sliceStart(4,0,0,iw,0),
1747 : // sliceLength(4,convSize/2-1,convSize/2-1,1,polInUse);
1748 :
1749 0 : IPosition start(4, 0, 0, 0, 0),
1750 0 : pbSlice(4, twoDPB.shape()[0]-1, twoDPB.shape()[1]-1, polInUse, 1);
1751 0 : IPosition sliceStart(4,0,0,iw,0),
1752 0 : sliceLength(4,convFunc.shape()[0]-1,convFunc.shape()[1]-1,1,polInUse);
1753 :
1754 0 : convFunc(Slicer(sliceStart,sliceLength)).nonDegenerate()
1755 0 : =(twoDPB.getSlice(start, pbSlice, true));
1756 :
1757 0 : IPosition shp(twoDPBSq.shape());
1758 :
1759 : //UNUSED: Int bufSize=convWeights.shape()[0], Org=shp[0]/2;
1760 : // IPosition sqStart(4, Org-bufSize/2, Org-bufSize/2, 0, 0),
1761 : // pbSqSlice(4, bufSize-1, bufSize-1, polInUse, 1);
1762 : // IPosition sqSliceStart(4,0,0,iw,0),
1763 : // sqSliceLength(4,bufSize-1,bufSize-1,1,polInUse);
1764 :
1765 0 : IPosition sqStart(4, 0, 0, 0, 0),
1766 0 : pbSqSlice(4, shp[0]-1, shp[1]-1, polInUse, 1);
1767 0 : IPosition sqSliceStart(4,0,0,iw,0),
1768 0 : sqSliceLength(4,shp[0]-1,shp[1]-1,1,polInUse);
1769 :
1770 0 : convWeights(Slicer(sqSliceStart,sqSliceLength)).nonDegenerate()
1771 0 : =(twoDPBSq.getSlice(sqStart, pbSqSlice, true));
1772 :
1773 0 : }
1774 0 : }
1775 :
1776 : {
1777 0 : Complex cpeak = max(convFunc);
1778 0 : convFunc/=cpeak;
1779 : // cout << "#### max(convFunc) = " << cpeak << endl;
1780 0 : cpeak=max(convWeights);
1781 : // cout << "#### max(convWeights) = " << cpeak << endl;
1782 0 : convWeights/=cpeak;
1783 : // cout << "#### max(convWeights) = " << max(convWeights) << endl;
1784 : }
1785 : //
1786 : // Find the convolution function support size. No assumption
1787 : // about the symmetry of the conv. func. can be made (except that
1788 : // they are same for all poln. planes).
1789 : //
1790 0 : Int lastPASlot=PAIndex;
1791 0 : for(Int iw=0;iw<wConvSize;iw++)
1792 0 : for(Int ipol=0;ipol<polInUse;ipol++)
1793 0 : convSupport(iw,ipol,lastPASlot)=-1;
1794 : //
1795 : // !!!Caution: This assumes that the support size at each Poln. is
1796 : // the same, starting from the center pixel (in pixel
1797 : // co-ordinates). For large pointing offsets, this might not be
1798 : // true.
1799 : //
1800 : // Int ConvFuncOrigin=convSize/4; // Conv. Func. is half that size of convSize
1801 0 : Int ConvFuncOrigin=convFunc.shape()[0]/2; // Conv. Func. is half that size of convSize
1802 0 : IPosition ndx(4,ConvFuncOrigin,0,0,0);
1803 : // Cube<Int> convWtSupport(convSupport.shape());
1804 0 : convWtSupport.resize(convSupport.shape(),true);
1805 0 : Int maxConvWtSupport=0, supportBuffer;
1806 0 : for (Int iw=0;iw<wConvSize;iw++)
1807 : {
1808 0 : Bool found=false;
1809 : Float threshold;
1810 : Int R;
1811 0 : ndx(2) = iw;
1812 :
1813 0 : ndx(0)=ndx(1)=ConvFuncOrigin;
1814 0 : ndx(2) = iw;
1815 : // Complex maxVal = max(convFunc);
1816 0 : threshold = abs(convFunc(ndx))*THRESHOLD;
1817 : //
1818 : // Find the support size of the conv. function in pixels
1819 : //
1820 : Int wtR;
1821 0 : found =findSupport(convWeights,threshold,ConvFuncOrigin,wtR);
1822 0 : found = findSupport(convFunc,threshold,ConvFuncOrigin,R);
1823 :
1824 : // R *=2.5;
1825 : //
1826 : // Set the support size for each W-plane and for all
1827 : // Pol-planes. Assuming that support size for all Pol-planes
1828 : // is same.
1829 : //
1830 0 : if(found)
1831 : {
1832 : // Int maxR=R;//max(ndx(0),ndx(1));
1833 0 : for(Int ipol=0;ipol<polInUse;ipol++)
1834 : {
1835 0 : convSupport(iw,ipol,lastPASlot)=Int(R/Float(convSampling));
1836 0 : convSupport(iw,ipol,lastPASlot)=Int(0.5+Float(R)/Float(convSampling))+1;
1837 : // convSupport(iw,ipol,lastPASlot) += (convSupport(iw,ipol,lastPASlot)+1)%2;
1838 0 : convWtSupport(iw,ipol,lastPASlot)=Int(R*CONVWTSIZEFACTOR/Float(convSampling));
1839 0 : convWtSupport(iw,ipol,lastPASlot)=Int(0.5+Float(R)*CONVWTSIZEFACTOR/Float(convSampling))+1;
1840 : // convWtSupport(iw,ipol,lastPASlot)=Int(wtR/Float(convSampling));
1841 : // convWtSupport(iw,ipol,lastPASlot) += (convWtSupport(iw,ipol,lastPASlot)+1)%2;
1842 0 : if ((lastPASlot == 0) || (maxConvSupport == -1))
1843 0 : if (convSupport(iw,ipol,lastPASlot) > maxConvSupport)
1844 0 : maxConvSupport = convSupport(iw,ipol,lastPASlot);
1845 0 : maxConvWtSupport=convWtSupport(iw,ipol,lastPASlot);
1846 : }
1847 : }
1848 : }
1849 :
1850 0 : if(convSupport(0,0,lastPASlot)<1)
1851 : logIO() << "Convolution function is misbehaved - support seems to be zero"
1852 0 : << LogIO::EXCEPTION;
1853 :
1854 0 : logIO() << LogOrigin("nPBWProjectFT", "makeConvFunction")
1855 : << "Re-sizing the convolution functions"
1856 0 : << LogIO::POST;
1857 :
1858 : {
1859 0 : supportBuffer = OVERSAMPLING;
1860 0 : Int bot=ConvFuncOrigin-convSampling*maxConvSupport-supportBuffer,//-convSampling/2,
1861 0 : top=ConvFuncOrigin+convSampling*maxConvSupport+supportBuffer;//+convSampling/2;
1862 0 : bot = max(0,bot);
1863 0 : top = min(top, convFunc.shape()(0)-1);
1864 : {
1865 0 : Array<Complex> tmp;
1866 0 : IPosition blc(4,bot,bot,0,0), trc(4,top,top,wConvSize-1,polInUse-1);
1867 :
1868 0 : tmp = convFunc(blc,trc);
1869 0 : (*convFuncCache[lastPASlot]).resize(tmp.shape());
1870 0 : (*convFuncCache[lastPASlot]) = tmp;
1871 0 : convFunc.reference(*convFuncCache[lastPASlot]);
1872 0 : }
1873 :
1874 0 : supportBuffer = (Int)(OVERSAMPLING*CONVWTSIZEFACTOR);
1875 0 : bot=ConvFuncOrigin-convSampling*maxConvWtSupport-supportBuffer;
1876 0 : top=ConvFuncOrigin+convSampling*maxConvWtSupport+supportBuffer;
1877 0 : bot=max(0,bot);
1878 0 : top=min(top,convWeights.shape()(0)-1);
1879 : {
1880 0 : Array<Complex> tmp;
1881 0 : IPosition blc(4,bot,bot,0,0), trc(4,top,top,wConvSize-1,polInUse-1);
1882 :
1883 0 : tmp = convWeights(blc,trc);
1884 0 : (*convWeightsCache[lastPASlot]).resize(tmp.shape());
1885 0 : (*convWeightsCache[lastPASlot]) = tmp;
1886 0 : convWeights.reference(*convWeightsCache[lastPASlot]);
1887 0 : }
1888 : }
1889 :
1890 : //
1891 : // Normalize such that plane 0 sums to 1 (when jumping in steps of
1892 : // convSampling). This is really not necessary here since we do
1893 : // the normalizing by the area more accurately in the gridder
1894 : // (fpbwproj.f).
1895 : //
1896 0 : ndx(2)=ndx(3)=0;
1897 :
1898 :
1899 0 : Complex pbSum=0.0;
1900 0 : IPosition peakPix(ndx);
1901 :
1902 0 : Int Nx = convFunc.shape()(0), Ny=convFunc.shape()(1);
1903 :
1904 0 : for(Int nw=0;nw<wConvSize;nw++)
1905 0 : for(Int np=0;np<polInUse;np++)
1906 : {
1907 0 : ndx(2) = nw; ndx(3)=np;
1908 : {
1909 : //
1910 : // Locate the pixel with the peak value. That's the
1911 : // origin in pixel co-ordinates.
1912 : //
1913 0 : Float peak=0;
1914 0 : peakPix = 0;
1915 0 : for(ndx(1)=0;ndx(1)<convFunc.shape()(1);ndx(1)++)
1916 0 : for(ndx(0)=0;ndx(0)<convFunc.shape()(0);ndx(0)++)
1917 0 : if (abs(convFunc(ndx)) > peak) {peakPix = ndx;peak=abs(convFunc(ndx));}
1918 : }
1919 :
1920 0 : ConvFuncOrigin = peakPix(0);
1921 : // ConvFuncOrigin = convFunc.shape()(0)/2+1;
1922 : // Int thisConvSupport=convSampling*convSupport(nw,np,lastPASlot);
1923 0 : Int thisConvSupport=convSupport(nw,np,lastPASlot);
1924 0 : pbSum=0.0;
1925 :
1926 0 : for(Int iy=-thisConvSupport;iy<thisConvSupport;iy++)
1927 0 : for(Int ix=-thisConvSupport;ix<thisConvSupport;ix++)
1928 : {
1929 0 : ndx(0)=ix*convSampling+ConvFuncOrigin;
1930 0 : ndx(1)=iy*convSampling+ConvFuncOrigin;
1931 0 : pbSum += real(convFunc(ndx));
1932 : }
1933 : /*
1934 : for(Int iy=0;iy<Ny;iy++)
1935 : for(Int ix=0;ix<Nx;ix++)
1936 : {
1937 : ndx(0)=ix;ndx(1)=iy;
1938 : pbSum += convFunc(ndx);
1939 : }
1940 : */
1941 0 : if(pbSum>0.0)
1942 : {
1943 : //
1944 : // Normalize each Poln. plane by the area under its convfunc.
1945 : //
1946 0 : Nx = convFunc.shape()(0), Ny = convFunc.shape()(1);
1947 0 : for (ndx(1)=0;ndx(1)<Ny;ndx(1)++)
1948 0 : for (ndx(0)=0;ndx(0)<Nx;ndx(0)++)
1949 : {
1950 0 : convFunc(ndx) /= pbSum;
1951 : }
1952 :
1953 0 : Nx = convWeights.shape()(0); Ny = convWeights.shape()(1);
1954 0 : for (ndx(1)=0; ndx(1)<Ny; ndx(1)++)
1955 0 : for (ndx(0)=0; ndx(0)<Nx; ndx(0)++)
1956 : {
1957 0 : convWeights(ndx) /= pbSum*pbSum;
1958 : // if ((ndx(0)==Nx/2+1) && (ndx(1)==Ny/2+1))
1959 : // {
1960 : // convWeights(ndx)=1.0;
1961 : // cout << ndx << " " << convWeights(ndx) << endl;
1962 : // }
1963 : // else
1964 : // convWeights(ndx)=0.0;
1965 : }
1966 : }
1967 : else
1968 0 : throw(SynthesisFTMachineError("Convolution function integral is not positive"));
1969 :
1970 0 : Vector<Float> maxVal(convWeights.shape()(2));
1971 0 : Vector<IPosition> posMax(convWeights.shape()(2));
1972 0 : SynthesisUtils::findLatticeMax(convWeights,maxVal,posMax);
1973 : // cout << "convWeights: " << maxVal << " " << posMax << endl;
1974 0 : }
1975 :
1976 0 : }
1977 : //
1978 : //------------------------------------------------------------------------------
1979 : //
1980 0 : void nPBWProjectFT::initializeToVis(ImageInterface<Complex>& iimage,
1981 : const VisBuffer& vb)
1982 : {
1983 0 : image=&iimage;
1984 :
1985 0 : ok();
1986 :
1987 0 : init();
1988 0 : makingPSF = false;
1989 0 : initMaps(vb);
1990 :
1991 0 : findConvFunction(*image, vb);
1992 : //
1993 : // Initialize the maps for polarization and channel. These maps
1994 : // translate visibility indices into image indices
1995 : //
1996 :
1997 0 : nx = image->shape()(0);
1998 0 : ny = image->shape()(1);
1999 0 : npol = image->shape()(2);
2000 0 : nchan = image->shape()(3);
2001 :
2002 0 : if(image->shape().product()>cachesize) isTiled=true;
2003 0 : else isTiled=false;
2004 : //
2005 : // If we are memory-based then read the image in and create an
2006 : // ArrayLattice otherwise just use the PagedImage
2007 : //
2008 :
2009 0 : isTiled=false;
2010 :
2011 0 : if(isTiled){
2012 0 : lattice=CountedPtr<Lattice<Complex> > (image, false);
2013 : }
2014 : else
2015 : {
2016 0 : IPosition gridShape(4, nx, ny, npol, nchan);
2017 0 : griddedData.resize(gridShape);
2018 0 : griddedData=Complex(0.0);
2019 :
2020 0 : IPosition stride(4, 1);
2021 0 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
2022 0 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
2023 0 : IPosition trc(blc+image->shape()-stride);
2024 :
2025 0 : IPosition start(4, 0);
2026 0 : griddedData(blc, trc) = image->getSlice(start, image->shape());
2027 :
2028 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
2029 0 : lattice=arrayLattice;
2030 0 : }
2031 :
2032 : //AlwaysAssert(lattice, AipsError);
2033 :
2034 0 : logIO() << LogIO::DEBUGGING << "Starting FFT of image" << LogIO::POST;
2035 :
2036 0 : Vector<Float> sincConv(nx);
2037 0 : Float centerX=nx/2;
2038 0 : for (Int ix=0;ix<nx;ix++)
2039 : {
2040 0 : Float x=C::pi*Float(ix-centerX)/(Float(nx)*Float(convSampling));
2041 0 : if(ix==centerX) sincConv(ix)=1.0;
2042 0 : else sincConv(ix)=sin(x)/x;
2043 : }
2044 :
2045 0 : Vector<Complex> correction(nx);
2046 : //
2047 : // Do the Grid-correction
2048 : //
2049 : {
2050 0 : normalizeAvgPB();
2051 :
2052 0 : IPosition cursorShape(4, nx, 1, 1, 1);
2053 0 : IPosition axisPath(4, 0, 1, 2, 3);
2054 0 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
2055 0 : LatticeIterator<Complex> lix(*lattice, lsx);
2056 :
2057 0 : verifyShapes(avgPB.shape(), image->shape());
2058 0 : Array<Float> avgBuf; avgPB.get(avgBuf);
2059 0 : if (max(avgBuf) < 1e-04)
2060 0 : throw(AipsError("Normalization by PB requested but either PB not found in the cache "
2061 0 : "or is ill-formed."));
2062 :
2063 0 : LatticeStepper lpb(avgPB.shape(),cursorShape,axisPath);
2064 0 : LatticeIterator<Float> lipb(avgPB, lpb);
2065 :
2066 0 : Vector<Complex> griddedVis;
2067 : //
2068 : // Grid correct in anticipation of the convolution by the
2069 : // convFunc. Each polarization plane is corrected by the
2070 : // appropraite primary beam.
2071 : //
2072 0 : for(lix.reset(),lipb.reset();!lix.atEnd();lix++,lipb++)
2073 : {
2074 0 : Int iy=lix.position()(1);
2075 0 : gridder->correctX1D(correction,iy);
2076 0 : griddedVis = lix.rwVectorCursor();
2077 :
2078 0 : Vector<Float> PBCorrection(lipb.rwVectorCursor().shape());
2079 0 : PBCorrection = lipb.rwVectorCursor();
2080 0 : for(int ix=0;ix<nx;ix++)
2081 : {
2082 : // PBCorrection(ix) = (FUNC(PBCorrection(ix)))/(sincConv(ix)*sincConv(iy));
2083 :
2084 : //
2085 : // This is with PS functions included
2086 : //
2087 : // if (doPBCorrection)
2088 : // {
2089 : // PBCorrection(ix) = FUNC(PBCorrection(ix))/(sincConv(ix)*sincConv(iy));
2090 : // //PBCorrection(ix) = FUNC(PBCorrection(ix))*(sincConv(ix)*sincConv(iy));
2091 : // if ((abs(PBCorrection(ix)*correction(ix))) >= pbLimit_p)
2092 : // {lix.rwVectorCursor()(ix) /= (PBCorrection(ix))*correction(ix);}
2093 : // else
2094 : // {lix.rwVectorCursor()(ix) *= (sincConv(ix)*sincConv(iy));}
2095 : // }
2096 : // else
2097 : // lix.rwVectorCursor()(ix) /= (correction(ix)/(sincConv(ix)*sincConv(iy)));
2098 : //
2099 : // This without the PS functions
2100 : //
2101 0 : if (doPBCorrection)
2102 : {
2103 : // PBCorrection(ix) = FUNC(PBCorrection(ix))/(sincConv(ix)*sincConv(iy));
2104 0 : PBCorrection(ix) = FUNC(PBCorrection(ix))*(sincConv(ix)*sincConv(iy));
2105 : // PBCorrection(ix) = (PBCorrection(ix))*(sincConv(ix)*sincConv(iy));
2106 0 : if ((abs(PBCorrection(ix))) >= pbLimit_p)
2107 0 : {lix.rwVectorCursor()(ix) /= (PBCorrection(ix));}
2108 : else
2109 0 : {lix.rwVectorCursor()(ix) *= (sincConv(ix)*sincConv(iy));}
2110 : }
2111 : else
2112 0 : lix.rwVectorCursor()(ix) /= (1.0/(sincConv(ix)*sincConv(iy)));
2113 : }
2114 0 : }
2115 0 : }
2116 : // {
2117 : // ostringstream name;
2118 : // cout << image->shape() << endl;
2119 : // name << "theModel.im";
2120 : // PagedImage<Float> tmp(image->shape(), image->coordinates(), name);
2121 : // Array<Complex> buf;
2122 : // Bool isRef = lattice->get(buf);
2123 : // cout << "The model max. = " << max(buf) << endl;
2124 : // LatticeExpr<Float> le(abs((*lattice)));
2125 : // tmp.copyData(le);
2126 : // }
2127 : //
2128 : // Now do the FFT2D in place
2129 : //
2130 : // {
2131 : // Array<Complex> buf;
2132 : // Bool isRef = lattice->get(buf);
2133 : // }
2134 0 : LatticeFFT::cfft2d(*lattice);
2135 :
2136 0 : logIO() << LogIO::DEBUGGING << "Finished FFT" << LogIO::POST;
2137 0 : }
2138 : //
2139 : //---------------------------------------------------------------
2140 : //
2141 0 : void nPBWProjectFT::initializeToVis(ImageInterface<Complex>& iimage,
2142 : const VisBuffer& vb,
2143 : Array<Complex>& griddedVis,
2144 : Vector<Double>& uvscale)
2145 : {
2146 0 : initializeToVis(iimage, vb);
2147 0 : griddedVis.assign(griddedData); //using the copy for storage
2148 0 : uvscale.assign(uvScale);
2149 0 : }
2150 : //
2151 : //---------------------------------------------------------------
2152 : //
2153 0 : void nPBWProjectFT::finalizeToVis()
2154 : {
2155 0 : logIO() << "##########finalizeToVis()###########" << LogIO::DEBUGGING << LogIO::POST;
2156 0 : if(isTiled)
2157 : {
2158 0 : logIO() << LogOrigin("nPBWProjectFT", "finalizeToVis") << LogIO::NORMAL;
2159 :
2160 0 : AlwaysAssert(imageCache, AipsError);
2161 0 : AlwaysAssert(image, AipsError);
2162 0 : ostringstream o;
2163 0 : imageCache->flush();
2164 0 : imageCache->showCacheStatistics(o);
2165 0 : logIO() << o.str() << LogIO::POST;
2166 0 : }
2167 0 : if(pointingToImage) delete pointingToImage; pointingToImage=0;
2168 0 : }
2169 : //
2170 : //---------------------------------------------------------------
2171 : //
2172 : // Initialize the FFT to the Sky. Here we have to setup and
2173 : // initialize the grid.
2174 : //
2175 0 : void nPBWProjectFT::initializeToSky(ImageInterface<Complex>& iimage,
2176 : Matrix<Float>& weight,
2177 : const VisBuffer& vb)
2178 : {
2179 0 : logIO() << "#########initializeToSky()##########" << LogIO::DEBUGGING << LogIO::POST;
2180 :
2181 : // image always points to the image
2182 0 : image=&iimage;
2183 :
2184 0 : init();
2185 0 : initMaps(vb);
2186 :
2187 : // Initialize the maps for polarization and channel. These maps
2188 : // translate visibility indices into image indices
2189 :
2190 0 : nx = image->shape()(0);
2191 0 : ny = image->shape()(1);
2192 0 : npol = image->shape()(2);
2193 0 : nchan = image->shape()(3);
2194 :
2195 0 : if(image->shape().product()>cachesize) isTiled=true;
2196 0 : else isTiled=false;
2197 :
2198 :
2199 0 : sumWeight=0.0;
2200 0 : weight.resize(sumWeight.shape());
2201 0 : weight=0.0;
2202 : //
2203 : // Initialize for in memory or to disk gridding. lattice will
2204 : // point to the appropriate Lattice, either the ArrayLattice for
2205 : // in memory gridding or to the image for to disk gridding.
2206 : //
2207 0 : if(isTiled)
2208 : {
2209 0 : imageCache->flush();
2210 0 : image->set(Complex(0.0));
2211 0 : lattice=CountedPtr<Lattice<Complex> > (image, false);
2212 : }
2213 : else
2214 : {
2215 0 : IPosition gridShape(4, nx, ny, npol, nchan);
2216 0 : griddedData.resize(gridShape);
2217 0 : griddedData=Complex(0.0);
2218 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
2219 0 : lattice=arrayLattice;
2220 0 : }
2221 0 : PAIndex = -1;
2222 0 : }
2223 : //
2224 : //---------------------------------------------------------------
2225 : //
2226 0 : void nPBWProjectFT::finalizeToSky()
2227 : {
2228 : //
2229 : // Now we flush the cache and report statistics For memory based,
2230 : // we don't write anything out yet.
2231 : //
2232 0 : logIO() << "#########finalizeToSky()#########" << LogIO::DEBUGGING << LogIO::POST;
2233 0 : if(isTiled)
2234 : {
2235 0 : logIO() << LogOrigin("nPBWProjectFT", "finalizeToSky") << LogIO::NORMAL;
2236 :
2237 0 : AlwaysAssert(image, AipsError);
2238 0 : AlwaysAssert(imageCache, AipsError);
2239 0 : imageCache->flush();
2240 0 : ostringstream o;
2241 0 : imageCache->showCacheStatistics(o);
2242 0 : logIO() << o.str() << LogIO::POST;
2243 0 : }
2244 0 : if(pointingToImage) delete pointingToImage; pointingToImage=0;
2245 0 : PAIndex = -1;
2246 :
2247 0 : paChangeDetector.reset();
2248 0 : cfCache.finalize();
2249 0 : convFuncCacheReady=true;
2250 0 : }
2251 : //
2252 : //---------------------------------------------------------------
2253 : //
2254 0 : Array<Complex>* nPBWProjectFT::getDataPointer(const IPosition& centerLoc2D,
2255 : Bool readonly)
2256 : {
2257 : Array<Complex>* result;
2258 : // Is tiled: get tiles and set up offsets
2259 0 : centerLoc(0)=centerLoc2D(0);
2260 0 : centerLoc(1)=centerLoc2D(1);
2261 0 : result=&imageCache->tile(offsetLoc, centerLoc, readonly);
2262 0 : gridder->setOffset(IPosition(2, offsetLoc(0), offsetLoc(1)));
2263 0 : return result;
2264 : }
2265 :
2266 : #define NEED_UNDERSCORES
2267 : #if defined(NEED_UNDERSCORES)
2268 : #define gpbwproj gpbwproj_
2269 : #define dpbwproj dpbwproj_
2270 : #define dpbwgrad dpbwgrad_
2271 : #endif
2272 :
2273 : extern "C" {
2274 : void gpbwproj(Double *uvw,
2275 : Double *dphase,
2276 : const Complex *values,
2277 : Int *nvispol,
2278 : Int *nvischan,
2279 : Int *dopsf,
2280 : const Int *flag,
2281 : const Int *rflag,
2282 : const Float *weight,
2283 : Int *nrow,
2284 : Int *rownum,
2285 : Double *scale,
2286 : Double *offset,
2287 : Complex *grid,
2288 : Int *nx,
2289 : Int *ny,
2290 : Int *npol,
2291 : Int *nchan,
2292 : const Double *freq,
2293 : const Double *c,
2294 : Int *support,
2295 : Int *convsize,
2296 : Int *sampling,
2297 : Int *wconvsize,
2298 : Complex *convfunc,
2299 : Int *chanmap,
2300 : Int *polmap,
2301 : Int *polused,
2302 : Double *sumwt,
2303 : Int *ant1,
2304 : Int *ant2,
2305 : Int *nant,
2306 : Int *scanno,
2307 : Double *sigma,
2308 : Float *raoff,
2309 : Float *decoff,
2310 : Double *area,
2311 : Int *doGrad,
2312 : Int *doPointingCorrection,
2313 : Int *nPA,
2314 : Int *paIndex,
2315 : Int *CFMap,
2316 : Int *ConjCFMap,
2317 : Double *currentCFPA, Double *actualPA,Double *cfRefFreq_p);
2318 : void dpbwproj(Double *uvw,
2319 : Double *dphase,
2320 : Complex *values,
2321 : Int *nvispol,
2322 : Int *nvischan,
2323 : const Int *flag,
2324 : const Int *rflag,
2325 : Int *nrow,
2326 : Int *rownum,
2327 : Double *scale,
2328 : Double *offset,
2329 : const Complex *grid,
2330 : Int *nx,
2331 : Int *ny,
2332 : Int *npol,
2333 : Int *nchan,
2334 : const Double *freq,
2335 : const Double *c,
2336 : Int *support,
2337 : Int *convsize,
2338 : Int *sampling,
2339 : Int *wconvsize,
2340 : Complex *convfunc,
2341 : Int *chanmap,
2342 : Int *polmap,
2343 : Int *polused,
2344 : Int *ant1,
2345 : Int *ant2,
2346 : Int *nant,
2347 : Int *scanno,
2348 : Double *sigma,
2349 : Float *raoff, Float *decoff,
2350 : Double *area,
2351 : Int *dograd,
2352 : Int *doPointingCorrection,
2353 : Int *nPA,
2354 : Int *paIndex,
2355 : Int *CFMap,
2356 : Int *ConjCFMap,
2357 : Double *currentCFPA, Double *actualPA, Double *cfRefFreq_p);
2358 : void dpbwgrad(Double *uvw,
2359 : Double *dphase,
2360 : Complex *values,
2361 : Int *nvispol,
2362 : Int *nvischan,
2363 : Complex *gazvalues,
2364 : Complex *gelvalues,
2365 : Int *doconj,
2366 : const Int *flag,
2367 : const Int *rflag,
2368 : Int *nrow,
2369 : Int *rownum,
2370 : Double *scale,
2371 : Double *offset,
2372 : const Complex *grid,
2373 : Int *nx,
2374 : Int *ny,
2375 : Int *npol,
2376 : Int *nchan,
2377 : const Double *freq,
2378 : const Double *c,
2379 : Int *support,
2380 : Int *convsize,
2381 : Int *sampling,
2382 : Int *wconvsize,
2383 : Complex *convfunc,
2384 : Int *chanmap,
2385 : Int *polmap,
2386 : Int *polused,
2387 : Int *ant1,
2388 : Int *ant2,
2389 : Int *nant,
2390 : Int *scanno,
2391 : Double *sigma,
2392 : Float *raoff, Float *decoff,
2393 : Double *area,
2394 : Int *dograd,
2395 : Int *doPointingCorrection,
2396 : Int *nPA,
2397 : Int *paIndex,
2398 : Int *CFMap,
2399 : Int *ConjCFMap,
2400 : Double *currentCFPA, Double *actualPA, Double *cfRefFreq_p);
2401 : }
2402 : //
2403 : //----------------------------------------------------------------------
2404 : //
2405 0 : void nPBWProjectFT::runFortranGet(Matrix<Double>& uvw,Vector<Double>& dphase,
2406 : Cube<Complex>& visdata,
2407 : IPosition& s,
2408 : // Cube<Complex>& gradVisAzData,
2409 : // Cube<Complex>& gradVisElData,
2410 : // IPosition& gradS,
2411 : Int& /*Conj*/,
2412 : Cube<Int>& flags,Vector<Int>& rowFlags,
2413 : Int& rownr,Vector<Double>& actualOffset,
2414 : Array<Complex>* dataPtr,
2415 : Int& aNx, Int& aNy, Int& npol, Int& nchan,
2416 : VisBuffer& vb,Int& Nant_p, Int& scanNo,
2417 : Double& sigma,
2418 : Array<Float>& l_off,
2419 : Array<Float>& m_off,
2420 : Double area,
2421 : Int& doGrad,
2422 : Int paIndex)
2423 : {
2424 : enum whichGetStorage {RAOFF,DECOFF,UVW,DPHASE,VISDATA,GRADVISAZ,GRADVISEL,
2425 : FLAGS,ROWFLAGS,UVSCALE,ACTUALOFFSET,DATAPTR,VBFREQ,
2426 : CONVSUPPORT,CONVFUNC,CHANMAP,POLMAP,VBANT1,VBANT2,CONJCFMAP,CFMAP};
2427 0 : Vector<Bool> deleteThem(21);
2428 :
2429 : Double *uvw_p, *dphase_p, *actualOffset_p, *vb_freq_p, *uvScale_p;
2430 : Complex *visdata_p, *dataPtr_p, *f_convFunc_p;
2431 : Int *flags_p, *rowFlags_p, *chanMap_p, *polMap_p, *convSupport_p, *vb_ant1_p, *vb_ant2_p,
2432 : *ConjCFMap_p, *CFMap_p;
2433 : Float *l_off_p, *m_off_p;
2434 : Double actualPA;
2435 :
2436 0 : Vector<Int> ConjCFMap, CFMap;
2437 : /*
2438 : ConjCFMap = CFMap = polMap;
2439 : CFMap = makeConjPolMap(vb);
2440 : */
2441 : Int N;
2442 0 : actualPA = getVBPA(vb);
2443 :
2444 0 : N=polMap.nelements();
2445 0 : CFMap = polMap; ConjCFMap = polMap;
2446 0 : for(Int i=0;i<N;i++) CFMap[i] = polMap[N-i-1];
2447 :
2448 0 : Array<Complex> rotatedConvFunc;
2449 : // SynthesisUtils::rotateComplexArray(logIO(), convFunc, convFuncCS_p,
2450 : // rotatedConvFunc,(currentCFPA-actualPA),"CUBIC");
2451 0 : SynthesisUtils::rotateComplexArray(logIO(), convFunc, convFuncCS_p,
2452 : rotatedConvFunc,0.0,"LINEAR");
2453 : // SynthesisUtils::rotateComplexArray(logIO(), convFunc, convFuncCS_p,
2454 : // rotatedConvFunc,(currentCFPA-actualPA),"LINEAR");
2455 :
2456 0 : ConjCFMap = polMap;
2457 0 : makeCFPolMap(vb,CFMap);
2458 0 : makeConjPolMap(vb,CFMap,ConjCFMap);
2459 :
2460 :
2461 0 : ConjCFMap_p = ConjCFMap.getStorage(deleteThem(CONJCFMAP));
2462 0 : CFMap_p = CFMap.getStorage(deleteThem(CFMAP));
2463 :
2464 0 : uvw_p = uvw.getStorage(deleteThem(UVW));
2465 0 : dphase_p = dphase.getStorage(deleteThem(DPHASE));
2466 0 : visdata_p = visdata.getStorage(deleteThem(VISDATA));
2467 : // gradVisAzData_p = gradVisAzData.getStorage(deleteThem(GRADVISAZ));
2468 : // gradVisElData_p = gradVisElData.getStorage(deleteThem(GRADVISEL));
2469 0 : flags_p = flags.getStorage(deleteThem(FLAGS));
2470 0 : rowFlags_p = rowFlags.getStorage(deleteThem(ROWFLAGS));
2471 0 : uvScale_p = uvScale.getStorage(deleteThem(UVSCALE));
2472 0 : actualOffset_p = actualOffset.getStorage(deleteThem(ACTUALOFFSET));
2473 0 : dataPtr_p = dataPtr->getStorage(deleteThem(DATAPTR));
2474 0 : vb_freq_p = vb.frequency().getStorage(deleteThem(VBFREQ));
2475 0 : convSupport_p = convSupport.getStorage(deleteThem(CONVSUPPORT));
2476 : // f_convFunc_p = convFunc.getStorage(deleteThem(CONVFUNC));
2477 0 : f_convFunc_p = rotatedConvFunc.getStorage(deleteThem(CONVFUNC));
2478 0 : chanMap_p = chanMap.getStorage(deleteThem(CHANMAP));
2479 0 : polMap_p = polMap.getStorage(deleteThem(POLMAP));
2480 0 : vb_ant1_p = vb.antenna1().getStorage(deleteThem(VBANT1));
2481 0 : vb_ant2_p = vb.antenna2().getStorage(deleteThem(VBANT2));
2482 0 : l_off_p = l_off.getStorage(deleteThem(RAOFF));
2483 0 : m_off_p = m_off.getStorage(deleteThem(DECOFF));
2484 :
2485 0 : Int npa=convSupport.shape()(2),actualConvSize;
2486 0 : Int paIndex_Fortran = paIndex;
2487 0 : actualConvSize = convFunc.shape()(0);
2488 :
2489 0 : IPosition shp=convSupport.shape();
2490 :
2491 0 : dpbwproj(uvw_p,
2492 : dphase_p,
2493 : // vb.modelVisCube().getStorage(del),
2494 : visdata_p,
2495 0 : &s.asVector()(0),
2496 0 : &s.asVector()(1),
2497 : // gradVisAzData_p,
2498 : // gradVisElData_p,
2499 : // &gradS(0),
2500 : // &gradS(1),
2501 : // &Conj,
2502 : flags_p,
2503 : rowFlags_p,
2504 0 : &s.asVector()(2),
2505 : &rownr,
2506 : uvScale_p,
2507 : actualOffset_p,
2508 : dataPtr_p,
2509 : &aNx,
2510 : &aNy,
2511 : &npol,
2512 : &nchan,
2513 : vb_freq_p,
2514 : &C::c,
2515 : convSupport_p,
2516 : &actualConvSize,
2517 : &convSampling,
2518 : &wConvSize,
2519 : f_convFunc_p,
2520 : chanMap_p,
2521 : polMap_p,
2522 : &polInUse,
2523 : vb_ant1_p,
2524 : vb_ant2_p,
2525 : &Nant_p,
2526 : &scanNo,
2527 : &sigma,
2528 : l_off_p, m_off_p,
2529 : &area,
2530 : &doGrad,
2531 : &doPointing,
2532 : &npa,
2533 : &paIndex_Fortran,
2534 : CFMap_p,
2535 : ConjCFMap_p,
2536 : ¤tCFPA
2537 : ,&actualPA,&cfRefFreq_p
2538 : );
2539 :
2540 0 : ConjCFMap.freeStorage((const Int *&)ConjCFMap_p,deleteThem(CONJCFMAP));
2541 0 : CFMap.freeStorage((const Int *&)CFMap_p,deleteThem(CFMAP));
2542 :
2543 0 : l_off.freeStorage((const Float*&)l_off_p,deleteThem(RAOFF));
2544 0 : m_off.freeStorage((const Float*&)m_off_p,deleteThem(DECOFF));
2545 0 : uvw.freeStorage((const Double*&)uvw_p,deleteThem(UVW));
2546 0 : dphase.freeStorage((const Double*&)dphase_p,deleteThem(DPHASE));
2547 0 : visdata.putStorage(visdata_p,deleteThem(VISDATA));
2548 0 : flags.freeStorage((const Int*&) flags_p,deleteThem(FLAGS));
2549 0 : rowFlags.freeStorage((const Int *&)rowFlags_p,deleteThem(ROWFLAGS));
2550 0 : actualOffset.freeStorage((const Double*&)actualOffset_p,deleteThem(ACTUALOFFSET));
2551 0 : dataPtr->freeStorage((const Complex *&)dataPtr_p,deleteThem(DATAPTR));
2552 0 : uvScale.freeStorage((const Double*&) uvScale_p,deleteThem(UVSCALE));
2553 0 : vb.frequency().freeStorage((const Double*&)vb_freq_p,deleteThem(VBFREQ));
2554 0 : convSupport.freeStorage((const Int*&)convSupport_p,deleteThem(CONVSUPPORT));
2555 0 : convFunc.freeStorage((const Complex *&)f_convFunc_p,deleteThem(CONVFUNC));
2556 0 : chanMap.freeStorage((const Int*&)chanMap_p,deleteThem(CHANMAP));
2557 0 : polMap.freeStorage((const Int*&) polMap_p,deleteThem(POLMAP));
2558 0 : vb.antenna1().freeStorage((const Int*&) vb_ant1_p,deleteThem(VBANT1));
2559 0 : vb.antenna2().freeStorage((const Int*&) vb_ant2_p,deleteThem(VBANT2));
2560 0 : }
2561 : //
2562 : //----------------------------------------------------------------------
2563 : //
2564 0 : void nPBWProjectFT::runFortranGetGrad(Matrix<Double>& uvw,Vector<Double>& dphase,
2565 : Cube<Complex>& visdata,
2566 : IPosition& s,
2567 : Cube<Complex>& gradVisAzData,
2568 : Cube<Complex>& gradVisElData,
2569 : // IPosition& gradS,
2570 : Int& Conj,
2571 : Cube<Int>& flags,Vector<Int>& rowFlags,
2572 : Int& rownr,Vector<Double>& actualOffset,
2573 : Array<Complex>* dataPtr,
2574 : Int& aNx, Int& aNy, Int& npol, Int& nchan,
2575 : VisBuffer& vb,Int& Nant_p, Int& scanNo,
2576 : Double& sigma,
2577 : Array<Float>& l_off,
2578 : Array<Float>& m_off,
2579 : Double area,
2580 : Int& doGrad,
2581 : Int paIndex)
2582 : {
2583 : enum whichGetStorage {RAOFF,DECOFF,UVW,DPHASE,VISDATA,GRADVISAZ,GRADVISEL,
2584 : FLAGS,ROWFLAGS,UVSCALE,ACTUALOFFSET,DATAPTR,VBFREQ,
2585 : CONVSUPPORT,CONVFUNC,CHANMAP,POLMAP,VBANT1,VBANT2,CONJCFMAP,CFMAP};
2586 0 : Vector<Bool> deleteThem(21);
2587 :
2588 : Double *uvw_p, *dphase_p, *actualOffset_p, *vb_freq_p, *uvScale_p;
2589 : Complex *visdata_p, *dataPtr_p, *f_convFunc_p;
2590 : Complex *gradVisAzData_p, *gradVisElData_p;
2591 : Int *flags_p, *rowFlags_p, *chanMap_p, *polMap_p, *convSupport_p, *vb_ant1_p, *vb_ant2_p,
2592 : *ConjCFMap_p, *CFMap_p;
2593 : Float *l_off_p, *m_off_p;
2594 : Double actualPA;
2595 :
2596 0 : Vector<Int> ConjCFMap, CFMap;
2597 0 : actualPA = getVBPA(vb);
2598 0 : ConjCFMap = polMap;
2599 0 : makeCFPolMap(vb,CFMap);
2600 0 : makeConjPolMap(vb,CFMap,ConjCFMap);
2601 :
2602 0 : Array<Complex> rotatedConvFunc;
2603 : // SynthesisUtils::rotateComplexArray(logIO(), convFunc, convFuncCS_p,
2604 : // rotatedConvFunc,(currentCFPA-actualPA),"LINEAR");
2605 0 : SynthesisUtils::rotateComplexArray(logIO(), convFunc, convFuncCS_p,
2606 : rotatedConvFunc,0.0);
2607 : // SynthesisUtils::rotateComplexArray(logIO(), convFunc, convFuncCS_p,
2608 : // rotatedConvFunc,(currentCFPA-actualPA),"LINEAR");
2609 :
2610 0 : ConjCFMap_p = ConjCFMap.getStorage(deleteThem(CONJCFMAP));
2611 0 : CFMap_p = CFMap.getStorage(deleteThem(CFMAP));
2612 :
2613 0 : uvw_p = uvw.getStorage(deleteThem(UVW));
2614 0 : dphase_p = dphase.getStorage(deleteThem(DPHASE));
2615 0 : visdata_p = visdata.getStorage(deleteThem(VISDATA));
2616 0 : gradVisAzData_p = gradVisAzData.getStorage(deleteThem(GRADVISAZ));
2617 0 : gradVisElData_p = gradVisElData.getStorage(deleteThem(GRADVISEL));
2618 0 : flags_p = flags.getStorage(deleteThem(FLAGS));
2619 0 : rowFlags_p = rowFlags.getStorage(deleteThem(ROWFLAGS));
2620 0 : uvScale_p = uvScale.getStorage(deleteThem(UVSCALE));
2621 0 : actualOffset_p = actualOffset.getStorage(deleteThem(ACTUALOFFSET));
2622 0 : dataPtr_p = dataPtr->getStorage(deleteThem(DATAPTR));
2623 0 : vb_freq_p = vb.frequency().getStorage(deleteThem(VBFREQ));
2624 0 : convSupport_p = convSupport.getStorage(deleteThem(CONVSUPPORT));
2625 : // f_convFunc_p = convFunc.getStorage(deleteThem(CONVFUNC));
2626 0 : f_convFunc_p = rotatedConvFunc.getStorage(deleteThem(CONVFUNC));
2627 0 : chanMap_p = chanMap.getStorage(deleteThem(CHANMAP));
2628 0 : polMap_p = polMap.getStorage(deleteThem(POLMAP));
2629 0 : vb_ant1_p = vb.antenna1().getStorage(deleteThem(VBANT1));
2630 0 : vb_ant2_p = vb.antenna2().getStorage(deleteThem(VBANT2));
2631 0 : l_off_p = l_off.getStorage(deleteThem(RAOFF));
2632 0 : m_off_p = m_off.getStorage(deleteThem(DECOFF));
2633 :
2634 0 : Int npa=convSupport.shape()(2),actualConvSize;
2635 0 : Int paIndex_Fortran = paIndex;
2636 0 : actualConvSize = convFunc.shape()(0);
2637 :
2638 0 : IPosition shp=convSupport.shape();
2639 :
2640 0 : dpbwgrad(uvw_p,
2641 : dphase_p,
2642 : // vb.modelVisCube().getStorage(del),
2643 : visdata_p,
2644 0 : &s.asVector()(0),
2645 0 : &s.asVector()(1),
2646 : gradVisAzData_p,
2647 : gradVisElData_p,
2648 : // &gradS(0),
2649 : // &gradS(1),
2650 : &Conj,
2651 : flags_p,
2652 : rowFlags_p,
2653 0 : &s.asVector()(2),
2654 : &rownr,
2655 : uvScale_p,
2656 : actualOffset_p,
2657 : dataPtr_p,
2658 : &aNx,
2659 : &aNy,
2660 : &npol,
2661 : &nchan,
2662 : vb_freq_p,
2663 : &C::c,
2664 : convSupport_p,
2665 : &actualConvSize,
2666 : &convSampling,
2667 : &wConvSize,
2668 : f_convFunc_p,
2669 : chanMap_p,
2670 : polMap_p,
2671 : &polInUse,
2672 : vb_ant1_p,
2673 : vb_ant2_p,
2674 : &Nant_p,
2675 : &scanNo,
2676 : &sigma,
2677 : l_off_p, m_off_p,
2678 : &area,
2679 : &doGrad,
2680 : &doPointing,
2681 : &npa,
2682 : &paIndex_Fortran,
2683 : CFMap_p,
2684 : ConjCFMap_p,
2685 : ¤tCFPA
2686 : ,&actualPA,&cfRefFreq_p
2687 : );
2688 :
2689 0 : ConjCFMap.freeStorage((const Int *&)ConjCFMap_p,deleteThem(CONJCFMAP));
2690 0 : CFMap.freeStorage((const Int *&)CFMap_p,deleteThem(CFMAP));
2691 :
2692 0 : l_off.freeStorage((const Float*&)l_off_p,deleteThem(RAOFF));
2693 0 : m_off.freeStorage((const Float*&)m_off_p,deleteThem(DECOFF));
2694 0 : uvw.freeStorage((const Double*&)uvw_p,deleteThem(UVW));
2695 0 : dphase.freeStorage((const Double*&)dphase_p,deleteThem(DPHASE));
2696 0 : visdata.putStorage(visdata_p,deleteThem(VISDATA));
2697 0 : gradVisAzData.putStorage(gradVisAzData_p,deleteThem(GRADVISAZ));
2698 0 : gradVisElData.putStorage(gradVisElData_p,deleteThem(GRADVISEL));
2699 0 : flags.freeStorage((const Int*&) flags_p,deleteThem(FLAGS));
2700 0 : rowFlags.freeStorage((const Int *&)rowFlags_p,deleteThem(ROWFLAGS));
2701 0 : actualOffset.freeStorage((const Double*&)actualOffset_p,deleteThem(ACTUALOFFSET));
2702 0 : dataPtr->freeStorage((const Complex *&)dataPtr_p,deleteThem(DATAPTR));
2703 0 : uvScale.freeStorage((const Double*&) uvScale_p,deleteThem(UVSCALE));
2704 0 : vb.frequency().freeStorage((const Double*&)vb_freq_p,deleteThem(VBFREQ));
2705 0 : convSupport.freeStorage((const Int*&)convSupport_p,deleteThem(CONVSUPPORT));
2706 0 : convFunc.freeStorage((const Complex *&)f_convFunc_p,deleteThem(CONVFUNC));
2707 0 : chanMap.freeStorage((const Int*&)chanMap_p,deleteThem(CHANMAP));
2708 0 : polMap.freeStorage((const Int*&) polMap_p,deleteThem(POLMAP));
2709 0 : vb.antenna1().freeStorage((const Int*&) vb_ant1_p,deleteThem(VBANT1));
2710 0 : vb.antenna2().freeStorage((const Int*&) vb_ant2_p,deleteThem(VBANT2));
2711 0 : }
2712 : //
2713 : //----------------------------------------------------------------------
2714 : //
2715 0 : void nPBWProjectFT::runFortranPut(Matrix<Double>& uvw,Vector<Double>& dphase,
2716 : const Complex& visdata,
2717 : IPosition& s,
2718 : // Cube<Complex>& gradVisAzData,
2719 : // Cube<Complex>& gradVisElData,
2720 : // IPosition& gradS,
2721 : Int& /*Conj*/,
2722 : Cube<Int>& flags,Vector<Int>& rowFlags,
2723 : const Matrix<Float>& weight,
2724 : Int& rownr,Vector<Double>& actualOffset,
2725 : Array<Complex>& dataPtr,
2726 : Int& aNx, Int& aNy, Int& npol, Int& nchan,
2727 : const VisBuffer& vb,Int& Nant_p, Int& scanNo,
2728 : Double& sigma,
2729 : Array<Float>& l_off,
2730 : Array<Float>& m_off,
2731 : Matrix<Double>& sumWeight,
2732 : Double& area,
2733 : Int& doGrad,
2734 : Int& doPSF,
2735 : Int paIndex)
2736 : {
2737 : enum whichGetStorage {RAOFF,DECOFF,UVW,DPHASE,VISDATA,GRADVISAZ,GRADVISEL,
2738 : FLAGS,ROWFLAGS,UVSCALE,ACTUALOFFSET,DATAPTR,VBFREQ,
2739 : CONVSUPPORT,CONVFUNC,CHANMAP,POLMAP,VBANT1,VBANT2,WEIGHT,
2740 : SUMWEIGHT,CONJCFMAP,CFMAP};
2741 0 : Vector<Bool> deleteThem(23);
2742 :
2743 : Double *uvw_p, *dphase_p, *actualOffset_p, *vb_freq_p, *uvScale_p;
2744 : Complex *dataPtr_p, *f_convFunc_p;
2745 : // Complex *gradVisAzData_p, *gradVisElData_p;
2746 : Int *flags_p, *rowFlags_p, *chanMap_p, *polMap_p, *convSupport_p, *vb_ant1_p, *vb_ant2_p,
2747 : *ConjCFMap_p, *CFMap_p;
2748 : Float *l_off_p, *m_off_p;
2749 : Float *weight_p;Double *sumwt_p;
2750 : Double actualPA;
2751 0 : const Complex *visdata_p=&visdata;
2752 :
2753 0 : Vector<Int> ConjCFMap, CFMap;
2754 0 : actualPA = getVBPA(vb);
2755 0 : ConjCFMap = polMap;
2756 :
2757 0 : Array<Complex> rotatedConvFunc;
2758 : // SynthesisUtils::rotateComplexArray(logIO(), convFunc, convFuncCS_p,
2759 : // rotatedConvFunc,(currentCFPA-actualPA),"LINEAR");
2760 0 : SynthesisUtils::rotateComplexArray(logIO(), convFunc, convFuncCS_p,
2761 : rotatedConvFunc,0.0,"LINEAR");
2762 :
2763 : /*
2764 : CFMap = polMap; ConjCFMap = polMap;
2765 : CFMap = makeConjPolMap(vb);
2766 : */
2767 0 : makeCFPolMap(vb,CFMap);
2768 0 : makeConjPolMap(vb,CFMap,ConjCFMap);
2769 :
2770 0 : ConjCFMap_p = ConjCFMap.getStorage(deleteThem(CONJCFMAP));
2771 0 : CFMap_p = CFMap.getStorage(deleteThem(CFMAP));
2772 :
2773 0 : uvw_p = uvw.getStorage(deleteThem(UVW));
2774 0 : dphase_p = dphase.getStorage(deleteThem(DPHASE));
2775 : // visdata_p = visdata.getStorage(deleteThem(VISDATA));
2776 : // gradVisAzData_p = gradVisAzData.getStorage(deleteThem(GRADVISAZ));
2777 : // gradVisElData_p = gradVisElData.getStorage(deleteThem(GRADVISEL));
2778 0 : flags_p = flags.getStorage(deleteThem(FLAGS));
2779 0 : rowFlags_p = rowFlags.getStorage(deleteThem(ROWFLAGS));
2780 0 : uvScale_p = uvScale.getStorage(deleteThem(UVSCALE));
2781 0 : actualOffset_p = actualOffset.getStorage(deleteThem(ACTUALOFFSET));
2782 0 : dataPtr_p = dataPtr.getStorage(deleteThem(DATAPTR));
2783 0 : vb_freq_p = (Double *)(vb.frequency().getStorage(deleteThem(VBFREQ)));
2784 0 : convSupport_p = convSupport.getStorage(deleteThem(CONVSUPPORT));
2785 : // f_convFunc_p = convFunc.getStorage(deleteThem(CONVFUNC));
2786 0 : f_convFunc_p = rotatedConvFunc.getStorage(deleteThem(CONVFUNC));
2787 0 : chanMap_p = chanMap.getStorage(deleteThem(CHANMAP));
2788 0 : polMap_p = polMap.getStorage(deleteThem(POLMAP));
2789 0 : vb_ant1_p = (Int *)(vb.antenna1().getStorage(deleteThem(VBANT1)));
2790 0 : vb_ant2_p = (Int *)(vb.antenna2().getStorage(deleteThem(VBANT2)));
2791 0 : l_off_p = l_off.getStorage(deleteThem(RAOFF));
2792 0 : m_off_p = m_off.getStorage(deleteThem(DECOFF));
2793 0 : weight_p = (Float *)(weight.getStorage(deleteThem(WEIGHT)));
2794 0 : sumwt_p = sumWeight.getStorage(deleteThem(SUMWEIGHT));
2795 :
2796 :
2797 0 : Int npa=convSupport.shape()(2),actualConvSize;
2798 0 : Int paIndex_Fortran = paIndex;
2799 0 : actualConvSize = convFunc.shape()(0);
2800 :
2801 0 : IPosition shp=convSupport.shape();
2802 :
2803 0 : gpbwproj(uvw_p,
2804 : dphase_p,
2805 : // vb.modelVisCube().getStorage(del),
2806 : visdata_p,
2807 0 : &s.asVector()(0),
2808 0 : &s.asVector()(1),
2809 : // gradVisAzData_p,
2810 : // gradVisElData_p,
2811 : // &gradS(0),
2812 : // &gradS(1),
2813 : // &Conj,
2814 : &doPSF,
2815 : flags_p,
2816 : rowFlags_p,
2817 : weight_p,
2818 0 : &s.asVector()(2),
2819 : &rownr,
2820 : uvScale_p,
2821 : actualOffset_p,
2822 : dataPtr_p,
2823 : &aNx,
2824 : &aNy,
2825 : &npol,
2826 : &nchan,
2827 : vb_freq_p,
2828 : &C::c,
2829 : convSupport_p,
2830 : &actualConvSize,
2831 : &convSampling,
2832 : &wConvSize,
2833 : f_convFunc_p,
2834 : chanMap_p,
2835 : polMap_p,
2836 : &polInUse,
2837 : sumwt_p,
2838 : vb_ant1_p,
2839 : vb_ant2_p,
2840 : &Nant_p,
2841 : &scanNo,
2842 : &sigma,
2843 : l_off_p, m_off_p,
2844 : &area,
2845 : &doGrad,
2846 : &doPointing,
2847 : &npa,
2848 : &paIndex_Fortran,
2849 : CFMap_p,
2850 : ConjCFMap_p,
2851 : ¤tCFPA
2852 : ,&actualPA,&cfRefFreq_p
2853 : );
2854 :
2855 0 : ConjCFMap.freeStorage((const Int *&)ConjCFMap_p,deleteThem(CONJCFMAP));
2856 0 : CFMap.freeStorage((const Int *&)CFMap_p,deleteThem(CFMAP));
2857 :
2858 0 : l_off.freeStorage((const Float*&)l_off_p,deleteThem(RAOFF));
2859 0 : m_off.freeStorage((const Float*&)m_off_p,deleteThem(DECOFF));
2860 0 : uvw.freeStorage((const Double*&)uvw_p,deleteThem(UVW));
2861 0 : dphase.freeStorage((const Double*&)dphase_p,deleteThem(DPHASE));
2862 : // visdata.putStorage(visdata_p,deleteThem(VISDATA));
2863 : // gradVisAzData.putStorage(gradVisAzData_p,deleteThem(GRADVISAZ));
2864 : // gradVisElData.putStorage(gradVisElData_p,deleteThem(GRADVISEL));
2865 0 : flags.freeStorage((const Int*&) flags_p,deleteThem(FLAGS));
2866 0 : rowFlags.freeStorage((const Int *&)rowFlags_p,deleteThem(ROWFLAGS));
2867 0 : actualOffset.freeStorage((const Double*&)actualOffset_p,deleteThem(ACTUALOFFSET));
2868 0 : dataPtr.freeStorage((const Complex *&)dataPtr_p,deleteThem(DATAPTR));
2869 0 : uvScale.freeStorage((const Double*&) uvScale_p,deleteThem(UVSCALE));
2870 0 : vb.frequency().freeStorage((const Double*&)vb_freq_p,deleteThem(VBFREQ));
2871 0 : convSupport.freeStorage((const Int*&)convSupport_p,deleteThem(CONVSUPPORT));
2872 0 : convFunc.freeStorage((const Complex *&)f_convFunc_p,deleteThem(CONVFUNC));
2873 0 : chanMap.freeStorage((const Int*&)chanMap_p,deleteThem(CHANMAP));
2874 0 : polMap.freeStorage((const Int*&) polMap_p,deleteThem(POLMAP));
2875 0 : vb.antenna1().freeStorage((const Int*&) vb_ant1_p,deleteThem(VBANT1));
2876 0 : vb.antenna2().freeStorage((const Int*&) vb_ant2_p,deleteThem(VBANT2));
2877 0 : weight.freeStorage((const Float*&)weight_p,deleteThem(WEIGHT));
2878 0 : sumWeight.putStorage(sumwt_p,deleteThem(SUMWEIGHT));
2879 0 : }
2880 : //
2881 : //---------------------------------------------------------------
2882 : //
2883 0 : void nPBWProjectFT::put(const VisBuffer& vb, Int row, Bool dopsf,
2884 : FTMachine::Type type)
2885 : {
2886 : // Take care of translation of Bools to Integer
2887 0 : makingPSF=dopsf;
2888 :
2889 0 : findConvFunction(*image, vb);
2890 :
2891 :
2892 : const Matrix<Float> *imagingweight;
2893 0 : imagingweight=&(vb.imagingWeight());
2894 :
2895 : const Cube<Complex> *data;
2896 0 : if(type==FTMachine::MODEL)
2897 0 : data=&(vb.modelVisCube());
2898 0 : else if(type==FTMachine::CORRECTED)
2899 0 : data=&(vb.correctedVisCube());
2900 : else
2901 0 : data=&(vb.visCube());
2902 :
2903 : Bool isCopy;
2904 0 : const casacore::Complex *datStorage=data->getStorage(isCopy);
2905 0 : Int NAnt = 0;
2906 :
2907 0 : if (doPointing) NAnt = findPointingOffsets(vb,l_offsets,m_offsets,true);
2908 :
2909 : //
2910 : // If row is -1 then we pass through all rows
2911 : //
2912 : Int startRow, endRow, nRow;
2913 0 : if (row==-1)
2914 : {
2915 0 : nRow=vb.nRow();
2916 0 : startRow=0;
2917 0 : endRow=nRow-1;
2918 : }
2919 : else
2920 : {
2921 0 : nRow=1;
2922 0 : startRow=row;
2923 0 : endRow=row;
2924 : }
2925 : //
2926 : // Get the uvws in a form that Fortran can use and do that
2927 : // necessary phase rotation. On a Pentium Pro 200 MHz when null,
2928 : // this step takes about 50us per uvw point. This is just barely
2929 : // noticeable for Stokes I continuum and irrelevant for other
2930 : // cases.
2931 : //
2932 0 : Matrix<Double> uvw(3, vb.uvw().nelements());
2933 0 : uvw=0.0;
2934 0 : Vector<Double> dphase(vb.uvw().nelements());
2935 0 : dphase=0.0;
2936 : //NEGATING to correct for an image inversion problem
2937 0 : for (Int i=startRow;i<=endRow;i++)
2938 : {
2939 0 : for (Int idim=0;idim<2;idim++) uvw(idim,i)=-vb.uvw()(i)(idim);
2940 0 : uvw(2,i)=vb.uvw()(i)(2);
2941 : }
2942 :
2943 0 : rotateUVW(uvw, dphase, vb);
2944 0 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
2945 :
2946 : // This is the convention for dphase
2947 0 : dphase*=-1.0;
2948 :
2949 0 : Cube<Int> flags(vb.flagCube().shape());
2950 0 : flags=0;
2951 0 : flags(vb.flagCube())=true;
2952 :
2953 0 : Vector<Int> rowFlags(vb.nRow());
2954 0 : rowFlags=0;
2955 0 : rowFlags(vb.flagRow())=true;
2956 0 : if(!usezero_p)
2957 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
2958 0 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
2959 : //Check if ms has changed then cache new spw and chan selection
2960 0 : if(vb.newMS())
2961 0 : matchAllSpwChans(vb);
2962 :
2963 : //Here we redo the match or use previous match
2964 :
2965 : //Channel matching for the actual spectral window of buffer
2966 0 : if(doConversion_p[vb.spectralWindow()])
2967 0 : matchChannel(vb.spectralWindow(), vb);
2968 : else
2969 : {
2970 0 : chanMap.resize();
2971 0 : chanMap=multiChanMap_p[vb.spectralWindow()];
2972 : }
2973 :
2974 0 : if(isTiled)
2975 : {// Tiled Version
2976 0 : Double invLambdaC=vb.frequency()(0)/C::c;
2977 0 : Vector<Double> uvLambda(2);
2978 0 : Vector<Int> centerLoc2D(2);
2979 0 : centerLoc2D=0;
2980 : //
2981 : // Loop over all rows
2982 : //
2983 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
2984 : {
2985 : // Calculate uvw for this row at the center frequency
2986 0 : uvLambda(0)=uvw(0,rownr)*invLambdaC;
2987 0 : uvLambda(1)=uvw(1,rownr)*invLambdaC;
2988 0 : centerLoc2D=gridder->location(centerLoc2D, uvLambda);
2989 : //
2990 : // Is this point on the grid?
2991 : //
2992 0 : if(gridder->onGrid(centerLoc2D))
2993 : {
2994 : // Get the tile
2995 0 : Array<Complex>* dataPtr=getDataPointer(centerLoc2D, false);
2996 0 : Int aNx=dataPtr->shape()(0);
2997 0 : Int aNy=dataPtr->shape()(1);
2998 : //
2999 : // Now use FORTRAN to do the gridding. Remember to
3000 : // ensure that the shape and offsets of the tile are
3001 : // accounted for.
3002 : //
3003 0 : Vector<Double> actualOffset(3);
3004 0 : for (Int i=0;i<2;i++) actualOffset(i)=uvOffset(i)-Double(offsetLoc(i));
3005 :
3006 0 : actualOffset(2)=uvOffset(2);
3007 0 : IPosition s(flags.shape());
3008 : //
3009 : // Now pass all the information down to a FORTRAN routine to
3010 : // do the work
3011 : //
3012 0 : Int Conj=0,doPSF;
3013 0 : Int ScanNo=0,doGrad=0;
3014 0 : Double area=1.0;
3015 :
3016 0 : Int tmpPAI=PAIndex+1;
3017 0 : if (dopsf) doPSF=1; else doPSF=0;
3018 0 : runFortranPut(uvw,dphase,*datStorage,s,Conj,flags,rowFlags,
3019 : *imagingweight,rownr,actualOffset,
3020 0 : *dataPtr,aNx,aNy,npol,nchan,vb,NAnt,ScanNo,sigma,
3021 0 : l_offsets,m_offsets,sumWeight,area,doGrad,doPSF,tmpPAI);
3022 0 : }
3023 : }
3024 0 : }
3025 : else
3026 : {//Non-tiled version
3027 0 : IPosition s(flags.shape());
3028 :
3029 0 : Int Conj=0,doPSF=0;
3030 0 : Int ScanNo=0,doGrad=0;Double area=1.0;
3031 :
3032 0 : if (dopsf) doPSF=1;
3033 :
3034 0 : Int tmpPAI=PAIndex+1;
3035 0 : runFortranPut(uvw,dphase,*datStorage,s,Conj,flags,rowFlags,
3036 : *imagingweight,
3037 0 : row,uvOffset,
3038 0 : griddedData,nx,ny,npol,nchan,vb,NAnt,ScanNo,sigma,
3039 0 : l_offsets,m_offsets,sumWeight,area,doGrad,doPSF,tmpPAI);
3040 0 : }
3041 :
3042 0 : data->freeStorage(datStorage, isCopy);
3043 0 : }
3044 : //
3045 : //----------------------------------------------------------------------
3046 : //
3047 0 : void nPBWProjectFT::initVisBuffer(VisBuffer& vb, Type whichVBColumn)
3048 : {
3049 0 : if (whichVBColumn == FTMachine::MODEL) vb.modelVisCube()=Complex(0.0,0.0);
3050 0 : else if (whichVBColumn == FTMachine::OBSERVED) vb.visCube()=Complex(0.0,0.0);
3051 0 : }
3052 : //
3053 : //----------------------------------------------------------------------
3054 : //
3055 0 : void nPBWProjectFT::initVisBuffer(VisBuffer& vb, Type whichVBColumn, Int row)
3056 : {
3057 0 : if (whichVBColumn == FTMachine::MODEL)
3058 0 : vb.modelVisCube().xyPlane(row)=Complex(0.0,0.0);
3059 0 : else if (whichVBColumn == FTMachine::OBSERVED)
3060 0 : vb.visCube().xyPlane(row)=Complex(0.0,0.0);
3061 0 : }
3062 : //
3063 : //---------------------------------------------------------------
3064 : //
3065 : // Predict the coherences as well as their derivatives w.r.t. the
3066 : // pointing offsets.
3067 : //
3068 0 : void nPBWProjectFT::nget(VisBuffer& vb,
3069 : // These offsets should be appropriate for the VB
3070 : Array<Float>& l_off, Array<Float>& m_off,
3071 : Cube<Complex>& Mout,
3072 : Cube<Complex>& dMout1,
3073 : Cube<Complex>& dMout2,
3074 : Int Conj, Int doGrad)
3075 : {
3076 : Int startRow, endRow, nRow;
3077 0 : nRow=vb.nRow();
3078 0 : startRow=0;
3079 0 : endRow=nRow-1;
3080 :
3081 0 : Mout = dMout1 = dMout2 = Complex(0,0);
3082 :
3083 0 : findConvFunction(*image, vb);
3084 0 : Int NAnt=0;
3085 0 : if (bandID_p == -1) bandID_p=getVisParams(vb);
3086 0 : if (doPointing)
3087 0 : NAnt = findPointingOffsets(vb,l_offsets,m_offsets,false);
3088 :
3089 0 : l_offsets=l_off;
3090 0 : m_offsets=m_off;
3091 0 : Matrix<Double> uvw(3, vb.uvw().nelements());
3092 0 : uvw=0.0;
3093 0 : Vector<Double> dphase(vb.uvw().nelements());
3094 0 : dphase=0.0;
3095 : //NEGATING to correct for an image inversion problem
3096 0 : for (Int i=startRow;i<=endRow;i++)
3097 : {
3098 0 : for (Int idim=0;idim<2;idim++) uvw(idim,i)=-vb.uvw()(i)(idim);
3099 0 : uvw(2,i)=vb.uvw()(i)(2);
3100 : }
3101 :
3102 0 : rotateUVW(uvw, dphase, vb);
3103 0 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
3104 :
3105 : // This is the convention for dphase
3106 0 : dphase*=-1.0;
3107 :
3108 0 : Cube<Int> flags(vb.flagCube().shape());
3109 0 : flags=0;
3110 0 : flags(vb.flagCube())=true;
3111 :
3112 : //Check if ms has changed then cache new spw and chan selection
3113 0 : if(vb.newMS())
3114 0 : matchAllSpwChans(vb);
3115 :
3116 : //Here we redo the match or use previous match
3117 : //
3118 : //Channel matching for the actual spectral window of buffer
3119 : //
3120 0 : if(doConversion_p[vb.spectralWindow()])
3121 0 : matchChannel(vb.spectralWindow(), vb);
3122 : else
3123 : {
3124 0 : chanMap.resize();
3125 0 : chanMap=multiChanMap_p[vb.spectralWindow()];
3126 : }
3127 :
3128 0 : Vector<Int> rowFlags(vb.nRow());
3129 0 : rowFlags=0;
3130 0 : rowFlags(vb.flagRow())=true;
3131 0 : if(!usezero_p)
3132 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
3133 0 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
3134 :
3135 0 : IPosition s,gradS;
3136 0 : Cube<Complex> visdata,gradVisAzData,gradVisElData;
3137 : //
3138 : // visdata now references the Mout data structure rather than to the internal VB storeage.
3139 : //
3140 0 : visdata.reference(Mout);
3141 :
3142 0 : if (doGrad)
3143 : {
3144 : // The following should reference some slice of dMout?
3145 0 : gradVisAzData.reference(dMout1);
3146 0 : gradVisElData.reference(dMout2);
3147 : }
3148 : //
3149 : // Begin the actual de-gridding.
3150 : //
3151 0 : if(isTiled)
3152 : {
3153 0 : logIO() << "nPBWProjectFT::nget(): The sky model is tiled" << LogIO::NORMAL << LogIO::POST;
3154 0 : Double invLambdaC=vb.frequency()(0)/C::c;
3155 0 : Vector<Double> uvLambda(2);
3156 0 : Vector<Int> centerLoc2D(2);
3157 0 : centerLoc2D=0;
3158 :
3159 : // Loop over all rows
3160 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
3161 : {
3162 :
3163 : // Calculate uvw for this row at the center frequency
3164 0 : uvLambda(0)=uvw(0, rownr)*invLambdaC;
3165 0 : uvLambda(1)=uvw(1, rownr)*invLambdaC;
3166 0 : centerLoc2D=gridder->location(centerLoc2D, uvLambda);
3167 :
3168 : // Is this point on the grid?
3169 0 : if(gridder->onGrid(centerLoc2D))
3170 : {
3171 :
3172 : // Get the tile
3173 0 : Array<Complex>* dataPtr=getDataPointer(centerLoc2D, true);
3174 0 : gridder->setOffset(IPosition(2, offsetLoc(0), offsetLoc(1)));
3175 0 : Int aNx=dataPtr->shape()(0);
3176 0 : Int aNy=dataPtr->shape()(1);
3177 :
3178 : // Now use FORTRAN to do the gridding. Remember to
3179 : // ensure that the shape and offsets of the tile are
3180 : // accounted for.
3181 :
3182 0 : Vector<Double> actualOffset(3);
3183 0 : for (Int i=0;i<2;i++)
3184 0 : actualOffset(i)=uvOffset(i)-Double(offsetLoc(i));
3185 :
3186 0 : actualOffset(2)=uvOffset(2);
3187 0 : IPosition s(vb.modelVisCube().shape());
3188 :
3189 0 : Int ScanNo=0, tmpPAI;
3190 0 : Double area=1.0;
3191 0 : tmpPAI = PAIndex+1;
3192 0 : runFortranGetGrad(uvw,dphase,visdata,s,
3193 : gradVisAzData,gradVisElData,
3194 : Conj,flags,rowFlags,rownr,
3195 0 : actualOffset,dataPtr,aNx,aNy,npol,nchan,vb,NAnt,ScanNo,sigma,
3196 0 : l_offsets,m_offsets,area,doGrad,tmpPAI);
3197 0 : }
3198 : }
3199 0 : }
3200 : else
3201 : {
3202 0 : IPosition s(vb.modelVisCube().shape());
3203 0 : Int ScanNo=0, tmpPAI, trow=-1;
3204 0 : Double area=1.0;
3205 0 : tmpPAI = PAIndex+1;
3206 0 : runFortranGetGrad(uvw,dphase,visdata/*vb.modelVisCube()*/,s,
3207 : gradVisAzData, gradVisElData,
3208 : Conj,flags,rowFlags,trow,
3209 0 : uvOffset,&griddedData,nx,ny,npol,nchan,vb,NAnt,ScanNo,sigma,
3210 0 : l_offsets,m_offsets,area,doGrad,tmpPAI);
3211 0 : }
3212 :
3213 0 : }
3214 0 : void nPBWProjectFT::get(VisBuffer& vb,
3215 : VisBuffer& gradVBAz,
3216 : VisBuffer& gradVBEl,
3217 : Cube<Float>& pointingOffsets,
3218 : Int row, // default row=-1
3219 : Type whichVBColumn, // default whichVBColumn = FTMachine::MODEL
3220 : Type whichGradVBColumn,// default whichGradVBColumn = FTMachine::MODEL
3221 : Int Conj, Int doGrad) // default Conj=0, doGrad=1
3222 : {
3223 : // If row is -1 then we pass through all rows
3224 : Int startRow, endRow, nRow;
3225 0 : if (row==-1)
3226 : {
3227 0 : nRow=vb.nRow();
3228 0 : startRow=0;
3229 0 : endRow=nRow-1;
3230 0 : initVisBuffer(vb,whichVBColumn);
3231 0 : if (doGrad)
3232 : {
3233 0 : initVisBuffer(gradVBAz, whichGradVBColumn);
3234 0 : initVisBuffer(gradVBEl, whichGradVBColumn);
3235 : }
3236 : }
3237 : else
3238 : {
3239 0 : nRow=1;
3240 0 : startRow=row;
3241 0 : endRow=row;
3242 0 : initVisBuffer(vb,whichVBColumn,row);
3243 0 : if (doGrad)
3244 : {
3245 0 : initVisBuffer(gradVBAz, whichGradVBColumn,row);
3246 0 : initVisBuffer(gradVBEl, whichGradVBColumn,row);
3247 : }
3248 : }
3249 :
3250 0 : findConvFunction(*image, vb);
3251 :
3252 0 : if (bandID_p == -1) bandID_p=getVisParams(vb);
3253 0 : Int NAnt=0;
3254 0 : if (doPointing)
3255 0 : NAnt = findPointingOffsets(vb,pointingOffsets,l_offsets,m_offsets,false);
3256 :
3257 0 : Matrix<Double> uvw(3, vb.uvw().nelements());
3258 0 : uvw=0.0;
3259 0 : Vector<Double> dphase(vb.uvw().nelements());
3260 0 : dphase=0.0;
3261 : //NEGATING to correct for an image inversion problem
3262 0 : for (Int i=startRow;i<=endRow;i++)
3263 : {
3264 0 : for (Int idim=0;idim<2;idim++) uvw(idim,i)=-vb.uvw()(i)(idim);
3265 0 : uvw(2,i)=vb.uvw()(i)(2);
3266 : }
3267 :
3268 0 : rotateUVW(uvw, dphase, vb);
3269 0 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
3270 :
3271 : // This is the convention for dphase
3272 0 : dphase*=-1.0;
3273 :
3274 :
3275 0 : Cube<Int> flags(vb.flagCube().shape());
3276 0 : flags=0;
3277 0 : flags(vb.flagCube())=true;
3278 : //
3279 : //Check if ms has changed then cache new spw and chan selection
3280 : //
3281 0 : if(vb.newMS()) matchAllSpwChans(vb);
3282 :
3283 : //Here we redo the match or use previous match
3284 : //
3285 : //Channel matching for the actual spectral window of buffer
3286 : //
3287 0 : if(doConversion_p[vb.spectralWindow()])
3288 0 : matchChannel(vb.spectralWindow(), vb);
3289 : else
3290 : {
3291 0 : chanMap.resize();
3292 0 : chanMap=multiChanMap_p[vb.spectralWindow()];
3293 : }
3294 :
3295 0 : Vector<Int> rowFlags(vb.nRow());
3296 0 : rowFlags=0;
3297 0 : rowFlags(vb.flagRow())=true;
3298 0 : if(!usezero_p)
3299 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
3300 0 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
3301 :
3302 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
3303 0 : if (vb.antenna1()(rownr) != vb.antenna2()(rownr))
3304 0 : rowFlags(rownr) = (vb.flagRow()(rownr)==true);
3305 :
3306 0 : IPosition s,gradS;
3307 0 : Cube<Complex> visdata,gradVisAzData,gradVisElData;
3308 0 : if (whichVBColumn == FTMachine::MODEL)
3309 : {
3310 0 : s = vb.modelVisCube().shape();
3311 0 : visdata.reference(vb.modelVisCube());
3312 : }
3313 0 : else if (whichVBColumn == FTMachine::OBSERVED)
3314 : {
3315 0 : s = vb.visCube().shape();
3316 0 : visdata.reference(vb.visCube());
3317 : }
3318 :
3319 0 : if (doGrad)
3320 : {
3321 0 : if (whichGradVBColumn == FTMachine::MODEL)
3322 : {
3323 : // gradS = gradVBAz.modelVisCube().shape();
3324 0 : gradVisAzData.reference(gradVBAz.modelVisCube());
3325 0 : gradVisElData.reference(gradVBEl.modelVisCube());
3326 : }
3327 0 : else if (whichGradVBColumn == FTMachine::OBSERVED)
3328 : {
3329 : // gradS = gradVBAz.visCube().shape();
3330 0 : gradVisAzData.reference(gradVBAz.visCube());
3331 0 : gradVisElData.reference(gradVBEl.visCube());
3332 : }
3333 : }
3334 :
3335 0 : if(isTiled)
3336 : {
3337 0 : Double invLambdaC=vb.frequency()(0)/C::c;
3338 0 : Vector<Double> uvLambda(2);
3339 0 : Vector<Int> centerLoc2D(2);
3340 0 : centerLoc2D=0;
3341 :
3342 : // Loop over all rows
3343 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
3344 : {
3345 :
3346 : // Calculate uvw for this row at the center frequency
3347 0 : uvLambda(0)=uvw(0, rownr)*invLambdaC;
3348 0 : uvLambda(1)=uvw(1, rownr)*invLambdaC;
3349 0 : centerLoc2D=gridder->location(centerLoc2D, uvLambda);
3350 :
3351 : // Is this point on the grid?
3352 0 : if(gridder->onGrid(centerLoc2D))
3353 : {
3354 :
3355 : // Get the tile
3356 0 : Array<Complex>* dataPtr=getDataPointer(centerLoc2D, true);
3357 0 : gridder->setOffset(IPosition(2, offsetLoc(0), offsetLoc(1)));
3358 0 : Int aNx=dataPtr->shape()(0);
3359 0 : Int aNy=dataPtr->shape()(1);
3360 :
3361 : // Now use FORTRAN to do the gridding. Remember to
3362 : // ensure that the shape and offsets of the tile are
3363 : // accounted for.
3364 :
3365 0 : Vector<Double> actualOffset(3);
3366 0 : for (Int i=0;i<2;i++)
3367 0 : actualOffset(i)=uvOffset(i)-Double(offsetLoc(i));
3368 :
3369 0 : actualOffset(2)=uvOffset(2);
3370 0 : IPosition s(vb.modelVisCube().shape());
3371 :
3372 0 : Int ScanNo=0, tmpPAI;
3373 0 : Double area=1.0;
3374 0 : tmpPAI = PAIndex+1;
3375 0 : runFortranGetGrad(uvw,dphase,visdata,s,
3376 : gradVisAzData,gradVisElData,
3377 : Conj,flags,rowFlags,rownr,
3378 0 : actualOffset,dataPtr,aNx,aNy,npol,nchan,vb,NAnt,ScanNo,sigma,
3379 0 : l_offsets,m_offsets,area,doGrad,tmpPAI);
3380 0 : }
3381 : }
3382 0 : }
3383 : else
3384 : {
3385 :
3386 0 : IPosition s(vb.modelVisCube().shape());
3387 0 : Int ScanNo=0, tmpPAI;
3388 0 : Double area=1.0;
3389 :
3390 0 : tmpPAI = PAIndex+1;
3391 :
3392 0 : runFortranGetGrad(uvw,dphase,visdata/*vb.modelVisCube()*/,s,
3393 : gradVisAzData, gradVisElData,
3394 : Conj,flags,rowFlags,row,
3395 0 : uvOffset,&griddedData,nx,ny,npol,nchan,vb,NAnt,ScanNo,sigma,
3396 0 : l_offsets,m_offsets,area,doGrad,tmpPAI);
3397 : // runFortranGet(uvw,dphase,vb.modelVisCube(),s,Conj,flags,rowFlags,row,
3398 : // uvOffset,&griddedData,nx,ny,npol,nchan,vb,NAnt,ScanNo,sigma,
3399 : // l_offsets,m_offsets,area,doGrad,tmpPAI);
3400 0 : }
3401 0 : }
3402 : //
3403 : //---------------------------------------------------------------
3404 : //
3405 0 : void nPBWProjectFT::get(VisBuffer& vb, Int row)
3406 : {
3407 : // If row is -1 then we pass through all rows
3408 : Int startRow, endRow, nRow;
3409 0 : if (row==-1)
3410 : {
3411 0 : nRow=vb.nRow();
3412 0 : startRow=0;
3413 0 : endRow=nRow-1;
3414 0 : vb.modelVisCube()=Complex(0.0,0.0);
3415 : }
3416 : else
3417 : {
3418 0 : nRow=1;
3419 0 : startRow=row;
3420 0 : endRow=row;
3421 0 : vb.modelVisCube().xyPlane(row)=Complex(0.0,0.0);
3422 : }
3423 :
3424 0 : findConvFunction(*image, vb);
3425 :
3426 0 : if (bandID_p == -1) bandID_p=getVisParams(vb);
3427 0 : Int NAnt=0;
3428 0 : if (doPointing) NAnt = findPointingOffsets(vb,l_offsets,m_offsets,true);
3429 :
3430 : // Get the uvws in a form that Fortran can use
3431 0 : Matrix<Double> uvw(3, vb.uvw().nelements());
3432 0 : uvw=0.0;
3433 0 : Vector<Double> dphase(vb.uvw().nelements());
3434 0 : dphase=0.0;
3435 : //NEGATING to correct for an image inversion problem
3436 0 : for (Int i=startRow;i<=endRow;i++)
3437 : {
3438 0 : for (Int idim=0;idim<2;idim++) uvw(idim,i)=-vb.uvw()(i)(idim);
3439 0 : uvw(2,i)=vb.uvw()(i)(2);
3440 : }
3441 :
3442 0 : rotateUVW(uvw, dphase, vb);
3443 0 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
3444 :
3445 : // This is the convention for dphase
3446 0 : dphase*=-1.0;
3447 :
3448 :
3449 0 : Cube<Int> flags(vb.flagCube().shape());
3450 0 : flags=0;
3451 0 : flags(vb.flagCube())=true;
3452 :
3453 : //Check if ms has changed then cache new spw and chan selection
3454 0 : if(vb.newMS())
3455 0 : matchAllSpwChans(vb);
3456 :
3457 : //Here we redo the match or use previous match
3458 : //
3459 : //Channel matching for the actual spectral window of buffer
3460 : //
3461 0 : if(doConversion_p[vb.spectralWindow()])
3462 0 : matchChannel(vb.spectralWindow(), vb);
3463 : else
3464 : {
3465 0 : chanMap.resize();
3466 0 : chanMap=multiChanMap_p[vb.spectralWindow()];
3467 : }
3468 :
3469 0 : Vector<Int> rowFlags(vb.nRow());
3470 0 : rowFlags=0;
3471 0 : rowFlags(vb.flagRow())=true;
3472 :
3473 0 : if(!usezero_p)
3474 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
3475 0 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
3476 :
3477 0 : if(isTiled)
3478 : {
3479 :
3480 0 : Double invLambdaC=vb.frequency()(0)/C::c;
3481 0 : Vector<Double> uvLambda(2);
3482 0 : Vector<Int> centerLoc2D(2);
3483 0 : centerLoc2D=0;
3484 :
3485 : // Loop over all rows
3486 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
3487 : {
3488 :
3489 : // Calculate uvw for this row at the center frequency
3490 0 : uvLambda(0)=uvw(0, rownr)*invLambdaC;
3491 0 : uvLambda(1)=uvw(1, rownr)*invLambdaC;
3492 0 : centerLoc2D=gridder->location(centerLoc2D, uvLambda);
3493 :
3494 : // Is this point on the grid?
3495 0 : if(gridder->onGrid(centerLoc2D))
3496 : {
3497 :
3498 : // Get the tile
3499 0 : Array<Complex>* dataPtr=getDataPointer(centerLoc2D, true);
3500 0 : gridder->setOffset(IPosition(2, offsetLoc(0), offsetLoc(1)));
3501 0 : Int aNx=dataPtr->shape()(0);
3502 0 : Int aNy=dataPtr->shape()(1);
3503 :
3504 : // Now use FORTRAN to do the gridding. Remember to
3505 : // ensure that the shape and offsets of the tile are
3506 : // accounted for.
3507 :
3508 0 : Vector<Double> actualOffset(3);
3509 0 : for (Int i=0;i<2;i++)
3510 0 : actualOffset(i)=uvOffset(i)-Double(offsetLoc(i));
3511 :
3512 0 : actualOffset(2)=uvOffset(2);
3513 0 : IPosition s(vb.modelVisCube().shape());
3514 :
3515 0 : Int Conj=0,doGrad=0,ScanNo=0;
3516 0 : Double area=1.0;
3517 :
3518 0 : runFortranGet(uvw,dphase,vb.modelVisCube(),s,Conj,flags,rowFlags,rownr,
3519 0 : actualOffset,dataPtr,aNx,aNy,npol,nchan,vb,NAnt,ScanNo,sigma,
3520 0 : l_offsets,m_offsets,area,doGrad,PAIndex+1);
3521 0 : }
3522 : }
3523 0 : }
3524 : else
3525 : {
3526 :
3527 0 : IPosition s(vb.modelVisCube().shape());
3528 0 : Int Conj=0,doGrad=0,ScanNo=0;
3529 0 : Double area=1.0;
3530 :
3531 0 : runFortranGet(uvw,dphase,vb.modelVisCube(),s,Conj,flags,rowFlags,row,
3532 0 : uvOffset,&griddedData,nx,ny,npol,nchan,vb,NAnt,ScanNo,sigma,
3533 0 : l_offsets,m_offsets,area,doGrad,PAIndex+1);
3534 : /*
3535 : static int junk=0;
3536 : if (junk==4)
3537 : {
3538 : cout << "Time = " << vb.time()/1e9 << endl;
3539 : for(Int i=0;i<vb.modelVisCube().shape()(2);i++)
3540 : cout << "PBWP: Residual: " << i
3541 : << " " << vb.modelVisCube()(0,0,i)
3542 : << " " << vb.modelVisCube()(3,0,i)
3543 : << " " << vb.visCube()(0,0,i)
3544 : << " " << vb.visCube()(3,0,i)
3545 : << " " << vb.flag()(0,i)
3546 : << " " << vb.antenna1()(i)<< "-" << vb.antenna2()(i)
3547 : << " " << vb.flagRow()(i)
3548 : << " " << vb.flagCube()(0,0,i)
3549 : << " " << vb.flagCube()(3,0,i)
3550 : << endl;
3551 : }
3552 : junk++;
3553 : */
3554 0 : }
3555 0 : }
3556 : //
3557 : //---------------------------------------------------------------
3558 : //
3559 0 : void nPBWProjectFT::get(VisBuffer& vb, Cube<Complex>& modelVis,
3560 : Array<Complex>& griddedVis, Vector<Double>& /*scale*/,
3561 : Int row)
3562 : {
3563 0 : Int nX=griddedVis.shape()(0);
3564 0 : Int nY=griddedVis.shape()(1);
3565 0 : Vector<Double> offset(2);
3566 0 : offset(0)=Double(nX)/2.0;
3567 0 : offset(1)=Double(nY)/2.0;
3568 : // If row is -1 then we pass through all rows
3569 : Int startRow, endRow, nRow;
3570 0 : if (row==-1)
3571 : {
3572 0 : nRow=vb.nRow();
3573 0 : startRow=0;
3574 0 : endRow=nRow-1;
3575 0 : modelVis.set(Complex(0.0,0.0));
3576 : }
3577 : else
3578 : {
3579 0 : nRow=1;
3580 0 : startRow=row;
3581 0 : endRow=row;
3582 0 : modelVis.xyPlane(row)=Complex(0.0,0.0);
3583 : }
3584 :
3585 0 : Int NAnt=0;
3586 :
3587 0 : if (doPointing)
3588 0 : NAnt = findPointingOffsets(vb,l_offsets, m_offsets,true);
3589 :
3590 :
3591 : //
3592 : // Get the uvws in a form that Fortran can use
3593 : //
3594 0 : Matrix<Double> uvw(3, vb.uvw().nelements());
3595 0 : uvw=0.0;
3596 0 : Vector<Double> dphase(vb.uvw().nelements());
3597 0 : dphase=0.0;
3598 : //
3599 : //NEGATING to correct for an image inversion problem
3600 : //
3601 0 : for (Int i=startRow;i<=endRow;i++)
3602 : {
3603 0 : for (Int idim=0;idim<2;idim++) uvw(idim,i)=-vb.uvw()(i)(idim);
3604 0 : uvw(2,i)=vb.uvw()(i)(2);
3605 : }
3606 :
3607 0 : rotateUVW(uvw, dphase, vb);
3608 0 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
3609 :
3610 : // This is the convention for dphase
3611 0 : dphase*=-1.0;
3612 :
3613 0 : Cube<Int> flags(vb.flagCube().shape());
3614 0 : flags=0;
3615 0 : flags(vb.flagCube())=true;
3616 :
3617 : //Check if ms has changed then cache new spw and chan selection
3618 0 : if(vb.newMS())
3619 0 : matchAllSpwChans(vb);
3620 :
3621 : //Channel matching for the actual spectral window of buffer
3622 0 : if(doConversion_p[vb.spectralWindow()])
3623 0 : matchChannel(vb.spectralWindow(), vb);
3624 : else
3625 : {
3626 0 : chanMap.resize();
3627 0 : chanMap=multiChanMap_p[vb.spectralWindow()];
3628 : }
3629 :
3630 0 : Vector<Int> rowFlags(vb.nRow());
3631 0 : rowFlags=0;
3632 0 : rowFlags(vb.flagRow())=true;
3633 0 : if(!usezero_p)
3634 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
3635 0 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
3636 :
3637 0 : IPosition s(modelVis.shape());
3638 0 : Int Conj=0,doGrad=0,ScanNo=0;
3639 0 : Double area=1.0;
3640 :
3641 0 : runFortranGet(uvw,dphase,vb.modelVisCube(),s,Conj,flags,rowFlags,row,
3642 0 : offset,&griddedVis,nx,ny,npol,nchan,vb,NAnt,ScanNo,sigma,
3643 0 : l_offsets,m_offsets,area,doGrad,PAIndex+1);
3644 0 : }
3645 : //
3646 : //---------------------------------------------------------------
3647 : //
3648 : // Finalize the FFT to the Sky. Here we actually do the FFT and
3649 : // return the resulting image
3650 0 : ImageInterface<Complex>& nPBWProjectFT::getImage(Matrix<Float>& weights,
3651 : Bool normalize)
3652 : {
3653 : //AlwaysAssert(lattice, AipsError);
3654 0 : AlwaysAssert(image, AipsError);
3655 :
3656 0 : logIO() << "#########getimage########" << LogIO::DEBUGGING << LogIO::POST;
3657 :
3658 0 : logIO() << LogOrigin("nPBWProjectFT", "getImage") << LogIO::NORMAL;
3659 :
3660 0 : weights.resize(sumWeight.shape());
3661 :
3662 0 : convertArray(weights, sumWeight);
3663 : //
3664 : // If the weights are all zero then we cannot normalize otherwise
3665 : // we don't care.
3666 : //
3667 0 : if(max(weights)==0.0)
3668 : {
3669 0 : if(normalize) logIO() << LogIO::SEVERE
3670 : << "No useful data in nPBWProjectFT: weights all zero"
3671 0 : << LogIO::POST;
3672 : else logIO() << LogIO::WARN << "No useful data in nPBWProjectFT: weights all zero"
3673 0 : << LogIO::POST;
3674 : }
3675 : else
3676 : {
3677 0 : const IPosition latticeShape = lattice->shape();
3678 :
3679 : logIO() << LogIO::DEBUGGING
3680 0 : << "Starting FFT and scaling of image" << LogIO::POST;
3681 : //
3682 : // x and y transforms (lattice has the gridded vis. Make the
3683 : // dirty images)
3684 : //
3685 0 : LatticeFFT::cfft2d(*lattice,false);
3686 :
3687 : //
3688 : // Apply the gridding correction
3689 : //
3690 : {
3691 0 : normalizeAvgPB();
3692 0 : Int inx = lattice->shape()(0);
3693 0 : Int iny = lattice->shape()(1);
3694 0 : Vector<Complex> correction(inx);
3695 :
3696 0 : Vector<Float> sincConv(nx);
3697 0 : Float centerX=nx/2;
3698 0 : for (Int ix=0;ix<nx;ix++)
3699 : {
3700 0 : Float x=C::pi*Float(ix-centerX)/(Float(nx)*Float(convSampling));
3701 0 : if(ix==centerX) sincConv(ix)=1.0;
3702 0 : else sincConv(ix)=sin(x)/x;
3703 : }
3704 :
3705 0 : IPosition cursorShape(4, inx, 1, 1, 1);
3706 0 : IPosition axisPath(4, 0, 1, 2, 3);
3707 0 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
3708 0 : LatticeIterator<Complex> lix(*lattice, lsx);
3709 :
3710 0 : LatticeStepper lavgpb(avgPB.shape(),cursorShape,axisPath);
3711 0 : LatticeIterator<Float> liavgpb(avgPB, lavgpb);
3712 :
3713 0 : for(lix.reset(),liavgpb.reset();
3714 0 : !lix.atEnd();
3715 0 : lix++,liavgpb++)
3716 : {
3717 0 : Int pol=lix.position()(2);
3718 0 : Int chan=lix.position()(3);
3719 :
3720 0 : if(weights(pol, chan)>0.0)
3721 : {
3722 : /*
3723 : Int iy=lix.position()(1);
3724 : gridder->correctX1D(correction,iy);
3725 :
3726 : Vector<Float> PBCorrection(liavgpb.rwVectorCursor().shape()),
3727 : avgPBVec(liavgpb.rwVectorCursor().shape());
3728 :
3729 : PBCorrection = liavgpb.rwVectorCursor();
3730 : avgPBVec = liavgpb.rwVectorCursor();
3731 :
3732 : for(int i=0;i<PBCorrection.shape();i++)
3733 : {
3734 : //
3735 : // This with the PS functions
3736 : //
3737 : // PBCorrection(i)=FUNC(avgPBVec(i))*sincConv(i)*sincConv(iy);
3738 : // if ((abs(PBCorrection(i)*correction(i))) >= pbLimit_p)
3739 : // lix.rwVectorCursor()(i) /= PBCorrection(i)*correction(i);
3740 : // else if (!makingPSF)
3741 : // lix.rwVectorCursor()(i) /= correction(i)*sincConv(i)*sincConv(iy);
3742 : //
3743 : // This without the PS functions
3744 : //
3745 : PBCorrection(i)=FUNC(avgPBVec(i))*sincConv(i)*sincConv(iy);
3746 : if ((abs(PBCorrection(i))) >= pbLimit_p)
3747 : lix.rwVectorCursor()(i) /= PBCorrection(i);
3748 : else if (!makingPSF)
3749 : lix.rwVectorCursor()(i) /= sincConv(i)*sincConv(iy);
3750 : }
3751 : */
3752 0 : if(normalize)
3753 : {
3754 0 : Complex rnorm(Float(inx)*Float(iny)/weights(pol,chan));
3755 0 : lix.rwCursor()*=rnorm;
3756 : }
3757 : else
3758 : {
3759 0 : Complex rnorm(Float(inx)*Float(iny));
3760 0 : lix.rwCursor()*=rnorm;
3761 : }
3762 : }
3763 : else
3764 0 : lix.woCursor()=0.0;
3765 : }
3766 0 : }
3767 :
3768 :
3769 0 : if(!isTiled)
3770 : {
3771 : //
3772 : // Check the section from the image BEFORE converting to a lattice
3773 : //
3774 0 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
3775 0 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
3776 0 : IPosition stride(4, 1);
3777 0 : IPosition trc(blc+image->shape()-stride);
3778 : //
3779 : // Do the copy
3780 : //
3781 0 : image->put(griddedData(blc, trc));
3782 0 : griddedData.resize(IPosition(1,0));
3783 0 : }
3784 0 : }
3785 :
3786 0 : return *image;
3787 : }
3788 : //
3789 : //---------------------------------------------------------------
3790 : //
3791 : // Get weight image
3792 0 : void nPBWProjectFT::getWeightImage(ImageInterface<Float>& weightImage,
3793 : Matrix<Float>& weights)
3794 : {
3795 :
3796 0 : logIO() << LogOrigin("nPBWProjectFT", "getWeightImage") << LogIO::NORMAL;
3797 :
3798 0 : weights.resize(sumWeight.shape());
3799 0 : convertArray(weights,sumWeight);
3800 :
3801 0 : const IPosition latticeShape = weightImage.shape();
3802 :
3803 0 : Int nx=latticeShape(0);
3804 0 : Int ny=latticeShape(1);
3805 :
3806 0 : IPosition loc(2, 0);
3807 0 : IPosition cursorShape(4, nx, ny, 1, 1);
3808 0 : IPosition axisPath(4, 0, 1, 2, 3);
3809 0 : LatticeStepper lsx(latticeShape, cursorShape, axisPath);
3810 0 : LatticeIterator<Float> lix(weightImage, lsx);
3811 0 : LatticeIterator<Float> liy(avgPB,lsx);
3812 0 : for(lix.reset();!lix.atEnd();lix++)
3813 : {
3814 : // UNUSED Int pol=lix.position()(2);
3815 : // UNUSED Int chan=lix.position()(3);
3816 : // lix.rwCursor()=weights(pol,chan);
3817 0 : lix.rwCursor()=liy.rwCursor();
3818 : }
3819 0 : }
3820 : //
3821 : //---------------------------------------------------------------
3822 : //
3823 0 : Bool nPBWProjectFT::toRecord(String&, //error,
3824 : RecordInterface& outRec, Bool withImage, const String diskimage)
3825 : {
3826 : // Save the current nPBWProjectFT object to an output state record
3827 0 : Bool retval = true;
3828 0 : Double cacheVal=(Double) cachesize;
3829 0 : outRec.define("cache", cacheVal);
3830 0 : outRec.define("tile", tilesize);
3831 :
3832 0 : Vector<Double> phaseValue(2);
3833 0 : String phaseUnit;
3834 0 : phaseValue=mTangent_p.getAngle().getValue();
3835 0 : phaseUnit= mTangent_p.getAngle().getUnit();
3836 0 : outRec.define("phasevalue", phaseValue);
3837 0 : outRec.define("phaseunit", phaseUnit);
3838 :
3839 0 : Vector<Double> dirValue(3);
3840 0 : String dirUnit;
3841 0 : dirValue=mLocation_p.get("m").getValue();
3842 0 : dirUnit=mLocation_p.get("m").getUnit();
3843 0 : outRec.define("dirvalue", dirValue);
3844 0 : outRec.define("dirunit", dirUnit);
3845 :
3846 0 : outRec.define("padding", padding_p);
3847 0 : outRec.define("maxdataval", maxAbsData);
3848 :
3849 0 : Vector<Int> center_loc(4), offset_loc(4);
3850 0 : for (Int k=0; k<4 ; k++)
3851 : {
3852 0 : center_loc(k)=centerLoc(k);
3853 0 : offset_loc(k)=offsetLoc(k);
3854 : }
3855 0 : outRec.define("centerloc", center_loc);
3856 0 : outRec.define("offsetloc", offset_loc);
3857 0 : outRec.define("sumofweights", sumWeight);
3858 0 : if(withImage && image ){
3859 0 : if(diskimage==""){
3860 0 : ImageInterface<Complex>& tempimage(*image);
3861 0 : Record imageContainer;
3862 0 : String error;
3863 0 : retval = (retval || tempimage.toRecord(error, imageContainer));
3864 0 : outRec.defineRecord("image", imageContainer);
3865 0 : }
3866 : else{
3867 0 : PagedImage<Complex> tempImage(image->shape(), image->coordinates(), diskimage);
3868 0 : tempImage.copyData(*image);
3869 0 : outRec.define("diskimage", diskimage);
3870 0 : }
3871 : }
3872 :
3873 0 : return retval;
3874 :
3875 : //
3876 0 : }
3877 : //---------------------------------------------------------------
3878 : //
3879 0 : Bool nPBWProjectFT::fromRecord(String&, //error,
3880 : const RecordInterface& inRec)
3881 : {
3882 0 : Bool retval = true;
3883 0 : imageCache=0; lattice=0; arrayLattice=0;
3884 : Double cacheVal;
3885 0 : inRec.get("cache", cacheVal);
3886 0 : cachesize=(Long) cacheVal;
3887 0 : inRec.get("tile", tilesize);
3888 :
3889 0 : Vector<Double> phaseValue(2);
3890 0 : inRec.get("phasevalue",phaseValue);
3891 0 : String phaseUnit;
3892 0 : inRec.get("phaseunit",phaseUnit);
3893 0 : Quantity val1(phaseValue(0), phaseUnit);
3894 0 : Quantity val2(phaseValue(1), phaseUnit);
3895 0 : MDirection phasecenter(val1, val2);
3896 :
3897 0 : mTangent_p=phasecenter;
3898 : // This should be passed down too but the tangent plane is
3899 : // expected to be specified in all meaningful cases.
3900 0 : tangentSpecified_p=true;
3901 0 : Vector<Double> dirValue(3);
3902 0 : String dirUnit;
3903 0 : inRec.get("dirvalue", dirValue);
3904 0 : inRec.get("dirunit", dirUnit);
3905 0 : MVPosition dummyMVPos(dirValue(0), dirValue(1), dirValue(2));
3906 0 : MPosition mLocation(dummyMVPos, MPosition::ITRF);
3907 0 : mLocation_p=mLocation;
3908 :
3909 0 : inRec.get("padding", padding_p);
3910 0 : inRec.get("maxdataval", maxAbsData);
3911 :
3912 0 : Vector<Int> center_loc(4), offset_loc(4);
3913 0 : inRec.get("centerloc", center_loc);
3914 0 : inRec.get("offsetloc", offset_loc);
3915 0 : uInt ndim4 = 4;
3916 0 : centerLoc=IPosition(ndim4, center_loc(0), center_loc(1), center_loc(2),
3917 0 : center_loc(3));
3918 0 : offsetLoc=IPosition(ndim4, offset_loc(0), offset_loc(1), offset_loc(2),
3919 0 : offset_loc(3));
3920 0 : inRec.get("sumofweights", sumWeight);
3921 0 : if(inRec.isDefined("image"))
3922 : {
3923 0 : Record imageAsRec=inRec.asRecord("image");
3924 0 : if(!image) image= new TempImage<Complex>();
3925 :
3926 0 : String error;
3927 0 : retval = (retval || image->fromRecord(error, imageAsRec));
3928 :
3929 : // Might be changing the shape of sumWeight
3930 0 : init();
3931 :
3932 0 : if(isTiled)
3933 0 : lattice=CountedPtr<Lattice<Complex> > (image, false);
3934 : else
3935 : {
3936 : //
3937 : // Make the grid the correct shape and turn it into an
3938 : // array lattice Check the section from the image BEFORE
3939 : // converting to a lattice
3940 : //
3941 0 : IPosition gridShape(4, nx, ny, npol, nchan);
3942 0 : griddedData.resize(gridShape);
3943 0 : griddedData=Complex(0.0);
3944 0 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
3945 0 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
3946 0 : IPosition start(4, 0);
3947 0 : IPosition stride(4, 1);
3948 0 : IPosition trc(blc+image->shape()-stride);
3949 0 : griddedData(blc, trc) = image->getSlice(start, image->shape());
3950 :
3951 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
3952 0 : lattice=arrayLattice;
3953 0 : }
3954 :
3955 0 : AlwaysAssert(image, AipsError);
3956 0 : }
3957 0 : else if(inRec.isDefined("diskimage")){
3958 0 : String diskim;
3959 0 : inRec.get("diskimage", diskim);
3960 0 : image=new PagedImage<Complex>(diskim);
3961 0 : }
3962 :
3963 :
3964 0 : return retval;
3965 0 : }
3966 : //
3967 : //---------------------------------------------------------------
3968 : //
3969 0 : void nPBWProjectFT::ok()
3970 : {
3971 0 : AlwaysAssert(image, AipsError);
3972 0 : }
3973 : //----------------------------------------------------------------------
3974 : //
3975 : // Make a plain straightforward honest-to-god image. This returns a
3976 : // complex image, without conversion to Stokes. The polarization
3977 : // representation is that required for the visibilities.
3978 : //
3979 0 : void nPBWProjectFT::makeImage(FTMachine::Type type,
3980 : VisSet& vs,
3981 : ImageInterface<Complex>& theImage,
3982 : Matrix<Float>& weight)
3983 : {
3984 0 : logIO() << LogOrigin("nPBWProjectFT", "makeImage") << LogIO::NORMAL;
3985 :
3986 0 : if(type==FTMachine::COVERAGE)
3987 : logIO() << "Type COVERAGE not defined for Fourier transforms"
3988 0 : << LogIO::EXCEPTION;
3989 :
3990 :
3991 : // Initialize the gradients
3992 0 : ROVisIter& vi(vs.iter());
3993 :
3994 : // Loop over all visibilities and pixels
3995 0 : VisBuffer vb(vi);
3996 :
3997 : // Initialize put (i.e. transform to Sky) for this model
3998 0 : vi.origin();
3999 :
4000 0 : if(vb.polFrame()==MSIter::Linear)
4001 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::LINEAR);
4002 : else
4003 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::CIRCULAR);
4004 :
4005 0 : initializeToSky(theImage,weight,vb);
4006 :
4007 : //
4008 : // Loop over the visibilities, putting VisBuffers
4009 : //
4010 0 : paChangeDetector.reset();
4011 :
4012 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk())
4013 : {
4014 0 : for (vi.origin(); vi.more(); vi++)
4015 : {
4016 0 : if (type==FTMachine::PSF) makingPSF=true;
4017 0 : findConvFunction(theImage,vb);
4018 :
4019 0 : switch(type)
4020 : {
4021 0 : case FTMachine::RESIDUAL:
4022 0 : vb.visCube()=vb.correctedVisCube();
4023 0 : vb.visCube()-=vb.modelVisCube();
4024 0 : put(vb, -1, false);
4025 0 : break;
4026 0 : case FTMachine::MODEL:
4027 0 : vb.visCube()=vb.modelVisCube();
4028 0 : put(vb, -1, false);
4029 0 : break;
4030 0 : case FTMachine::CORRECTED:
4031 0 : vb.visCube()=vb.correctedVisCube();
4032 0 : put(vb, -1, false);
4033 0 : break;
4034 0 : case FTMachine::PSF:
4035 0 : vb.visCube()=Complex(1.0,0.0);
4036 0 : makingPSF = true;
4037 0 : put(vb, -1, true);
4038 0 : break;
4039 0 : case FTMachine::OBSERVED:
4040 : default:
4041 0 : put(vb, -1, false);
4042 0 : break;
4043 : }
4044 : }
4045 : }
4046 0 : finalizeToSky();
4047 : // Normalize by dividing out weights, etc.
4048 0 : getImage(weight, true);
4049 0 : }
4050 :
4051 0 : void nPBWProjectFT::setPAIncrement(const Quantity& paIncrement)
4052 : {
4053 0 : paChangeDetector.setTolerance(paIncrement);
4054 0 : rotateAperture_p = true;
4055 0 : if (paIncrement.getValue("rad") < 0)
4056 0 : rotateAperture_p = false;
4057 0 : logIO() << "Setting PA increment to " << paIncrement.getValue("deg") << " deg" << endl;
4058 0 : }
4059 :
4060 0 : Bool nPBWProjectFT::verifyShapes(IPosition pbShape, IPosition skyShape)
4061 : {
4062 0 : if ((pbShape(0) != skyShape(0)) && // X-axis
4063 0 : (pbShape(1) != skyShape(1)) && // Y-axis
4064 0 : (pbShape(2) != skyShape(2))) // Poln-axis
4065 : {
4066 0 : throw(AipsError("Sky and/or polarization shape of the avgPB "
4067 0 : "and the sky model do not match."));
4068 : return false;
4069 : }
4070 0 : return true;
4071 :
4072 : }
4073 : //
4074 : // The functions are not azimuthally symmetric - so compute the
4075 : // support size carefully. Collect every PixInc pixel along a
4076 : // quarter circle of radius R (four-fold azimuthal symmetry is
4077 : // assumed), and check if any pixel in this list is above the
4078 : // threshold. If TRUE, then R is the support radius. Else decrease
4079 : // R and check again.
4080 : //
4081 : //
4082 0 : Bool nPBWProjectFT::findSupport(Array<Complex>& func, Float& threshold,Int& origin, Int& R)
4083 : {
4084 : Double NSteps;
4085 0 : Int PixInc=1;
4086 0 : Vector<Complex> vals;
4087 0 : IPosition ndx(4,origin,0,0,0);
4088 0 : Bool found=false;
4089 0 : IPosition cfShape=func.shape();
4090 0 : for(R=convSize/4;R>1;R--)
4091 : {
4092 0 : NSteps = 90*R/PixInc; //Check every PixInc pixel along a
4093 : //circle of radious R
4094 0 : vals.resize((Int)(NSteps+0.5));
4095 0 : vals=0;
4096 0 : for(Int th=0;th<NSteps;th++)
4097 : {
4098 0 : ndx(0)=(int)(origin + R*sin(2.0*M_PI*th*PixInc/R));
4099 0 : ndx(1)=(int)(origin + R*cos(2.0*M_PI*th*PixInc/R));
4100 :
4101 0 : if ((ndx(0) < cfShape(0)) && (ndx(1) < cfShape(1)))
4102 0 : vals(th)=convFunc(ndx);
4103 : }
4104 0 : if (max(abs(vals)) > threshold)
4105 0 : {found=true;break;}
4106 : }
4107 0 : return found;
4108 0 : }
4109 : /*
4110 : void nPBWProjectFT::makeAveragePB(const VisBuffer& vb,
4111 : const ImageInterface<Complex>& image,
4112 : Int& polInUse,
4113 : //TempImage<Float>& thesquintPB,
4114 : TempImage<Float>& theavgPB)
4115 : {
4116 : if (!resetPBs) return;
4117 :
4118 : Vector<Float> paList;
4119 : //
4120 : // Compute a list of PA present in the data with an increment of 10deg
4121 : // for PA-averaging of the average PB
4122 : //
4123 : {
4124 : Record time;
4125 : MSRange msr(*ms_p);
4126 : time = msr.range(MSS::TIMES);
4127 : Vector<Double> times=time.asArrayDouble("times");
4128 :
4129 : Float pa0=vb.feed_pa(times[0])[0],pa1, dpa;
4130 : paList.resize(1);
4131 :
4132 : paList[0]=pa0;
4133 :
4134 : for(uInt ipa=0;ipa<times.nelements();ipa++)
4135 : {
4136 : pa1=vb.feed_pa(times[ipa])[0];
4137 : dpa = abs(pa1-pa0);
4138 : if (dpa >= 10.0/57.2956)
4139 : {
4140 : pa0=pa1;
4141 : paList.resize(paList.nelements()+1,true);
4142 : paList[paList.nelements()-1] = pa0;
4143 : }
4144 : }
4145 : // paList.resize(1);
4146 : // paList[0]=pa0;
4147 : }
4148 : TempImage<Float> localPB;
4149 : localPB.setMaximumCacheSize(cachesize);
4150 : logIO() << LogOrigin("nPBWProjecFT","makeAveragePB")
4151 : << LogIO::NORMAL;
4152 :
4153 : localPB.resize(image.shape()); localPB.setCoordinateInfo(image.coordinates());
4154 : //
4155 : // If this is the first time, resize the average PB
4156 : //
4157 : if (resetPBs)
4158 : {
4159 : logIO() << LogOrigin("nPBWProjectFT", "makeAveragePB")
4160 : << "Computing the PA-averaged PBs. Averaging over "
4161 : << paList.nelements() << " PA values."
4162 : << LogIO::NORMAL
4163 : << LogIO::POST;
4164 : theavgPB.resize(localPB.shape());
4165 : theavgPB.setCoordinateInfo(localPB.coordinates());
4166 : theavgPB.set(0.0);
4167 : noOfPASteps = 0;
4168 : pbPeaks.resize(theavgPB.shape()(2));
4169 : pbPeaks.set(0.0);
4170 : resetPBs=false;
4171 : }
4172 : ProgressMeter pm(1.0, Double(paList.nelements()),
4173 : "Making average PB",
4174 : "", "", "", true);
4175 : //
4176 : // Make the Stokes PB
4177 : //
4178 : for(Int iPA=0;iPA<paList.nelements();iPA++)
4179 : {
4180 : localPB.set(1.0);
4181 : // vpSJ->applySquare(localPB,localPB,vb,1); // This is Stokes PB (not squinted)
4182 : {
4183 : VLACalcIlluminationConvFunc vlaPB;
4184 : vlaPB.setMaximumCacheSize(cachesize);
4185 : Vector<Float> pa(1);
4186 : pa[0] = paList(iPA);
4187 : // vlaPB.applyPB(localPB, vb, paList, bandID_p);
4188 : if (bandID_p == -1) bandID_p=getVisParams(vb);
4189 : vlaPB.applyPB(localPB, vb, pa, bandID_p);
4190 : }
4191 :
4192 : // IPosition twoDPBShape(localPB.shape());
4193 : IPosition twoDPBShape(localPB.shape());
4194 : // TempImage<Complex> localTwoDPB(twoDPBShape,localPB.coordinates());
4195 : // localTwoDPB.setMaximumCacheSize(cachesize);
4196 :
4197 : // Array<Float> raoff,decoff;
4198 : Float peak=0;
4199 : Int NAnt;
4200 : // NAnt=findPointingOffsets(vb,l_offsets,m_offsets,true);
4201 : // if (doPointing)
4202 : // NAnt=findPointingOffsets(vb,l_offsets,m_offsets,false);
4203 : // else
4204 : // NAnt = 1;
4205 :
4206 : // logIO() << LogOrigin("nPBWProjectFT", "makeAveragePB")
4207 : // << "Updating average PB @ PA = " << paList(iPA)*57.2956 << "deg"
4208 : // << LogIO::NORMAL
4209 : // << LogIO::POST;
4210 :
4211 : noOfPASteps++;
4212 : if (!doPointing) NAnt=1;
4213 : //
4214 : // For now, disable the "jittering" of the avg. PB by antenna pointing offsets.
4215 : // May not be required even for the very deep imaging - but thinking is a bit
4216 : // too clouded right now to delete the code below (SB)
4217 : //
4218 : NAnt=1;
4219 :
4220 : for(Int ant=0;ant<NAnt;ant++)
4221 : { //Ant loop
4222 : // {
4223 : // IPosition ndx(4,0,0,0,0);
4224 : // for(ndx(0)=0; ndx(0)<twoDPBShape(0); ndx(0)++)
4225 : // for(ndx(1)=0; ndx(1)<twoDPBShape(1); ndx(1)++)
4226 : // // for(ndx(2)=0; ndx(2)<polInUse; ndx(2)++)
4227 : // for(ndx(2)=0; ndx(2)<twoDPBShape(2); ndx(2)++)
4228 : // localTwoDPB.putAt(Complex((localPB(ndx)),0.0),ndx);
4229 : // }
4230 : //
4231 : // If antenna pointing errors are not applied, no shifting
4232 : // (which can be expensive) is required.
4233 : //
4234 : //
4235 : // Accumulate the shifted PBs
4236 : //
4237 : {
4238 : Bool isRefF,isRefC;
4239 : Array<Float> fbuf,cbuf;
4240 : // Array<Complex> cbuf;
4241 : isRefF=theavgPB.get(fbuf);
4242 : // isRefC=localTwoDPB.get(cbuf);
4243 : isRefC=localPB.get(cbuf);
4244 :
4245 : // Vector<Float> tmpPeaks(pbPeaks.nelements());
4246 : IPosition cs(cbuf.shape()),fs(fbuf.shape()),peakLoc(4,0,0,0,0);;
4247 : // if (makingPSF)
4248 : {
4249 : IPosition ndx(4,0,0,0,0),avgNDX(4,0,0,0,0);
4250 : for(ndx(2)=0,avgNDX(2)=0;ndx(2)<twoDPBShape(2);ndx(2)++,avgNDX(2)++)
4251 : {
4252 : for(ndx(0)=0,avgNDX(0)=0;ndx(0)<fs(0);ndx(0)++,avgNDX(0)++)
4253 : for(ndx(1)=0,avgNDX(1)=0;ndx(1)<fs(1);ndx(1)++,avgNDX(1)++)
4254 : {
4255 : Float val;
4256 : val = real(cbuf(ndx));
4257 : fbuf(avgNDX) += val;
4258 : if (fbuf(avgNDX) > peak) peak=fbuf(avgNDX);
4259 : }
4260 : }
4261 : }
4262 : if (!isRefF) theavgPB.put(fbuf);
4263 : // pbPeaks += tmpPeaks;
4264 : // cout << "Current PB peak = " << peak << endl;
4265 : }
4266 : }
4267 : pm.update(Double(iPA+1));
4268 : }
4269 : }
4270 : */
4271 :
4272 0 : void nPBWProjectFT::makeAntiAliasingOp(Vector<Complex>& op, const Int nx_l)
4273 : {
4274 0 : MathFunc<Float> sf(SPHEROIDAL);
4275 0 : if (op.nelements() != (uInt)nx_l)
4276 : {
4277 0 : op.resize(nx_l);
4278 0 : Int inner=nx_l/2, center=nx_l/2;
4279 : //
4280 : // The "complicated" equation below is worthy of a comment (as
4281 : // notes are called in program text).
4282 : //
4283 : // uvScale/nx == Size of a pixel size in the image in radians.
4284 : // Lets call it dx. HPBW is the HPBW for the antenna at the
4285 : // centre freq. in use. HPBW/dx == Pixel where the PB will be
4286 : // ~0.5x its peak value. ((2*N*HPBW)/dx) == the pixel where
4287 : // the N th. PB sidelobe will be (rougly speaking). When this
4288 : // value is equal to 3.0, the Spheroidal implemtation goes to
4289 : // zero!
4290 : //
4291 0 : Float dx=uvScale(0)*convSampling/nx;
4292 0 : Float MaxSideLobeNum = 3.0;
4293 0 : Float S=1.0*dx/(MaxSideLobeNum*2*HPBW),cfScale;
4294 :
4295 0 : cout << "UVSCALE = " << uvScale(0) << " " << convSampling << endl;
4296 0 : cout << "HPBW = " << HPBW
4297 0 : << " " << Diameter_p
4298 0 : << " " << uvScale(0)/nx
4299 0 : << " " << S
4300 0 : << " " << dx
4301 0 : << endl;
4302 :
4303 :
4304 0 : cfScale=S=6.0/inner;
4305 0 : for(Int ix=-inner;ix<inner;ix++) op(ix+center)=sf.value(abs((ix)*cfScale));
4306 : // for(Int ix=-inner;ix<inner;ix++)
4307 : // if (abs(op(ix+center)) > 1e-8)
4308 : // cout << "SF: " << ix
4309 : // << " " << (ix)*cfScale
4310 : // << " " << real(op(ix+center))
4311 : // << " " << imag(op(ix+center))
4312 : // << endl;
4313 : }
4314 0 : }
4315 :
4316 0 : void nPBWProjectFT::makeAntiAliasingCorrection(Vector<Complex>& correction,
4317 : const Vector<Complex>& op,
4318 : const Int nx_l)
4319 : {
4320 0 : if (correction.nelements() != (uInt)nx_l)
4321 : {
4322 0 : correction.resize(nx_l);
4323 0 : correction=0.0;
4324 0 : Int opLen=op.nelements(), orig=nx_l/2;
4325 0 : for(Int i=-opLen;i<opLen;i++)
4326 : {
4327 0 : correction(i+orig) += op(abs(i));
4328 : }
4329 0 : ArrayLattice<Complex> tmp(correction);
4330 0 : LatticeFFT::cfft(tmp,false);
4331 0 : correction=tmp.get();
4332 0 : }
4333 : // for(uInt i=0;i<correction.nelements();i++)
4334 : // cout << "FTSF: " << real(correction(i)) << " " << imag(correction(i)) << endl;
4335 0 : }
4336 :
4337 0 : void nPBWProjectFT::correctAntiAliasing(Lattice<Complex>& image)
4338 : {
4339 : // applyAntiAliasingOp(cf,2);
4340 0 : IPosition shape(image.shape());
4341 0 : IPosition ndx(shape);
4342 0 : ndx=0;
4343 0 : makeAntiAliasingCorrection(antiAliasingCorrection,
4344 0 : antiAliasingOp,shape(0));
4345 :
4346 0 : Complex tmp,val;
4347 0 : for(Int i=0;i<polInUse;i++)
4348 : {
4349 0 : ndx(2)=i;
4350 0 : for (Int iy=0;iy<shape(1);iy++)
4351 : {
4352 0 : for (Int ix=0;ix<shape(0);ix++)
4353 : {
4354 0 : ndx(0)=ix;
4355 0 : ndx(1)=iy;
4356 0 : tmp = image.getAt(ndx);
4357 0 : val=(antiAliasingCorrection(ix)*antiAliasingCorrection(iy));
4358 0 : if (abs(val) > 1e-5) tmp = tmp/val; else tmp=0.0;
4359 0 : image.putAt(tmp,ndx);
4360 : }
4361 : }
4362 : }
4363 0 : }
4364 :
4365 0 : void nPBWProjectFT::applyAntiAliasingOp(ImageInterface<Complex>& cf,
4366 : Vector<IPosition>& maxPos,
4367 : Int op,
4368 : Bool Square)
4369 : {
4370 : //
4371 : // First the spheroidal function
4372 : //
4373 0 : IPosition shape(cf.shape());
4374 0 : IPosition ndx(shape);
4375 0 : Vector<Double> refPixel=cf.coordinates().referencePixel();
4376 0 : ndx=0;
4377 0 : Int nx_l=shape(0),posX,posY;
4378 0 : makeAntiAliasingOp(antiAliasingOp, nx_l);
4379 :
4380 0 : Complex tmp,gain;
4381 :
4382 0 : for(Int i=0;i<polInUse;i++)
4383 : {
4384 0 : ndx(2)=i;
4385 0 : for (Int iy=-nx_l/2;iy<nx_l/2;iy++)
4386 : {
4387 0 : for (Int ix=-nx_l/2;ix<nx_l/2;ix++)
4388 : {
4389 0 : ndx(0)=ix+nx_l/2;
4390 0 : ndx(1)=iy+nx_l/2;
4391 0 : tmp = cf.getAt(ndx);
4392 0 : posX=ndx(0)+(Int)(refPixel(0)-maxPos(i)(0))+1;
4393 0 : posY=ndx(1)+(Int)(refPixel(1)-maxPos(i)(1))+1;
4394 0 : if ((posX > 0) && (posX < nx_l) &&
4395 0 : (posY > 0) && (posY < nx_l))
4396 0 : gain = antiAliasingOp(posX)*antiAliasingOp(posY);
4397 : else
4398 0 : if (op==2) gain = 1.0; else gain=0.0;
4399 0 : if (Square) gain *= gain;
4400 0 : switch (op)
4401 : {
4402 0 : case 0: tmp = tmp+gain;break;
4403 0 : case 1:
4404 : {
4405 0 : tmp = tmp*gain;
4406 0 : break;
4407 : }
4408 0 : case 2: tmp = tmp/gain;break;
4409 : }
4410 0 : cf.putAt(tmp,ndx);
4411 : }
4412 : }
4413 : }
4414 0 : };
4415 0 : void nPBWProjectFT::applyAntiAliasingOp(ImageInterface<Float>& cf,
4416 : Vector<IPosition>& maxPos,
4417 : Int op,
4418 : Bool Square)
4419 : {
4420 : //
4421 : // First the spheroidal function
4422 : //
4423 0 : IPosition shape(cf.shape());
4424 0 : IPosition ndx(shape);
4425 0 : Vector<Double> refPixel=cf.coordinates().referencePixel();
4426 0 : ndx=0;
4427 0 : Int nx_l=shape(0),posX,posY;
4428 0 : makeAntiAliasingOp(antiAliasingOp, nx_l);
4429 :
4430 : Float tmp,gain;
4431 :
4432 0 : for(Int i=0;i<polInUse;i++)
4433 : {
4434 0 : ndx(2)=i;
4435 0 : for (Int iy=-nx_l/2;iy<nx_l/2;iy++)
4436 : {
4437 0 : for (Int ix=-nx_l/2;ix<nx_l/2;ix++)
4438 : {
4439 0 : ndx(0)=ix+nx_l/2;
4440 0 : ndx(1)=iy+nx_l/2;
4441 0 : tmp = cf.getAt(ndx);
4442 0 : posX=ndx(0)+(Int)(refPixel(0)-maxPos(i)(0))+1;
4443 0 : posY=ndx(1)+(Int)(refPixel(1)-maxPos(i)(1))+1;
4444 0 : if ((posX > 0) && (posX < nx_l) &&
4445 0 : (posY > 0) && (posY < nx_l))
4446 0 : gain = real(antiAliasingOp(posX)*antiAliasingOp(posY));
4447 : else
4448 0 : if (op==2) gain = 1.0; else gain=0.0;
4449 0 : if (Square) gain *= gain;
4450 0 : switch (op)
4451 : {
4452 0 : case 0: tmp = tmp+gain;break;
4453 0 : case 1:
4454 : {
4455 0 : tmp = tmp*gain;
4456 0 : break;
4457 : }
4458 0 : case 2: tmp = tmp/gain;break;
4459 : }
4460 0 : cf.putAt(tmp,ndx);
4461 : }
4462 : }
4463 : }
4464 0 : };
4465 : } //# NAMESPACE CASA - END
4466 :
|
