Line data Source code
1 : // -*- C++ -*-
2 : //# AWProjectFT.cc: Implementation of AWProjectFT 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 <casacore/casa/Quanta/UnitMap.h>
30 : #include <casacore/casa/Quanta/MVTime.h>
31 : #include <casacore/casa/Quanta/UnitVal.h>
32 : #include <casacore/casa/Containers/Block.h>
33 : #include <casacore/casa/Containers/Record.h>
34 : #include <casacore/casa/Arrays/Array.h>
35 : #include <casacore/casa/OS/HostInfo.h>
36 : #include <sstream>
37 :
38 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
39 : #include <casacore/images/Images/ImageInterface.h>
40 :
41 : #include <synthesis/TransformMachines/StokesImageUtil.h>
42 : #include <synthesis/TransformMachines/SynthesisError.h>
43 : #include <synthesis/TransformMachines/AWProjectFT.h>
44 : #include <synthesis/TransformMachines/AWConvFunc.h>
45 : #include <synthesis/TransformMachines/CFStore2.h>
46 : #include <synthesis/MeasurementComponents/ExpCache.h>
47 : #include <synthesis/MeasurementComponents/CExp.h>
48 : #include <synthesis/TransformMachines/AWVisResampler.h>
49 : #include <synthesis/TransformMachines/VBStore.h>
50 :
51 : #include <casacore/scimath/Mathematics/FFTServer.h>
52 : #include <casacore/scimath/Mathematics/MathFunc.h>
53 : #include <casacore/measures/Measures/MeasTable.h>
54 : #include <iostream>
55 : #include <casacore/casa/OS/Timer.h>
56 :
57 : #include <synthesis/TransformMachines/ATerm.h>
58 : #include <synthesis/TransformMachines/NoOpATerm.h>
59 : #include <synthesis/TransformMachines/AWConvFunc.h>
60 : #include <synthesis/TransformMachines/EVLAAperture.h>
61 : #include <synthesis/TransformMachines/AWConvFuncEPJones.h>
62 :
63 : //#define CONVSIZE (1024*2)
64 : // #define OVERSAMPLING 2
65 : #define USETABLES 0 // If equal to 1, use tabulated exp() and
66 : // complex exp() functions.
67 : #define MAXPOINTINGERROR 250.0 // Max. pointing error in arcsec used to
68 : // determine the resolution of the
69 : // tabulated exp() function.
70 : #define DORES true
71 :
72 :
73 : using namespace casacore;
74 : namespace casa { //# NAMESPACE CASA - BEGIN
75 :
76 : #define NEED_UNDERSCORES
77 : extern "C"
78 : {
79 : //
80 : // The Gridding Convolution Function (GCF) used by the underlying
81 : // gridder written in FORTRAN.
82 : //
83 : // The arguments must all be pointers and the value of the GCF at
84 : // the given (u,v) point is returned in the weight variable. Making
85 : // this a function which returns a complex value (namely the weight)
86 : // has problems when called in FORTRAN - I (SB) don't understand
87 : // why.
88 : //
89 : #if defined(NEED_UNDERSCORES)
90 : #define nwcppeij nwcppeij_
91 : #endif
92 : //
93 : //---------------------------------------------------------------
94 : //
95 : IlluminationConvFunc awEij;
96 0 : void awcppeij(Double *griduvw, Double *area,
97 : Double *raoff1, Double *decoff1,
98 : Double *raoff2, Double *decoff2,
99 : Int *doGrad,
100 : Complex *weight,
101 : Complex *dweight1,
102 : Complex *dweight2,
103 : Double *currentCFPA)
104 : {
105 0 : Complex w,d1,d2;
106 0 : awEij.getValue(griduvw, raoff1, raoff2, decoff1, decoff2,
107 : area,doGrad,w,d1,d2,*currentCFPA);
108 0 : *weight = w;
109 0 : *dweight1 = d1;
110 0 : *dweight2 = d2;
111 0 : }
112 : }
113 :
114 0 : ATerm* AWProjectFT::createTelescopeATerm(const String& telescopeName, const Bool& isATermOn)
115 : {
116 0 : LogIO os(LogOrigin("AWProjectFT", "createTelescopeATerm",WHERE));
117 :
118 0 : if (!isATermOn) return new NoOpATerm();
119 :
120 : // MSObservationColumns msoc(ms.observation());
121 : // String ObsName=msoc.telescopeName()(0);
122 0 : if ((telescopeName == "EVLA") || (telescopeName == "VLA"))
123 0 : return new EVLAAperture();
124 : else
125 : {
126 0 : os << "Telescope name ('"+
127 0 : telescopeName+"') in the MS not recognized to create the telescope specific ATerm"
128 0 : << LogIO::WARN;
129 : }
130 :
131 0 : return NULL;
132 0 : }
133 :
134 0 : CountedPtr<ConvolutionFunction> AWProjectFT::makeCFObject(const String& telescopeName,
135 : const Bool aTermOn,
136 : const Bool psTermOn,
137 : const Bool wTermOn,
138 : const Bool mTermOn,
139 : const Bool wbAWP,
140 : const Bool conjBeams)
141 : {
142 0 : LogIO log_l(LogOrigin("AWProjectFT", "makeCFObject"));
143 0 : if ((psTermOn==false) && (aTermOn==false))
144 0 : log_l << "Both, psterm and aterm cannot be False. Please set one or both to True." << LogIO::EXCEPTION;
145 :
146 0 : CountedPtr<ATerm> apertureFunction = AWProjectFT::createTelescopeATerm(telescopeName, aTermOn);
147 0 : CountedPtr<PSTerm> psTerm = new PSTerm();
148 0 : CountedPtr<WTerm> wTerm = new WTerm();
149 :
150 : //
151 : // Selectively switch off CFTerms.
152 : //
153 0 : if (aTermOn == false) {apertureFunction->setOpCode(CFTerms::NOOP);}
154 0 : if (psTermOn == false) psTerm->setOpCode(CFTerms::NOOP);
155 0 : if (wTermOn == False) wTerm->setOpCode(CFTerms::NOOP);
156 : //
157 : // Construct the CF object with appropriate CFTerms.
158 : //
159 0 : CountedPtr<ConvolutionFunction> awConvFunc;
160 : // awConvFunc = new AWConvFunc(apertureFunction,psTerm,wTerm, !wbAWP);
161 : //if ((ftmName=="mawprojectft") || (mTermOn))
162 0 : if (mTermOn)
163 0 : awConvFunc = new AWConvFuncEPJones(apertureFunction,psTerm,wTerm,wbAWP,conjBeams);
164 : else
165 0 : awConvFunc = new AWConvFunc(apertureFunction,psTerm,wTerm,wbAWP, conjBeams);
166 0 : return awConvFunc;
167 0 : }
168 : //
169 : //---------------------------------------------------------------
170 : //
171 0 : AWProjectFT::AWProjectFT()
172 0 : : FTMachine(), padding_p(1.0), nWPlanes_p(1),
173 0 : imageCache(0), cachesize(0), tilesize(16),
174 0 : gridder(0), isTiled(false), arrayLattice( ), lattice( ),
175 0 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
176 0 : pointingToImage(0), usezero_p(false),
177 : // convFunc_p(), convWeights_p(),
178 0 : epJ_p(),
179 0 : doPBCorrection(true), conjBeams_p(true),/*cfCache_p(cfcache),*/ paChangeDetector(),
180 0 : rotateOTFPAIncr_p(0.1),
181 0 : Second("s"),Radian("rad"),Day("d"), pbNormalized_p(false), paNdxProcessed_p(),
182 0 : visResampler_p(), sensitivityPatternQualifier_p(-1),sensitivityPatternQualifierStr_p(""),
183 0 : rotatedConvFunc_p(),//cfs2_p(), cfwts2_p(),
184 0 : runTime1_p(0.0), previousSPWID_p(-1)
185 : {
186 : // convSize=0;
187 0 : tangentSpecified_p=false;
188 0 : lastIndex_p=0;
189 0 : paChangeDetector.reset();
190 0 : pbLimit_p=5e-2;
191 : //
192 : // Get various parameters from the visibilities.
193 : //
194 0 : doPointing=1;
195 :
196 0 : maxConvSupport=-1;
197 : //
198 : // Set up the Conv. Func. disk cache manager object.
199 : //
200 : // if (!cfCache_p.null()) delete &cfCache_p;
201 : // cfCache_p=cfcache;
202 0 : convSampling=OVERSAMPLING;
203 : //convSize=CONVSIZE;
204 0 : Long hostRAM = (HostInfo::memoryTotal(true)*1024); // In bytes
205 0 : hostRAM = hostRAM/(sizeof(Float)*2); // In complex pixels
206 0 : if (cachesize > hostRAM) cachesize=hostRAM;
207 0 : sigma=1.0;
208 0 : canComputeResiduals_p=DORES;
209 : // cfs2_p = &cfCache_p->memCache2_p[0];//new CFStore2;
210 : // cfwts2_p = &cfCache_p->memCacheWt2_p[0];//new CFStore2;
211 0 : pop_p->init();
212 0 : CFBuffer::initCFBStruct(cfbst_pub);
213 : // rotatedConvFunc_p.data=new Array<Complex>();
214 0 : }
215 : //
216 : //---------------------------------------------------------------
217 : //
218 0 : AWProjectFT::AWProjectFT(Int nWPlanes, Long icachesize,
219 : CountedPtr<CFCache>& cfcache,
220 : CountedPtr<ConvolutionFunction>& cf,
221 : CountedPtr<VisibilityResamplerBase>& visResampler,
222 : Bool applyPointingOffset,
223 : Bool doPBCorr,
224 : Int itilesize,
225 : Float pbLimit,
226 : Bool usezero,
227 : Bool conjBeams,
228 : Bool doublePrecGrid,
229 0 : PolOuterProduct::MuellerType muellerType)
230 0 : : FTMachine(cfcache,cf), padding_p(1.0), nWPlanes_p(nWPlanes),
231 0 : imageCache(0), cachesize(icachesize), tilesize(itilesize),
232 0 : gridder(0), isTiled(false), arrayLattice( ), lattice( ),
233 0 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
234 0 : pointingToImage(0), usezero_p(usezero),
235 : // convFunc_p(), convWeights_p(),
236 0 : epJ_p(),
237 0 : doPBCorrection(doPBCorr), conjBeams_p(conjBeams),
238 0 : /*cfCache_p(cfcache),*/ paChangeDetector(),
239 0 : rotateOTFPAIncr_p(0.1),
240 0 : Second("s"),Radian("rad"),Day("d"), pbNormalized_p(false),
241 0 : visResampler_p(visResampler), sensitivityPatternQualifier_p(-1),sensitivityPatternQualifierStr_p(""),
242 0 : rotatedConvFunc_p(), runTime1_p(0.0), previousSPWID_p(-1)
243 : {
244 : //convSize=0;
245 0 : tangentSpecified_p=false;
246 0 : lastIndex_p=0;
247 0 : paChangeDetector.reset();
248 0 : pbLimit_p=pbLimit;
249 : //
250 : // Get various parameters from the visibilities.
251 : //
252 0 : if (applyPointingOffset) doPointing=1; else doPointing=0;
253 :
254 0 : maxConvSupport=-1;
255 : //
256 : // Set up the Conv. Func. disk cache manager object.
257 : //
258 : // if (!cfCache_p.null()) delete &cfCache_p;
259 : // cfCache_p=cfcache;
260 0 : convSampling=OVERSAMPLING;
261 : //convSize=CONVSIZE;
262 0 : Long hostRAM = (HostInfo::memoryTotal(true)*1024); // In bytes
263 0 : hostRAM = hostRAM/(sizeof(Float)*2); // In complex pixels
264 0 : if (cachesize > hostRAM) cachesize=hostRAM;
265 0 : sigma=1.0;
266 0 : canComputeResiduals_p=DORES;
267 0 : if (!cfCache_p.null())
268 : {
269 0 : cfs2_p = CountedPtr<CFStore2>(&(cfCache_p->memCache2_p)[0],false);//new CFStore2;
270 0 : cfwts2_p = CountedPtr<CFStore2>(&cfCache_p->memCacheWt2_p[0],false);//new CFStore2;
271 : }
272 0 : pop_p->init();
273 0 : useDoubleGrid_p=doublePrecGrid;
274 0 : useDoubleGrid_p=doublePrecGrid;
275 0 : CFBuffer::initCFBStruct(cfbst_pub);
276 0 : muellerType_p = muellerType;
277 0 : }
278 : //
279 : //---------------------------------------------------------------
280 : //
281 0 : AWProjectFT::AWProjectFT(const RecordInterface& stateRec)
282 0 : : FTMachine(),Second("s"),Radian("rad"),Day("d"),visResampler_p()//, cfs2_p(), cfwts2_p()
283 : {
284 0 : LogIO log_l(LogOrigin("AWProjectFT", "AWProjectFT[R&D]"));
285 : //
286 : // Construct from the input state record
287 : //
288 0 : String error;
289 :
290 0 : if (!fromRecord(stateRec)) {
291 0 : log_l << "Failed to create " << name() << " object." << LogIO::EXCEPTION;
292 : };
293 0 : maxConvSupport=-1;
294 0 : convSampling=OVERSAMPLING;
295 0 : visResampler_p->init(useDoubleGrid_p);
296 : //convSize=CONVSIZE;
297 0 : canComputeResiduals_p=DORES;
298 0 : if (!cfCache_p.null())
299 : {
300 0 : cfs2_p = CountedPtr<CFStore2>(&cfCache_p->memCache2_p[0],false);//new CFStore2;
301 0 : cfwts2_p = CountedPtr<CFStore2>(&cfCache_p->memCacheWt2_p[0],false);//new CFStore2;
302 : }
303 0 : pop_p->init();
304 0 : }
305 : //
306 : //----------------------------------------------------------------------
307 : //
308 0 : AWProjectFT::AWProjectFT(const AWProjectFT& other):FTMachine()
309 : {
310 0 : operator=(other);
311 0 : }
312 : //
313 : //---------------------------------------------------------------
314 : //
315 : // This is nasty, we should use CountedPointers here.
316 0 : AWProjectFT::~AWProjectFT()
317 : {
318 0 : if(imageCache) delete imageCache; imageCache=0;
319 0 : if(gridder) delete gridder; gridder=0;
320 0 : }
321 : //
322 : //---------------------------------------------------------------
323 : //
324 0 : AWProjectFT& AWProjectFT::operator=(const AWProjectFT& other)
325 : {
326 0 : if(this!=&other)
327 : {
328 : //Do the base parameters
329 0 : FTMachine::operator=(other);
330 :
331 :
332 0 : padding_p=other.padding_p;
333 :
334 0 : nWPlanes_p=other.nWPlanes_p;
335 0 : imageCache=other.imageCache;
336 0 : cachesize=other.cachesize;
337 0 : tilesize=other.tilesize;
338 0 : cfRefFreq_p = other.cfRefFreq_p;
339 0 : if(other.gridder==0) gridder=0;
340 : else
341 : {
342 0 : uvScale.resize();
343 0 : uvOffset.resize();
344 0 : uvScale=other.uvScale;
345 0 : uvOffset=other.uvOffset;
346 0 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
347 0 : uvScale, uvOffset,
348 0 : "SF");
349 : }
350 :
351 0 : isTiled=other.isTiled;
352 0 : lattice=0;
353 0 : arrayLattice=0;
354 :
355 0 : maxAbsData=other.maxAbsData;
356 0 : centerLoc=other.centerLoc;
357 0 : offsetLoc=other.offsetLoc;
358 0 : pointingToImage=other.pointingToImage;
359 0 : usezero_p=other.usezero_p;
360 :
361 :
362 0 : padding_p=other.padding_p;
363 0 : nWPlanes_p=other.nWPlanes_p;
364 0 : imageCache=other.imageCache;
365 0 : cachesize=other.cachesize;
366 0 : tilesize=other.tilesize;
367 0 : isTiled=other.isTiled;
368 0 : maxAbsData=other.maxAbsData;
369 0 : centerLoc=other.centerLoc;
370 0 : offsetLoc=other.offsetLoc;
371 0 : pointingToImage=other.pointingToImage;
372 0 : usezero_p=other.usezero_p;
373 0 : doPBCorrection = other.doPBCorrection;
374 0 : maxConvSupport= other.maxConvSupport;
375 :
376 0 : epJ_p=other.epJ_p;
377 : //convSize=other.convSize;
378 0 : lastIndex_p=other.lastIndex_p;
379 0 : paChangeDetector=other.paChangeDetector;
380 0 : pbLimit_p=other.pbLimit_p;
381 : //
382 : // Get various parameters from the visibilities.
383 : //
384 0 : doPointing=other.doPointing;
385 :
386 0 : maxConvSupport=other.maxConvSupport;
387 : //
388 : // Set up the Conv. Func. disk cache manager object.
389 : //
390 0 : cfCache_p=other.cfCache_p;
391 0 : convSampling=other.convSampling;
392 : //convSize=other.convSize;
393 0 : cachesize=other.cachesize;
394 :
395 0 : currentCFPA=other.currentCFPA;
396 0 : lastPAUsedForWtImg = other.lastPAUsedForWtImg;
397 0 : avgPB_p = other.avgPB_p;
398 0 : avgPBSq_p = other.avgPBSq_p;
399 0 : convFuncCtor_p = other.convFuncCtor_p;
400 0 : pbNormalized_p = other.pbNormalized_p;
401 0 : sensitivityPatternQualifier_p = other.sensitivityPatternQualifier_p;
402 0 : sensitivityPatternQualifierStr_p = other.sensitivityPatternQualifierStr_p;
403 0 : visResampler_p=other.visResampler_p; // Copy the counted pointer
404 : // visResampler_p=other.visResampler_p->clone();
405 : // *visResampler_p = *other.visResampler_p; // Call the appropriate operator=()
406 :
407 0 : rotatedConvFunc_p = other.rotatedConvFunc_p;
408 0 : cfs2_p = other.cfs2_p;
409 0 : cfwts2_p = other.cfwts2_p;
410 0 : paNdxProcessed_p = other.paNdxProcessed_p;
411 0 : imRefFreq_p = other.imRefFreq_p;
412 0 : conjBeams_p = other.conjBeams_p;
413 0 : rotateOTFPAIncr_p=other.rotateOTFPAIncr_p;
414 0 : computePAIncr_p=other.computePAIncr_p;
415 0 : runTime1_p = other.runTime1_p;
416 0 : muellerType_p = other.muellerType_p;
417 0 : previousSPWID_p = other.previousSPWID_p;
418 : };
419 0 : return *this;
420 : };
421 : //
422 : //----------------------------------------------------------------------
423 : //
424 0 : void AWProjectFT::init()
425 : {
426 0 : LogIO log_l(LogOrigin("AWProjectFT", "init[R&D]"));
427 :
428 0 : nx = image->shape()(0);
429 0 : ny = image->shape()(1);
430 0 : npol = image->shape()(2);
431 0 : nchan = image->shape()(3);
432 :
433 0 : if(image->shape().product()>cachesize)
434 0 : isTiled=true;
435 : else
436 0 : isTiled=false;
437 :
438 0 : sumWeight.resize(npol, nchan);
439 0 : sumCFWeight.resize(npol, nchan);
440 :
441 0 : wConvSize=max(1, nWPlanes_p);
442 :
443 0 : CoordinateSystem cs=image->coordinates();
444 0 : uvScale.resize(3);
445 0 : uvScale=0.0;
446 0 : uvScale(0)=Float(nx)*cs.increment()(0);
447 0 : uvScale(1)=Float(ny)*cs.increment()(1);
448 0 : uvScale(2)=Float(wConvSize)*abs(cs.increment()(0));
449 :
450 0 : Int index= cs.findCoordinate(Coordinate::SPECTRAL);
451 0 : SpectralCoordinate spCS = cs.spectralCoordinate(index);
452 0 : imRefFreq_p = spCS.referenceValue()(0);
453 :
454 0 : uvOffset.resize(3);
455 0 : uvOffset(0)=nx/2;
456 0 : uvOffset(1)=ny/2;
457 0 : uvOffset(2)=0;
458 :
459 0 : if(gridder) delete gridder; gridder=0;
460 0 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
461 0 : uvScale, uvOffset,
462 0 : "SF");
463 :
464 : // Set up image cache needed for gridding.
465 0 : if(imageCache) delete imageCache; imageCache=0;
466 :
467 : // The tile size should be large enough that the
468 : // extended convolution function can fit easily
469 0 : if(isTiled)
470 : {
471 0 : Float tileOverlap=0.5;
472 0 : tilesize=min(256,tilesize);
473 0 : IPosition tileShape=IPosition(4,tilesize,tilesize,npol,nchan);
474 0 : Vector<Float> tileOverlapVec(4);
475 0 : tileOverlapVec=0.0;
476 0 : tileOverlapVec(0)=tileOverlap;
477 0 : tileOverlapVec(1)=tileOverlap;
478 : if (sizeof(long) < 4) // 32-bit machine
479 : {
480 : Int tmpCacheVal=static_cast<Int>(cachesize);
481 : imageCache=new LatticeCache <Complex> (*image, tmpCacheVal, tileShape,
482 : tileOverlapVec,
483 : (tileOverlap>0.0));
484 : }
485 : else // 64-bit machine
486 : {
487 0 : Long tmpCacheVal=cachesize;
488 0 : imageCache=new LatticeCache <Complex> (*image, tmpCacheVal, tileShape,
489 : tileOverlapVec,
490 0 : (tileOverlap>0.0));
491 : }
492 0 : }
493 :
494 : #if(USETABLES)
495 : Double StepSize;
496 : Int N=500000;
497 : StepSize = abs((((2*nx)/uvScale(0))/(sigma) +
498 : MAXPOINTINGERROR*1.745329E-02*(sigma)/3600.0))/N;
499 : if (!awEij.isReady())
500 : {
501 : log_l << "Making lookup table for exp function with a resolution of "
502 : << StepSize << " radians. "
503 : << "Memory used: " << sizeof(Float)*N/(1024.0*1024.0)<< " MB."
504 : << LogIO::NORMAL
505 : <<LogIO::POST;
506 :
507 : awEij.setSigma(sigma);
508 : awEij.initExpTable(N,StepSize);
509 : // ExpTab.build(N,StepSize);
510 :
511 : log_l << "Making lookup table for complex exp function with a resolution of "
512 : << 2*M_PI/N << " radians. "
513 : << "Memory used: " << 2*sizeof(Float)*N/(1024.0*1024.0) << " MB."
514 : << LogIO::NORMAL
515 : << LogIO::POST;
516 : awEij.initCExpTable(N);
517 : // CExpTab.build(N);
518 : }
519 : #endif
520 : // vpSJ->reset();
521 0 : paChangeDetector.reset();
522 0 : makingPSF = false;
523 :
524 : //cerr << "Current runTime = " << runTime << endl;
525 0 : }
526 : //
527 : //---------------------------------------------------------------
528 : //
529 0 : MDirection::Convert AWProjectFT::makeCoordinateMachine(const VisBuffer& vb,
530 : const MDirection::Types& From,
531 : const MDirection::Types& To,
532 : MEpoch& last)
533 : {
534 0 : LogIO log_l(LogOrigin("AWProjectFT","makeCoordinateMachine[R&D]"));
535 0 : Double time = getCurrentTimeStamp(vb);
536 :
537 0 : MEpoch epoch(Quantity(time,Second),MEpoch::TAI);
538 : // epoch = MEpoch(Quantity(time,Second),MEpoch::TAI);
539 : //
540 : // ...now make an object to hold the observatory position info...
541 : //
542 0 : MPosition pos;
543 0 : String ObsName=vb.msColumns().observation().telescopeName()(vb.arrayId());
544 :
545 0 : if (!MeasTable::Observatory(pos,ObsName))
546 0 : log_l << "Observatory position for "+ ObsName + " not found"
547 0 : << LogIO::EXCEPTION;
548 : //
549 : // ...now make a Frame object out of the observatory position and
550 : // time objects...
551 : //
552 0 : MeasFrame frame(epoch,pos);
553 : //
554 : // ...finally make the convert machine.
555 : //
556 0 : MDirection::Convert mac(MDirection::Ref(From,frame),
557 0 : MDirection::Ref(To,frame));
558 :
559 0 : MEpoch::Convert toLAST = MEpoch::Convert(MEpoch::Ref(MEpoch::TAI,frame),
560 0 : MEpoch::Ref(MEpoch::LAST,frame));
561 0 : last = toLAST(epoch);
562 :
563 0 : return mac;
564 0 : }
565 : //
566 : //---------------------------------------------------------------
567 : //
568 0 : int AWProjectFT::findPointingOffsets(const VisBuffer& vb,
569 : Array<Float> &l_off,
570 : Array<Float> &m_off,
571 : Bool Evaluate)
572 : {
573 0 : LogIO log_l(LogOrigin("AWProjectFT", "findPointingOffsets[R&D]"));
574 0 : Int NAnt = 0;
575 0 : MEpoch LAST;
576 0 : Double thisTime = getCurrentTimeStamp(vb);
577 : // Array<Float> pointingOffsets = epJ->nearest(thisTime);
578 0 : if (epJ_p.null()) return 0;
579 0 : Array<Float> pointingOffsets; epJ_p->nearest(thisTime,pointingOffsets);
580 0 : NAnt=pointingOffsets.shape()(2);
581 0 : l_off.resize(IPosition(3,1,1,NAnt)); // Poln x NChan x NAnt
582 0 : m_off.resize(IPosition(3,1,1,NAnt)); // Poln x NChan x NAnt
583 : // Can't figure out how to do the damn slicing of [Pol,NChan,NAnt,1] array
584 : // into [Pol,NChan,NAnt] array
585 : //
586 : // l_off = pointingOffsets(Slicer(IPosition(4,0,0,0,0),IPosition(4,1,1,NAnt+1,0)));
587 : // m_off = pointingOffsets(Slicer(IPosition(4,1,0,0,0),IPosition(4,1,1,NAnt+1,0)));
588 0 : IPosition tndx(3,0,0,0), sndx(4,0,0,0,0);
589 0 : for(tndx(2)=0;tndx(2)<NAnt; tndx(2)++,sndx(2)++)
590 : // for(Int j=0;j<NAnt;j++)
591 : {
592 : // l_off(IPosition(3,0,0,j)) = pointingOffsets(IPosition(4,0,0,j,0));
593 : // m_off(IPosition(3,0,0,j)) = pointingOffsets(IPosition(4,1,0,j,0));
594 :
595 0 : sndx(0)=0; l_off(tndx) = pointingOffsets(sndx);
596 0 : sndx(0)=2; m_off(tndx) = pointingOffsets(sndx);
597 : }
598 0 : return NAnt;
599 : if (!Evaluate) return NAnt;
600 :
601 : // cout << "AzEl Offsets: " << pointingOffsets << endl;
602 : //
603 : // Make a Coordinate Conversion Machine to go from (Az,El) to
604 : // (HA,Dec).
605 : //
606 : MDirection::Convert toAzEl = makeCoordinateMachine(vb,MDirection::HADEC,
607 : MDirection::AZEL,
608 : LAST);
609 : MDirection::Convert toHADec = makeCoordinateMachine(vb,MDirection::AZEL,
610 : MDirection::HADEC,
611 : LAST);
612 : //
613 : // ...and now hope that it all works and works correctly!!!
614 : //
615 : Quantity dAz(0,Radian),dEl(0,Radian);
616 : //
617 : // An array of shape [2,1,1]!
618 : //
619 : Array<Double> phaseDir = vb.msColumns().field().phaseDir().getColumn();
620 : Double RA0 = phaseDir(IPosition(3,0,0,0));
621 : Double Dec0 = phaseDir(IPosition(3,1,0,0));
622 : //
623 : // Compute reference (HA,Dec)
624 : //
625 : Double LST = LAST.get(Day).getValue();
626 : Double SDec0 = sin(Dec0), CDec0=cos(Dec0);
627 : LST -= floor(LST); // Extract the fractional day
628 : LST *= 2*C::pi;// Convert to Raidan
629 :
630 : Double HA0;
631 : HA0 = LST - RA0;
632 : Quantity QHA0(HA0,Radian), QDEC0(Dec0,Radian);
633 : //
634 : // Convert reference (HA,Dec) to reference (Az,El)
635 : //
636 : MDirection PhaseCenter(QHA0, QDEC0,MDirection::Ref(MDirection::HADEC));
637 : MDirection AzEl0 = toAzEl(PhaseCenter);
638 :
639 : MDirection tmpHADec = toHADec(AzEl0);
640 :
641 : Double Az0_Rad = AzEl0.getAngle(Radian).getValue()(0);
642 : Double El0_Rad = AzEl0.getAngle(Radian).getValue()(1);
643 :
644 : //
645 : // Convert the antenna pointing offsets from (Az,El)-->(RA,Dec)-->(l,m)
646 : //
647 :
648 : for(IPosition n(3,0,0,0);n(2)<=NAnt;n(2)++)
649 : {
650 : //
651 : // From (Az,El) -> (HA,Dec)
652 : //
653 : // Add (Az,El) offsets to the reference (Az,El)
654 : //
655 : dAz.setValue(l_off(n)+Az0_Rad); dEl.setValue(m_off(n)+El0_Rad);
656 : // dAz.setValue(0.0+Az0_Rad); dEl.setValue(0.0+El0_Rad);
657 : MDirection AzEl(dAz,dEl,MDirection::Ref(MDirection::AZEL));
658 : //
659 : // Convert offsetted (Az,El) to (HA,Dec) and then to (RA,Dec)
660 : //
661 : MDirection HADec = toHADec(AzEl);
662 : Double HA,Dec,RA, dRA;
663 : HA = HADec.getAngle(Radian).getValue()(0);
664 : Dec = HADec.getAngle(Radian).getValue()(1);
665 : RA = LST - HA;
666 : dRA = RA - RA0;
667 : //
668 : // Convert offsetted (RA,Dec) -> (l,m)
669 : //
670 : l_off(n) = sin(dRA)*cos(Dec);
671 : m_off(n) = sin(Dec)*CDec0-cos(Dec)*SDec0*cos(dRA);
672 : // cout << "FindOff: " << n(2) << " " << l_offsets(n) << " " << m_offsets(n) << endl;
673 :
674 : // cout << l_off(n) << " " << m_off(n) << " "
675 : // << " " << HA << " " << Dec
676 : // << " " << LST << " " << RA0 << " " << Dec0
677 : // << " " << RA << " " << Dec
678 : // << endl;
679 : }
680 :
681 : return NAnt+1;
682 0 : }
683 : //
684 : //---------------------------------------------------------------
685 : //
686 0 : int AWProjectFT::findPointingOffsets(const VisBuffer& vb,
687 : Cube<Float>& pointingOffsets,
688 : Array<Float> &l_off,
689 : Array<Float> &m_off,
690 : Bool Evaluate)
691 : {
692 0 : LogIO log_l(LogOrigin("AWProjectFT", "findPointingOffsets[R&D]"));
693 0 : Int NAnt = 0;
694 : //Float tmp;
695 : // TBD: adapt the following to VisCal mechanism:
696 0 : MEpoch LAST;
697 :
698 0 : NAnt=pointingOffsets.shape()(2);
699 0 : l_off.resize(IPosition(3,2,1,NAnt));
700 0 : m_off.resize(IPosition(3,2,1,NAnt));
701 0 : IPosition ndx(3,0,0,0),ndx1(3,0,0,0);
702 0 : for(ndx(2)=0;ndx(2)<NAnt;ndx(2)++)
703 : {
704 0 : ndx1=ndx;
705 0 : ndx(0)=0;ndx1(0)=0; //tmp=l_off(ndx) = pointingOffsets(ndx1);//Axis_0,Pol_0,Ant_i
706 0 : ndx(0)=1;ndx1(0)=1; //tmp=l_off(ndx) = pointingOffsets(ndx1);//Axis_0,Pol_1,Ant_i
707 0 : ndx(0)=0;ndx1(0)=2; //tmp=m_off(ndx) = pointingOffsets(ndx1);//Axis_1,Pol_0,Ant_i
708 0 : ndx(0)=1;ndx1(0)=3; //tmp=m_off(ndx) = pointingOffsets(ndx1);//Axis_1,Pol_1,Ant_i
709 : }
710 :
711 : // l_off = pointingOffsets(IPosition(3,0,0,0),IPosition(3,0,0,NAnt));
712 : // m_off = pointingOffsets(IPosition(3,1,0,0),IPosition(3,1,0,NAnt));
713 : /*
714 : IPosition shp(pointingOffsets.shape());
715 : IPosition shp1(l_off.shape()),shp2(m_off.shape());
716 : for(Int ii=0;ii<NAnt;ii++)
717 : {
718 : IPosition ndx(3,0,0,0);
719 : ndx(2)=ii;
720 : cout << "Pointing Offsets: " << ii << " "
721 : << l_off(ndx)*57.295*60.0 << " "
722 : << m_off(ndx)*57.295*60.0 << endl;
723 : }
724 : */
725 0 : return NAnt;
726 : if (!Evaluate) return NAnt;
727 :
728 : //
729 : // Make a Coordinate Conversion Machine to go from (Az,El) to
730 : // (HA,Dec).
731 : //
732 : MDirection::Convert toAzEl = makeCoordinateMachine(vb,MDirection::HADEC,
733 : MDirection::AZEL,
734 : LAST);
735 : MDirection::Convert toHADec = makeCoordinateMachine(vb,MDirection::AZEL,
736 : MDirection::HADEC,
737 : LAST);
738 : //
739 : // ...and now hope that it all works and works correctly!!!
740 : //
741 : Quantity dAz(0,Radian),dEl(0,Radian);
742 : //
743 : // An array of shape [2,1,1]!
744 : //
745 : Array<Double> phaseDir = vb.msColumns().field().phaseDir().getColumn();
746 : Double RA0 = phaseDir(IPosition(3,0,0,0));
747 : Double Dec0 = phaseDir(IPosition(3,1,0,0));
748 : //
749 : // Compute reference (HA,Dec)
750 : //
751 : Double LST = LAST.get(Day).getValue();
752 : Double SDec0 = sin(Dec0), CDec0=cos(Dec0);
753 : LST -= floor(LST); // Extract the fractional day
754 : LST *= 2*C::pi;// Convert to Raidan
755 :
756 : Double HA0;
757 : HA0 = LST - RA0;
758 : Quantity QHA0(HA0,Radian), QDEC0(Dec0,Radian);
759 : //
760 : // Convert reference (HA,Dec) to reference (Az,El)
761 : //
762 : MDirection PhaseCenter(QHA0, QDEC0,MDirection::Ref(MDirection::HADEC));
763 : MDirection AzEl0 = toAzEl(PhaseCenter);
764 :
765 : MDirection tmpHADec = toHADec(AzEl0);
766 :
767 : Double Az0_Rad = AzEl0.getAngle(Radian).getValue()(0);
768 : Double El0_Rad = AzEl0.getAngle(Radian).getValue()(1);
769 :
770 : //
771 : // Convert the antenna pointing offsets from (Az,El)-->(RA,Dec)-->(l,m)
772 : //
773 :
774 : for(IPosition n(3,0,0,0);n(2)<=NAnt;n(2)++)
775 : {
776 : //
777 : // From (Az,El) -> (HA,Dec)
778 : //
779 : // Add (Az,El) offsets to the reference (Az,El)
780 : //
781 : dAz.setValue(l_off(n)+Az0_Rad); dEl.setValue(m_off(n)+El0_Rad);
782 : // dAz.setValue(0.0+Az0_Rad); dEl.setValue(0.0+El0_Rad);
783 : MDirection AzEl(dAz,dEl,MDirection::Ref(MDirection::AZEL));
784 : //
785 : // Convert offsetted (Az,El) to (HA,Dec) and then to (RA,Dec)
786 : //
787 : MDirection HADec = toHADec(AzEl);
788 : Double HA,Dec,RA, dRA;
789 : HA = HADec.getAngle(Radian).getValue()(0);
790 : Dec = HADec.getAngle(Radian).getValue()(1);
791 : RA = LST - HA;
792 : dRA = RA - RA0;
793 : //
794 : // Convert offsetted (RA,Dec) -> (l,m)
795 : //
796 : l_off(n) = sin(dRA)*cos(Dec);
797 : m_off(n) = sin(Dec)*CDec0-cos(Dec)*SDec0*cos(dRA);
798 : // cout << "FindOff: " << n(2) << " " << l_offsets(n) << " " << m_offsets(n) << endl;
799 :
800 : // cout << l_off(n) << " " << m_off(n) << " "
801 : // << " " << HA << " " << Dec
802 : // << " " << LST << " " << RA0 << " " << Dec0
803 : // << " " << RA << " " << Dec
804 : // << endl;
805 : }
806 :
807 : return NAnt+1;
808 0 : }
809 : //
810 : //---------------------------------------------------------------
811 : //
812 0 : void AWProjectFT::makeSensitivityImage(const VisBuffer& vb,
813 : const ImageInterface<Complex>& imageTemplate,
814 : ImageInterface<Float>& sensitivityImage)
815 : {
816 0 : if (convFuncCtor_p->makeAverageResponse(vb, imageTemplate, sensitivityImage))
817 0 : cfCache_p->flush(sensitivityImage,sensitivityPatternQualifierStr_p);
818 0 : }
819 : //
820 : //---------------------------------------------------------------
821 : //
822 0 : void AWProjectFT::normalizeAvgPB(ImageInterface<Complex>& inImage,
823 : ImageInterface<Float>& outImage)
824 : {
825 0 : LogIO log_l(LogOrigin("AWProjectFT", "normalizeAvgPB[R&D]"));
826 0 : if (pbNormalized_p) return;
827 0 : IPosition inShape(inImage.shape()),ndx(4,0,0,0,0);
828 0 : Vector<Complex> peak(inShape(2));
829 :
830 0 : outImage.resize(inShape);
831 0 : outImage.setCoordinateInfo(inImage.coordinates());
832 :
833 : Bool isRefIn;
834 0 : Array<Complex> inBuf;
835 0 : Array<Float> outBuf;
836 :
837 0 : isRefIn = inImage.get(inBuf);
838 : //isRefOut = outImage.get(outBuf);
839 : log_l << "Normalizing the average PBs to unity"
840 0 : << LogIO::NORMAL << LogIO::POST;
841 : //
842 : // Normalize each plane of the inImage separately to unity.
843 : //
844 0 : Complex inMax = max(inBuf);
845 0 : if (abs(inMax)-1.0 > 1E-3)
846 : {
847 0 : for(ndx(3)=0;ndx(3)<inShape(3);ndx(3)++)
848 0 : for(ndx(2)=0;ndx(2)<inShape(2);ndx(2)++)
849 : {
850 0 : peak(ndx(2)) = 0;
851 0 : for(ndx(1)=0;ndx(1)<inShape(1);ndx(1)++)
852 0 : for(ndx(0)=0;ndx(0)<inShape(0);ndx(0)++)
853 0 : if (abs(inBuf(ndx)) > peak(ndx(2)))
854 0 : peak(ndx(2)) = inBuf(ndx);
855 :
856 0 : for(ndx(1)=0;ndx(1)<inShape(1);ndx(1)++)
857 0 : for(ndx(0)=0;ndx(0)<inShape(0);ndx(0)++)
858 : // avgPBBuf(ndx) *= (pbPeaks(ndx(2))/peak(ndx(2)));
859 0 : inBuf(ndx) /= peak(ndx(2));
860 : }
861 0 : if (isRefIn) inImage.put(inBuf);
862 : }
863 :
864 0 : ndx=0;
865 0 : for(ndx(1)=0;ndx(1)<inShape(1);ndx(1)++)
866 0 : for(ndx(0)=0;ndx(0)<inShape(0);ndx(0)++)
867 : {
868 0 : IPosition plane1(ndx);
869 0 : plane1=ndx;
870 0 : plane1(2)=1; // The other poln. plane
871 : // avgPBBuf(ndx) = (avgPBBuf(ndx) + avgPBBuf(plane1))/2.0;
872 0 : outBuf(ndx) = sqrt(real(inBuf(ndx) * inBuf(plane1)));
873 0 : }
874 : //
875 : // Rather convoluted way of copying Pol. plane-0 to Pol. plane-1!!!
876 : //
877 0 : for(ndx(1)=0;ndx(1)<inShape(1);ndx(1)++)
878 0 : for(ndx(0)=0;ndx(0)<inShape(0);ndx(0)++)
879 : {
880 0 : IPosition plane1(ndx);
881 0 : plane1=ndx;
882 0 : plane1(2)=1; // The other poln. plane
883 0 : outBuf(plane1) = real(outBuf(ndx));
884 0 : }
885 :
886 0 : pbNormalized_p = true;
887 0 : }
888 : //
889 : //---------------------------------------------------------------
890 : //
891 0 : void AWProjectFT::normalizeAvgPB()
892 : {
893 0 : LogIO log_l(LogOrigin("AWProjectFT", "normalizeAvgPB[R&D]"));
894 0 : if (pbNormalized_p) return;
895 : Bool isRefF;
896 0 : Array<Float> avgPBBuf;
897 0 : isRefF=avgPB_p->get(avgPBBuf);
898 : // Float pbMax = max(avgPBBuf);
899 : {
900 0 : pbPeaks.resize(avgPB_p->shape()(2),true);
901 : // if (makingPSF) pbPeaks = 1.0;
902 : // else pbPeaks /= (Float)noOfPASteps;
903 0 : pbPeaks = 1.0;
904 : log_l << "Normalizing the average PBs to " << 1.0
905 0 : << LogIO::NORMAL << LogIO::POST;
906 :
907 0 : IPosition avgPBShape(avgPB_p->shape()),ndx(4,0,0,0,0);
908 0 : Vector<Float> peak(avgPBShape(2));
909 :
910 :
911 0 : Float pbMax = max(avgPBBuf);
912 0 : if (fabs(pbMax-1.0) > 1E-3)
913 : {
914 : // avgPBBuf = avgPBBuf/noOfPASteps;
915 0 : for(ndx(3)=0;ndx(3)<avgPBShape(3);ndx(3)++)
916 0 : for(ndx(2)=0;ndx(2)<avgPBShape(2);ndx(2)++)
917 : {
918 0 : peak(ndx(2)) = 0;
919 0 : for(ndx(1)=0;ndx(1)<avgPBShape(1);ndx(1)++)
920 0 : for(ndx(0)=0;ndx(0)<avgPBShape(0);ndx(0)++)
921 0 : if (abs(avgPBBuf(ndx)) > peak(ndx(2)))
922 0 : peak(ndx(2)) = avgPBBuf(ndx);
923 :
924 0 : for(ndx(1)=0;ndx(1)<avgPBShape(1);ndx(1)++)
925 0 : for(ndx(0)=0;ndx(0)<avgPBShape(0);ndx(0)++)
926 : // avgPBBuf(ndx) *= (pbPeaks(ndx(2))/peak(ndx(2)));
927 0 : avgPBBuf(ndx) /= peak(ndx(2));
928 : }
929 0 : if (isRefF) avgPB_p->put(avgPBBuf);
930 : }
931 :
932 0 : ndx=0;
933 0 : for(ndx(1)=0;ndx(1)<avgPBShape(1);ndx(1)++)
934 0 : for(ndx(0)=0;ndx(0)<avgPBShape(0);ndx(0)++)
935 : {
936 0 : IPosition plane1(ndx);
937 0 : plane1=ndx;
938 0 : plane1(2)=1; // The other poln. plane
939 0 : avgPBBuf(ndx) = (avgPBBuf(ndx) + avgPBBuf(plane1))/2.0;
940 : // avgPBBuf(ndx) = (avgPBBuf(ndx) * avgPBBuf(plane1));
941 0 : }
942 0 : for(ndx(1)=0;ndx(1)<avgPBShape(1);ndx(1)++)
943 0 : for(ndx(0)=0;ndx(0)<avgPBShape(0);ndx(0)++)
944 : {
945 0 : IPosition plane1(ndx);
946 0 : plane1=ndx;
947 0 : plane1(2)=1; // The other poln. plane
948 0 : avgPBBuf(plane1) = avgPBBuf(ndx);
949 0 : }
950 0 : }
951 0 : pbNormalized_p = true;
952 0 : }
953 : //
954 : //---------------------------------------------------------------
955 0 : void AWProjectFT::makeCFPolMap(const VisBuffer& vb, const Vector<Int>& locCfStokes,
956 : Vector<Int>& polM)
957 : {
958 0 : LogIO log_l(LogOrigin("AWProjectFT", "makeCFPolMap[R&D]"));
959 0 : Vector<Int> msStokes = vb.corrType();
960 0 : Int nPol = msStokes.nelements();
961 0 : polM.resize(polMap.shape());
962 0 : polM = -1;
963 :
964 0 : for(Int i=0;i<nPol;i++)
965 0 : for(uInt j=0;j<locCfStokes.nelements();j++)
966 0 : if (locCfStokes(j) == msStokes(i))
967 0 : {polM(i) = j;break;}
968 0 : }
969 : //
970 : //---------------------------------------------------------------
971 : //
972 : // Given a polMap (mapping of which Visibility polarization is
973 : // gridded onto which grid plane), make a map of the conjugate
974 : // planes of the grid E.g, for Stokes-I and -V imaging, the two
975 : // planes of the uv-grid are [LL,RR]. For input VisBuffer
976 : // visibilites in order [RR,RL,LR,LL], polMap = [1,-1,-1,0]. The
977 : // conjugate map will be [0,-1,-1,1].
978 : //
979 0 : void AWProjectFT::makeConjPolMap(const VisBuffer& vb,
980 : const Vector<Int> cfPolMap,
981 : Vector<Int>& conjPolMap)
982 : {
983 0 : if (conjPolMap.nelements() > 0) return;
984 :
985 0 : LogIO log_l(LogOrigin("AWProjectFT", "makConjPolMap[R&D]"));
986 :
987 : //
988 : // All the Natak (Drama) below with slicers etc. is to extract the
989 : // Poln. info. for the first IF only (not much "information
990 : // hiding" for the code to slice arrays in a general fashion).
991 : //
992 : // Extract the shape of the array to be sliced.
993 : //
994 : log_l << "############....temp code!!!!!!!!!! "
995 0 : << SynthesisUtils::mapSpwIDToPolID(vb,selectedSpw_p(0))
996 0 : << LogIO::DEBUG1;
997 0 : Vector<Int> polIDs=SynthesisUtils::mapSpwIDToPolID(vb,selectedSpw_p(0));
998 0 : if (polIDs.nelements()==0)
999 0 : log_l << "Internal Error: Selected SPW did not map to any pol!" << LogIO::EXCEPTION;
1000 : //
1001 : // Use the first selected Spw to determine the pol. mapping in the
1002 : // VB. This implies that the code assumes that all selected SPWs
1003 : // have the same pol. setup.
1004 : //
1005 0 : Int polID=polIDs(0);
1006 0 : Array<Int> stokesForAllIFs = vb.msColumns().polarization().corrType().get(polID);
1007 0 : IPosition stokesShape(stokesForAllIFs.shape());
1008 0 : IPosition firstIFStart(stokesShape),firstIFLength(stokesShape);
1009 : //
1010 : // Set up the start and length IPositions to extract only the
1011 : // first column of the array. The following is required since the
1012 : // array could have only one column as well.
1013 : //
1014 0 : firstIFStart(0)=0;firstIFLength(0)=stokesShape(0);
1015 0 : for(uInt i=1;i<stokesShape.nelements();i++) {firstIFStart(i)=0;firstIFLength(i)=1;}
1016 : //
1017 : // Construct the slicer and produce the slice. .nonDegenerate
1018 : // required to ensure the result of slice is a pure vector.
1019 : //
1020 0 : Vector<Int> visStokes = stokesForAllIFs(Slicer(firstIFStart,firstIFLength)).nonDegenerate();
1021 :
1022 0 : conjPolMap = cfPolMap;
1023 :
1024 0 : Int i,j,N = cfPolMap.nelements();
1025 0 : for(i=0;i<N;i++)
1026 0 : if (cfPolMap[i] > -1)
1027 : {
1028 0 : if (visStokes[i] == Stokes::RR)
1029 : {
1030 0 : conjPolMap[i]=-1;
1031 0 : for(j=0;j<N;j++) if (visStokes[j] == Stokes::LL) break;
1032 0 : conjPolMap[i]=cfPolMap[j];
1033 : }
1034 0 : else if (visStokes[i] == Stokes::LL)
1035 : {
1036 0 : conjPolMap[i]=-1;
1037 0 : for(j=0;j<N;j++) if (visStokes[j] == Stokes::RR) break;
1038 0 : conjPolMap[i]=cfPolMap[j];
1039 : }
1040 0 : else if (visStokes[i] == Stokes::LR)
1041 : {
1042 0 : conjPolMap[i]=-1;
1043 0 : for(j=0;j<N;j++) if (visStokes[j] == Stokes::RL) break;
1044 0 : conjPolMap[i]=cfPolMap[j];
1045 : }
1046 0 : else if (visStokes[i] == Stokes::RL)
1047 : {
1048 0 : conjPolMap[i]=-1;
1049 0 : for(j=0;j<N;j++) if (visStokes[j] == Stokes::LR) break;
1050 0 : conjPolMap[i]=cfPolMap[j];
1051 : }
1052 : }
1053 0 : }
1054 : //
1055 : //---------------------------------------------------------------
1056 : // Convert the CFs in the supplied CFStore to
1057 : // A(nu)<Convolution>A(nu_*)
1058 : //
1059 0 : void AWProjectFT::makeWBCFWt(CFStore2& cfs, const Double imRefFreq)
1060 : {
1061 0 : LogIO log_l(LogOrigin("AWProjectFT", "makeWBCFWt[R&D]"));
1062 0 : log_l << "Converting WTCFs to wide-band versions" << LogIO::POST;
1063 :
1064 0 : Vector<Int> ant1List, ant2List;
1065 0 : Vector<Quantity> paList;
1066 0 : ant1List = cfs.getAnt1List();
1067 0 : ant2List = cfs.getAnt2List();
1068 0 : paList = cfs.getPAList();
1069 :
1070 0 : if (paNdxProcessed_p.nelements() == 0) {paNdxProcessed_p.resize(1); paNdxProcessed_p[0]=false;}
1071 0 : CountedPtr<CFBuffer> cfb_l, cfb_clone;
1072 0 : Quantity dPA(360.0,"deg");
1073 0 : for (uInt pa=0;pa<paList.nelements();pa++)
1074 0 : for (uint a1=0;a1<ant1List.nelements(); a1++)
1075 0 : for (uint a2=0;a2<ant2List.nelements(); a2++)
1076 : {
1077 0 : if (paNdxProcessed_p.nelements() < pa) {paNdxProcessed_p.resize(pa+1,true); paNdxProcessed_p[pa]=false;}
1078 0 : if (!paNdxProcessed_p[pa])
1079 : {
1080 0 : paNdxProcessed_p[pa]=true;
1081 0 : Vector<Double> wVals, fVals; PolMapType mVals, mNdx, conjMVals, conjMNdx;
1082 : Double fIncr, wIncr;
1083 0 : cfb_l = cfs.getCFBuffer(paList[pa], dPA, ant1List(a1), ant2List(a2));
1084 0 : cfb_l->getCoordList(fVals,wVals,mNdx, mVals, conjMNdx, conjMVals, fIncr, wIncr);
1085 :
1086 : // cerr << paList << " " << endl
1087 : // << wVals << " " << endl
1088 : // << fVals << " " << endl
1089 : // << conjMVals << " " << endl;
1090 :
1091 0 : cfb_clone=cfb_l->clone();
1092 0 : cfb_clone->getCoordList(fVals,wVals,mNdx, mVals, conjMNdx, conjMVals, fIncr, wIncr);
1093 : // cerr << "Clone: "
1094 : // << paList << " " << endl
1095 : // << wVals << " " << endl
1096 : // << fVals << " " << endl
1097 : // << conjMVals << " " << endl;
1098 :
1099 0 : CountedPtr<Array<Complex> > cfc_l, conjCFC_l, result_l;
1100 : //
1101 : // Damn! Convolver does not work for complex convolution!
1102 : // Convolver<Complex> convolver;
1103 :
1104 0 : for (Int nw=0;nw<(Int)wVals.nelements(); nw++)
1105 0 : for (Int nf=0;nf<(Int)fVals.nelements(); nf++)
1106 0 : for (Int ipol=0;ipol<(Int)conjMVals.nelements();ipol++)
1107 0 : for (Int mRow=0;mRow<(Int)conjMVals[ipol].nelements(); mRow++)
1108 : {
1109 : Double f, cf;
1110 0 : f=fVals[nf];
1111 0 : cf=sqrt(2*imRefFreq*imRefFreq - f*f);
1112 0 : Int conjNF = cfb_l->nearestFreqNdx(cf);
1113 : // cerr << "###PTRs: "
1114 : // << &(*(cfb_l->getCFCellPtr(nf,nw,conjMNdx[ipol][mRow]))) << " "
1115 : // << &(*(cfb_clone->getCFCellPtr(nf,nw,conjMNdx[ipol][mRow]))) << " "
1116 : // << &(*(cfb_clone->getCFCellPtr(conjNF,nw,conjMNdx[ipol][mRow])))
1117 : // << endl;
1118 0 : result_l = cfb_l->getCFCellPtr(nf,nw,conjMNdx[ipol][mRow])->getStorage();
1119 0 : cfc_l = cfb_clone->getCFCellPtr(nf,nw,conjMNdx[ipol][mRow])->getStorage();
1120 0 : conjCFC_l = cfb_clone->getCFCellPtr(conjNF,nw,conjMNdx[ipol][mRow])->getStorage();
1121 :
1122 : // convolver.setPsf(*cfc_l);
1123 : // convolver.linearConv(*result_l, *conjCFC_l);
1124 0 : SynthesisUtils::libreConvolver(*result_l,*conjCFC_l);
1125 :
1126 :
1127 : // {
1128 : // CountedPtr<CFCell> cfc_t, conjCFC_t;
1129 : // cfc_t = cfb_l->getCFCellPtr(nf,nw,conjMNdx[ipol][mRow]);
1130 : // conjCFC_t = cfb_l->getCFCellPtr(conjNF,nw,conjMNdx[ipol][mRow]);
1131 : // cout << "-------------------------------" << " " << f << " " << cf << " "
1132 : // << nf << " " << nw << " " << conjMNdx[ipol][mRow] << " "
1133 : // << conjNF << " " << nw << " " << conjMNdx[ipol][mRow] << " "
1134 : // <<endl;
1135 : // cfc_t->show("CFC: ",cout);
1136 : // conjCFC_t->show("conjCFC: ",cout);
1137 : // cout << "-------------------------------" << endl;
1138 : // }
1139 : }
1140 0 : }
1141 : }
1142 0 : }
1143 : //
1144 : // Locate a convlution function. It will be either in the cache
1145 : // (mem. or disk cache) or will be computed and cached for possible
1146 : // later use.
1147 : //
1148 0 : void AWProjectFT::findConvFunction(const ImageInterface<Complex>& image,
1149 : const VisBuffer& vb)
1150 : {
1151 0 : LogIO log_l(LogOrigin("AWProjectFT", "findConvFunction[R&D]"));
1152 0 : if (!paChangeDetector.changed(vb,0)) return;
1153 0 : Int cfSource=CFDefs::NOTCACHED;
1154 0 : CoordinateSystem ftcoords;
1155 : // Think of a generic call to get the key-values. And a
1156 : // overloadable method (or an externally supplied one?) to convert
1157 : // the values to key-ids. That will ensure that AWProjectFT
1158 : // remains the A-Projection algorithm implementation configurable
1159 : // by the behaviour of the supplied objects.
1160 0 : Float pa=getVBPA(vb);
1161 : //UUU// ok();
1162 0 : visResampler_p->setMaps(chanMap, polMap); //UUU Added here.
1163 0 : visResampler_p->setFreqMaps(expandedSpwFreqSel_p,expandedSpwConjFreqSel_p);
1164 :
1165 0 : lastPAUsedForWtImg = currentCFPA = pa;
1166 :
1167 0 : Vector<Double> pointingOffset(convFuncCtor_p->findPointingOffset(image,vb));
1168 0 : Float dPA = paChangeDetector.getParAngleTolerance().getValue("rad");
1169 0 : Quantity dPAQuant = Quantity(paChangeDetector.getParAngleTolerance());
1170 0 : cfSource = visResampler_p->makeVBRow2CFMap(*cfs2_p,*convFuncCtor_p, vb,
1171 : dPAQuant,
1172 0 : chanMap,polMap,pointingOffset);
1173 :
1174 0 : if (cfSource == CFDefs::NOTCACHED)
1175 : {
1176 0 : PolMapType polMat, polIndexMat, conjPolMat, conjPolIndexMat;
1177 0 : Vector<Int> visPolMap(vb.corrType());
1178 0 : polMat = pop_p->makePolMat(visPolMap,polMap);
1179 0 : polIndexMat = pop_p->makePol2CFMat(visPolMap,polMap);
1180 :
1181 : // polMat = pop_p->makePolMat(visPolMap,polMap);
1182 : // polIndexMat = pop_p->makePol2CFMat(visPolMap,polMap);
1183 :
1184 0 : conjPolMat = pop_p->makeConjPolMat(visPolMap,polMap);
1185 0 : conjPolIndexMat = pop_p->makeConjPol2CFMat(visPolMap,polMap);
1186 :
1187 0 : convFuncCtor_p->setPolMap(polMap);
1188 0 : convFuncCtor_p->setSpwSelection(spwChanSelFlag_p);
1189 0 : convFuncCtor_p->setSpwFreqSelection(spwFreqSel_p);
1190 :
1191 : // USEFUL DEBUG MESSAGE
1192 : //cerr << "Freq. selection: " << expandedSpwFreqSel_p << endl << expandedSpwConjFreqSel_p << endl;
1193 :
1194 :
1195 : // isDryRun = true;
1196 0 : Bool pleaseDoAlsoFillTheCF=!dryRun();
1197 0 : convFuncCtor_p->makeConvFunction(image,vb,wConvSize,
1198 0 : pop_p, pa, dPA, uvScale, uvOffset,spwFreqSel_p,
1199 : *cfs2_p, *cfwts2_p, pleaseDoAlsoFillTheCF);
1200 : // log_l << "Converting WTCFs to wide-band versions" << LogIO::POST;
1201 : // makeWBCFWt(*cfwts2_p, imRefFreq_p);
1202 :
1203 : // cfSource = visResampler_p->makeVBRow2CFMap(*cfs2_p,*convFuncCtor_p, vb,
1204 : // paChangeDetector.getParAngleTolerance(),
1205 : // chanMap,polMap);
1206 0 : }
1207 : //
1208 : // If one-time-operations in the CFCache not yet done, set the
1209 : // pol. index maps in the CFCache.
1210 : //
1211 0 : if (!cfCache_p->OTODone())
1212 : {
1213 0 : Vector<Int> visPolMap(vb.corrType());
1214 :
1215 0 : PolMapType polMat, conjPolMat;
1216 0 : polMat = pop_p->makePolMat(visPolMap,polMap);
1217 0 : conjPolMat = pop_p->makeConjPolMat(visPolMap,polMap);
1218 :
1219 0 : PolMapType pNdx, cpNdx;
1220 0 : pNdx = pop_p->makePol2CFMat(visPolMap,polMap);
1221 0 : cpNdx = pop_p->makeConjPol2CFMat(visPolMap,polMap);
1222 :
1223 0 : cfCache_p->initPolMaps(pNdx,cpNdx);
1224 :
1225 0 : cfs2_p->initMaps(vb,spwFreqSel_p,imRefFreq_p);
1226 0 : cfwts2_p->initMaps(vb,spwFreqSel_p,imRefFreq_p);
1227 0 : }
1228 : //
1229 : // Load the average PB (sensitivity pattern) from the cache. If
1230 : // not found in the cache, make one and cache it.
1231 : //
1232 0 : if (cfCache_p->loadAvgPB(avgPB_p,sensitivityPatternQualifierStr_p) == CFDefs::NOTCACHED)
1233 0 : makeSensitivityImage(vb,image,*avgPB_p);
1234 :
1235 0 : verifyShapes(avgPB_p->shape(), image.shape());
1236 :
1237 0 : if (paChangeDetector.changed(vb,0)) paChangeDetector.update(vb,0);
1238 : //
1239 : // Write some useful info. to the logger.
1240 : //
1241 0 : if (cfSource != CFDefs::MEMCACHE)
1242 : {
1243 0 : if (dryRun())
1244 : {
1245 : // PagedImage<Complex> thisGrid(lattice->shape(),image.coordinates(),
1246 : // cfCache_p->getCacheDir()+"/uvgrid.im");
1247 0 : PagedImage<Complex> thisGrid(image.shape(),image.coordinates(),
1248 0 : cfCache_p->getCacheDir()+"/uvgrid.im");
1249 0 : }
1250 :
1251 : // cfs2_p->makePersistent(cfCache_p->getCacheDir().c_str());
1252 : // cfwts2_p->makePersistent(cfCache_p->getCacheDir().c_str(),"","WT");
1253 :
1254 : // Save only the CF Cube for the current value of PA (not the
1255 : // entire CFStore -- CFs for PA values encountered earlier
1256 : // than current value have already need made persistent).
1257 0 : cfs2_p->makePersistent(cfCache_p->getCacheDir().c_str(),"","", Quantity(pa,"rad"),dPAQuant,0,0);
1258 0 : cfwts2_p->makePersistent(cfCache_p->getCacheDir().c_str(),"","WT",Quantity(pa,"rad"),dPAQuant,0,0);
1259 0 : Double memUsed=cfs2_p->memUsage();
1260 0 : String unit(" KB");
1261 0 : memUsed = (Int)(memUsed/1024.0+0.5);
1262 0 : if (memUsed > 1024) {memUsed /=1024; unit=" MB";}
1263 : log_l << "Convolution function memory footprint:"
1264 : << (Int)(memUsed) << unit << " out of a maximum of "
1265 0 : << HostInfo::memoryTotal(true)/1024 << " MB" << LogIO::POST;
1266 :
1267 : //
1268 : // Initialize any internal maps that may be used later for
1269 : // efficient access.
1270 : //
1271 0 : cfs2_p->initMaps(vb,spwFreqSel_p,imRefFreq_p);
1272 0 : cfwts2_p->initMaps(vb,spwFreqSel_p,imRefFreq_p);
1273 0 : }
1274 :
1275 : // cfs2_p->makePersistent("test.cf");
1276 : // cfwts2_p->makePersistent("test.wtcf");
1277 0 : }
1278 : //
1279 : //------------------------------------------------------------------------------
1280 : // Vectorized initializeToVis. See design related comments in the
1281 : // .h file.
1282 : //
1283 : //
1284 : // The functional goals here are:
1285 : //
1286 : // 1. If doPBCorrection==true (i.e. the input sky-image is flat-sky
1287 : // image), divide the image by the PB
1288 : //
1289 : // 2. Convert from Stokes to Feed (correlation) frame
1290 : //
1291 : // 3. FFT the image to produce gridded vis.
1292 : //
1293 : // 4. And since the same image buffer is used to accumulate the
1294 : // model, multiply the sky-image back with PB
1295 : //
1296 : // The call to non-vectorized version of initializeToVis() which is
1297 : // pure virtual and hence the local version is called, is only to do
1298 : // the FFT. FFT is *NOT* in place.
1299 : //
1300 : // These operations effecitvely maintains the model image as PB*Sky
1301 : // and supplies flat-sky model for prediction. This can probably be
1302 : // achieve with fewer operations and same memory buffers....but that
1303 : // for later (SB)
1304 0 : void AWProjectFT::initializeToVis(Block<CountedPtr<ImageInterface<Complex> > > & compImageVec,
1305 : PtrBlock<SubImage<Float> *> & modelImageVec,
1306 : PtrBlock<SubImage<Float> *>& weightImageVec,
1307 : PtrBlock<SubImage<Float> *>& /*fluxScaleVec*/,
1308 : Block<Matrix<Float> >& weightsVec,
1309 : const VisBuffer& vb)
1310 : {
1311 0 : LogIO log_p(LogOrigin("AWProjectFT","initToVis[V][R&D]"));
1312 : //
1313 : // Setting the image below is crucial since init() and
1314 : // initMaps(vb) below expect this to be set.
1315 : //
1316 0 : image=&(*compImageVec[0]);
1317 :
1318 : log_p << "Total flux in model image (before avgPB normalization): "
1319 0 : << sum((*(modelImageVec[0])).get())
1320 : // << " predicing SPW = " <<vb.spectralWindow()
1321 : // << " Pointing Offset = " << convFuncCtor_p->findPointingOffset(*(compImageVec[0]), vb)
1322 : // << " Qualifier String = " << sensitivityPatternQualifier_p
1323 0 : << LogIO::POST;
1324 0 : if(doPBCorrection)
1325 : {
1326 : // Make the sensitivity Image if applicable
1327 0 : init();
1328 0 : initMaps(vb);
1329 0 : findConvFunction(*(compImageVec[0]), vb); // Pure virtual -- call local version
1330 :
1331 0 : if (isDryRun) return;
1332 :
1333 : // Get the sensitivity Image
1334 0 : Matrix<Float> tempWts;
1335 0 : tempWts.resize();
1336 0 : getWeightImage(*(weightImageVec[0]), tempWts); // Pure virtual -- call local version
1337 :
1338 : // Normalize the model image by the sensitivity image only.
1339 : // No local implementation -- call FTMachine version
1340 :
1341 : // Divide by avgPB ///// PBWeight
1342 : //
1343 : //normalizeImage( *(modelImageVec[0]) , weightsVec[0], *(weightImageVec[0]),
1344 : // false, (Float)pbLimit_p, (Int)1);
1345 :
1346 : //
1347 : // Divide by sqrt(avgPB) ////// PBSQWeight
1348 : //
1349 0 : normalizeImage( *(modelImageVec[0]) , weightsVec[0], *(weightImageVec[0]),
1350 0 : false, (Float)pbLimit_p, (Int)4);
1351 0 : }
1352 : log_p << "Total flux in model image (after avgPB normalization): "
1353 0 : << sum((*(modelImageVec[0])).get()) << LogIO::POST;
1354 :
1355 : // Convert from Stokes planes to Correlation planes
1356 : // No local implementation -- call FTMachine version
1357 0 : stokesToCorrelation(*(modelImageVec[0]), *(compImageVec[0]));
1358 :
1359 : // Call initializeToVis
1360 0 : initializeToVis(*(compImageVec[0]), vb); // Pure virtual
1361 :
1362 : // Multiply the flat-sky model by the PB.
1363 : // No local implementation -- call FTMachine version
1364 0 : if(doPBCorrection)
1365 : // Multiply by avgPB /////// PBWeight
1366 : //normalizeImage( *(modelImageVec[0]) , weightsVec[0], *(weightImageVec[0]) , false, (Float)pbLimit_p, (Int)3);
1367 :
1368 : // Multiply by sqrt(avgPB) ///// PBSQWeight
1369 0 : normalizeImage( *(modelImageVec[0]) , weightsVec[0], *(weightImageVec[0]) , false, (Float)pbLimit_p, (Int)5);
1370 0 : }
1371 : //
1372 : //------------------------------------------------------------------------------
1373 : //
1374 0 : void AWProjectFT::initializeToVis(ImageInterface<Complex>& iimage,
1375 : const VisBuffer& vb)
1376 : {
1377 0 : LogIO log_l(LogOrigin("AWProjectFT", "initializeToVis[R&D]"));
1378 0 : image=&iimage;
1379 :
1380 0 : ok();
1381 :
1382 0 : init();
1383 0 : makingPSF = false;
1384 0 : initMaps(vb);
1385 : //visResampler_p->setMaps(chanMap, polMap);
1386 :
1387 0 : findConvFunction(*image, vb);
1388 0 : if (isDryRun) return;
1389 : //UUU if (!cfCache_p->avgPBReady() && doPBCorrection)
1390 : //UUU log_l << "Sensitivity pattern not found." << LogIO::EXCEPTION;
1391 : //UUU// if (!cfCache_p->avgPBReady(sensitivityPatternQualifierStr_p) && doPBCorrection)
1392 : //UUU// log_l << "Sensitivity pattern for term " << sensitivityPatternQualifierStr_p << " not found." << LogIO::EXCEPTION;
1393 :
1394 : //
1395 : // Initialize the maps for polarization and channel. These maps
1396 : // translate visibility indices into image indices
1397 : //
1398 :
1399 0 : nx = image->shape()(0);
1400 0 : ny = image->shape()(1);
1401 0 : npol = image->shape()(2);
1402 0 : nchan = image->shape()(3);
1403 :
1404 0 : if(image->shape().product()>cachesize) isTiled=true;
1405 0 : else isTiled=false;
1406 : //
1407 : // If we are memory-based then read the image in and create an
1408 : // ArrayLattice otherwise just use the PagedImage
1409 : //
1410 :
1411 0 : isTiled=false;
1412 :
1413 0 : if(isTiled){
1414 0 : lattice=CountedPtr<Lattice<Complex> > (image, false);
1415 : }
1416 : else
1417 : {
1418 0 : IPosition gridShape(4, nx, ny, npol, nchan);
1419 0 : griddedData.resize(gridShape);
1420 0 : griddedData=Complex(0.0);
1421 :
1422 0 : IPosition stride(4, 1);
1423 0 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
1424 0 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
1425 0 : IPosition trc(blc+image->shape()-stride);
1426 :
1427 0 : IPosition start(4, 0);
1428 0 : griddedData(blc, trc) = image->getSlice(start, image->shape());
1429 :
1430 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
1431 0 : lattice=arrayLattice;
1432 0 : }
1433 :
1434 : //AlwaysAssert(lattice, AipsError);
1435 :
1436 0 : log_l << LogIO::DEBUGGING << "Starting FFT of image" << LogIO::POST;
1437 :
1438 0 : Vector<Float> sincConv(nx);
1439 0 : Float centerX=nx/2;
1440 0 : for (Int ix=0;ix<nx;ix++)
1441 : {
1442 0 : Float x=C::pi*Float(ix-centerX)/(Float(nx)*Float(convSampling));
1443 0 : if(ix==centerX) sincConv(ix)=1.0;
1444 0 : else sincConv(ix)=sin(x)/x;
1445 : }
1446 :
1447 : // Vector<Complex> correction(nx);
1448 : //
1449 : // Do the Grid-correction
1450 : //
1451 : // if(0) //UUU
1452 0 : if(cfCache_p->avgPBReady()) //SB
1453 : {
1454 : // normalizeAvgPB();
1455 :
1456 0 : IPosition cursorShape(4, nx, 1, 1, 1);
1457 0 : IPosition axisPath(4, 0, 1, 2, 3);
1458 0 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
1459 0 : LatticeIterator<Complex> lix(*lattice, lsx);
1460 :
1461 0 : verifyShapes(avgPB_p->shape(), image->shape());
1462 0 : Array<Float> avgBuf; avgPB_p->get(avgBuf);
1463 : // If the total-power sensitivity pattern peak is too low, warn
1464 : // the user. This usually is indicative of a rat somewhere in
1465 : // the pipes upstream...
1466 0 : if ((sensitivityPatternQualifier_p==0) && (max(avgBuf) < 1e-04))
1467 : log_l << "Normalization by PB requested but either PB was not"
1468 : <<" found in the cache or is ill-formed. Peak = "
1469 0 : << max(avgBuf)// << " " << sensitivityPatternQualifier_p
1470 0 : << LogIO::WARN << LogIO::POST;
1471 :
1472 :
1473 0 : LatticeStepper lpb(avgPB_p->shape(),cursorShape,axisPath);
1474 0 : LatticeIterator<Float> lipb(*avgPB_p, lpb);
1475 :
1476 0 : Vector<Complex> griddedVis;
1477 : //
1478 : // Grid correct in anticipation of the convolution by the
1479 : // convFunc. Each polarization plane is corrected by the
1480 : // appropraite primary beam.
1481 : //
1482 0 : for(lix.reset(),lipb.reset();!lix.atEnd();lix++,lipb++)
1483 : {
1484 0 : Int iy=lix.position()(1);
1485 : // gridder->correctX1D(correction,iy);
1486 0 : griddedVis = lix.rwVectorCursor();
1487 :
1488 0 : Vector<Float> PBCorrection(lipb.rwVectorCursor().shape());
1489 0 : PBCorrection = lipb.rwVectorCursor();
1490 0 : for(int ix=0;ix<nx;ix++)
1491 : {
1492 0 : if (doPBCorrection)
1493 : {
1494 0 : PBCorrection(ix) = pbFunc(PBCorrection(ix),pbLimit_p)*(sincConv(ix)*sincConv(iy));
1495 0 : lix.rwVectorCursor()(ix) /= (PBCorrection(ix));
1496 : }
1497 : else
1498 0 : lix.rwVectorCursor()(ix) /= (1.0/(sincConv(ix)*sincConv(iy)));
1499 : }
1500 0 : }
1501 0 : }
1502 : //
1503 : // Now do the FFT2D in place
1504 : //
1505 :
1506 : // storeImg(String("img.before.mod"), *image);
1507 0 : LatticeFFT::cfft2d(*lattice);
1508 : // storeImg(String("img.after.mod"), *image);
1509 :
1510 0 : log_l << LogIO::DEBUGGING << "Finished FFT" << LogIO::POST;
1511 0 : }
1512 : //
1513 : //---------------------------------------------------------------
1514 : //
1515 0 : void AWProjectFT::initializeToVis(ImageInterface<Complex>& iimage,
1516 : const VisBuffer& vb,
1517 : Array<Complex>& griddedVis,
1518 : Vector<Double>& uvscale)
1519 : {
1520 0 : initializeToVis(iimage, vb);
1521 0 : griddedVis.assign(griddedData); //using the copy for storage
1522 0 : uvscale.assign(uvScale);
1523 0 : }
1524 : //
1525 : //---------------------------------------------------------------
1526 : //
1527 0 : void AWProjectFT::finalizeToVis()
1528 : {
1529 0 : LogIO log_l(LogOrigin("AWProjectFT", "finalizeToVis[R&D]"));
1530 :
1531 : // cerr << "De-gridding run time = " << visResampler_p->runTimeDG_p << endl;
1532 0 : visResampler_p->runTimeDG_p=0.0;
1533 :
1534 0 : if(!arrayLattice.null()) arrayLattice=0;
1535 0 : if(!lattice.null()) lattice=0;
1536 0 : griddedData.resize();
1537 :
1538 0 : if(isTiled)
1539 : {
1540 0 : AlwaysAssert(imageCache, AipsError);
1541 0 : AlwaysAssert(image, AipsError);
1542 0 : ostringstream o;
1543 0 : imageCache->flush();
1544 0 : imageCache->showCacheStatistics(o);
1545 0 : log_l << o.str() << LogIO::POST;
1546 0 : }
1547 0 : if(pointingToImage) delete pointingToImage; pointingToImage=0;
1548 0 : }
1549 : //
1550 : //---------------------------------------------------------------
1551 : //
1552 : // Initialize the FFT to the Sky. Here we have to setup and
1553 : // initialize the grid.
1554 : //
1555 0 : void AWProjectFT::initializeToSky(ImageInterface<Complex>& iimage,
1556 : Matrix<Float>& weight,
1557 : const VisBuffer& vb)
1558 : {
1559 0 : LogIO log_l(LogOrigin("AWProjectFT", "initializeToSky[R&D]"));
1560 :
1561 : // image always points to the image
1562 0 : image=&iimage;
1563 :
1564 0 : init();
1565 0 : initMaps(vb);
1566 0 : visResampler_p->setMaps(chanMap, polMap);
1567 0 : visResampler_p->setFreqMaps(expandedSpwFreqSel_p,expandedSpwConjFreqSel_p);
1568 :
1569 : // Initialize the maps for polarization and channel. These maps
1570 : // translate visibility indices into image indices
1571 :
1572 0 : nx = image->shape()(0);
1573 0 : ny = image->shape()(1);
1574 0 : npol = image->shape()(2);
1575 0 : nchan = image->shape()(3);
1576 :
1577 0 : if(image->shape().product()>cachesize) isTiled=true;
1578 0 : else isTiled=false;
1579 :
1580 0 : sumWeight=0.0;
1581 0 : sumCFWeight = 0.0;
1582 0 : weight.resize(sumWeight.shape());
1583 0 : weight=0.0;
1584 : //
1585 : // Initialize for in memory or to disk gridding. lattice will
1586 : // point to the appropriate Lattice, either the ArrayLattice for
1587 : // in memory gridding or to the image for to disk gridding.
1588 : //
1589 0 : if(isTiled)
1590 : {
1591 0 : imageCache->flush();
1592 0 : image->set(Complex(0.0));
1593 0 : lattice=CountedPtr<Lattice<Complex> > (image, false);
1594 : }
1595 : else
1596 : {
1597 0 : IPosition gridShape(4, nx, ny, npol, nchan);
1598 0 : if(!useDoubleGrid_p){
1599 0 : griddedData.resize(gridShape);
1600 0 : griddedData=Complex(0.0);
1601 : }
1602 : // arrayLattice = new ArrayLattice<Complex>(griddedData);
1603 : // lattice=arrayLattice;
1604 : else {
1605 : //griddedData.resize();
1606 0 : griddedData2.resize(gridShape);
1607 0 : griddedData2=DComplex(0.0);
1608 : }
1609 0 : }
1610 :
1611 : //cerr << "initializeToSky for grid" << endl;
1612 0 : if(useDoubleGrid_p)
1613 0 : visResampler_p->initializeToSky(griddedData2, sumWeight);
1614 : else
1615 0 : visResampler_p->initializeToSky(griddedData, sumWeight);
1616 0 : }
1617 : //
1618 : //---------------------------------------------------------------
1619 : //
1620 0 : void AWProjectFT::finalizeToSky()
1621 : {
1622 : //
1623 : // Now we flush the cache and report statistics For memory based,
1624 : // we don't write anything out yet.
1625 : //
1626 0 : LogIO log_l(LogOrigin("AWProjectFT", "finalizeToSky[R&D]"));
1627 :
1628 0 : if(isTiled)
1629 : {
1630 0 : AlwaysAssert(image, AipsError);
1631 0 : AlwaysAssert(imageCache, AipsError);
1632 0 : imageCache->flush();
1633 0 : ostringstream o;
1634 0 : imageCache->showCacheStatistics(o);
1635 0 : log_l << o.str() << LogIO::POST;
1636 0 : }
1637 0 : if(pointingToImage) delete pointingToImage; pointingToImage=0;
1638 :
1639 0 : paChangeDetector.reset();
1640 0 : cfCache_p->flush();
1641 0 : if(useDoubleGrid_p)
1642 0 : visResampler_p->finalizeToSky(griddedData2, sumWeight);
1643 : else
1644 0 : visResampler_p->finalizeToSky(griddedData, sumWeight);
1645 0 : }
1646 : //
1647 : //---------------------------------------------------------------
1648 : //
1649 0 : Array<Complex>* AWProjectFT::getDataPointer(const IPosition& centerLoc2D,
1650 : Bool readonly)
1651 : {
1652 : Array<Complex>* result;
1653 : // Is tiled: get tiles and set up offsets
1654 0 : centerLoc(0)=centerLoc2D(0);
1655 0 : centerLoc(1)=centerLoc2D(1);
1656 0 : result=&imageCache->tile(offsetLoc, centerLoc, readonly);
1657 0 : gridder->setOffset(IPosition(2, offsetLoc(0), offsetLoc(1)));
1658 0 : return result;
1659 : }
1660 :
1661 : // The following file has the runFORTRAN* stuff. Moving it to a
1662 : // separate file to reduce clutter and ultimately delete it.
1663 : #include "AWProjectFT.FORTRANSTUFF"
1664 : //
1665 : //---------------------------------------------------------------
1666 : //
1667 0 : void AWProjectFT::put(const VisBuffer& vb, Int row, Bool dopsf,
1668 : FTMachine::Type type)
1669 : {
1670 0 : LogIO log_l(LogOrigin("AWProjectFT", "put[R&D]"));
1671 : // Take care of translation of Bools to Integer
1672 0 : makingPSF=dopsf;
1673 :
1674 : try
1675 : {
1676 0 : findConvFunction(*image, vb);
1677 : }
1678 0 : catch(AipsError& x)
1679 : {
1680 0 : log_l << x.getMesg() << LogIO::WARN;
1681 0 : return;
1682 0 : }
1683 0 : if (isDryRun) return;
1684 :
1685 0 : Nant_p = vb.msColumns().antenna().nrow();
1686 :
1687 : const Matrix<Float> *imagingweight;
1688 0 : imagingweight=&(vb.imagingWeight());
1689 :
1690 0 : Cube<Complex> data;
1691 : //Fortran gridder need the flag as ints
1692 0 : Cube<Int> flags;
1693 0 : Matrix<Float> elWeight;
1694 0 : interpolateFrequencyTogrid(vb, *imagingweight,data, flags, elWeight, type);
1695 : // NAnt no longer appears to be used
1696 : //Int NAnt;
1697 : //if (doPointing) NAnt = findPointingOffsets(vb,l_offsets,m_offsets,true);
1698 : //NAnt=NAnt; // Dummy statement to supress complier warnings and will be used when pointing offsets are used.
1699 : //
1700 : // If row is -1 then we pass through all rows
1701 : //
1702 : Int startRow, endRow, nRow;
1703 0 : if (row==-1)
1704 : {
1705 0 : nRow=vb.nRow();
1706 0 : startRow=0;
1707 0 : endRow=nRow-1;
1708 : }
1709 : else
1710 : {
1711 0 : nRow=1;
1712 0 : startRow=row;
1713 0 : endRow=row;
1714 : }
1715 : //
1716 : // Get the uvws in a form that Fortran can use and do that
1717 : // necessary phase rotation. On a Pentium Pro 200 MHz when null,
1718 : // this step takes about 50us per uvw point. This is just barely
1719 : // noticeable for Stokes I continuum and irrelevant for other
1720 : // cases.
1721 : //
1722 0 : Matrix<Double> uvw(3, vb.uvw().nelements());
1723 0 : uvw=0.0;
1724 0 : Vector<Double> dphase(vb.uvw().nelements());
1725 0 : dphase=0.0;
1726 : //NEGATING to correct for an image inversion problem
1727 0 : for (Int i=startRow;i<=endRow;i++)
1728 : {
1729 0 : for (Int idim=0;idim<2;idim++) uvw(idim,i)=-vb.uvw()(i)(idim);
1730 0 : uvw(2,i)=vb.uvw()(i)(2);
1731 : }
1732 :
1733 0 : doUVWRotation_p=true;
1734 : //rotateUVW(uvw, dphase, vb);
1735 0 : girarUVW(uvw, dphase, vb);
1736 0 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
1737 :
1738 : // // This is the convention for dphase
1739 : // dphase*=-1.0;
1740 :
1741 : //Check if ms has changed then cache new spw and chan selection
1742 0 : if(vb.newMS())
1743 0 : matchAllSpwChans(vb);
1744 :
1745 : //Here we redo the match or use previous match
1746 :
1747 : //Channel matching for the actual spectral window of buffer
1748 0 : if(doConversion_p[vb.spectralWindow()])
1749 0 : matchChannel(vb.spectralWindow(), vb);
1750 : else
1751 : {
1752 0 : chanMap.resize();
1753 0 : chanMap=multiChanMap_p[vb.spectralWindow()];
1754 : }
1755 :
1756 0 : VBStore vbs;
1757 0 : Vector<Int> gridShape = griddedData2.shape().asVector();
1758 0 : setupVBStore(vbs,vb, elWeight,data,uvw,flags, dphase,dopsf,gridShape);
1759 :
1760 0 : if (useDoubleGrid_p)
1761 : {
1762 0 : resampleDataToGrid(griddedData2, vbs, vb, dopsf);//, *imagingweight, *data, uvw,flags,dphase,dopsf);
1763 : //if (!dopsf) storeArrayAsImage(String("cgrid_awp_d.im"), image->coordinates(), griddedData2);
1764 : }
1765 : else
1766 : {
1767 0 : resampleDataToGrid(griddedData, vbs, vb, dopsf);//, *imagingweight, *data, uvw,flags,dphase,dopsf);
1768 : // String msg("WTH");
1769 : // isVBNaN(vb,msg);
1770 : //if (!dopsf) storeArrayAsImage(String("cgrid_awp_s.im"), image->coordinates(), griddedData);
1771 : }
1772 0 : }
1773 : //
1774 : //-------------------------------------------------------------------------
1775 : // Gridding
1776 0 : void AWProjectFT::resampleDataToGrid(Array<Complex>& griddedData_l, VBStore& vbs,
1777 : const VisBuffer& /*vb*/, Bool& dopsf)
1778 : {
1779 0 : LogIO log_l(LogOrigin("AWProjectFT", "resampleDataToGrid[R&D]"));
1780 0 : visResampler_p->DataToGrid(griddedData_l, vbs, sumWeight, dopsf);
1781 : // cerr << "####SumWt(C): " << sumWeight << endl;
1782 0 : }
1783 : //
1784 : //-------------------------------------------------------------------------
1785 : // Gridding
1786 0 : void AWProjectFT::resampleDataToGrid(Array<DComplex>& griddedData_l, VBStore& vbs,
1787 : const VisBuffer& /*vb*/, Bool& dopsf)
1788 : {
1789 0 : LogIO log_l(LogOrigin("AWProjectFT", "resampleDataToGridD[R&D]"));
1790 0 : visResampler_p->DataToGrid(griddedData_l, vbs, sumWeight, dopsf);
1791 : // cerr << "####SumWt(DC): " << sumWeight << endl;
1792 0 : }
1793 : //
1794 : //---------------------------------------------------------------
1795 : //
1796 0 : void AWProjectFT::get(VisBuffer& vb, Int row)
1797 : {
1798 0 : LogIO log_l(LogOrigin("AWProjectFT", "get[R&D]"));
1799 : // If row is -1 then we pass through all rows
1800 :
1801 : //------COMMON FROM HERE
1802 : Int startRow, endRow, nRow;
1803 0 : if (row==-1)
1804 : {
1805 0 : nRow=vb.nRow();
1806 0 : startRow=0;
1807 0 : endRow=nRow-1;
1808 0 : vb.modelVisCube()=Complex(0.0,0.0);
1809 : }
1810 : else
1811 : {
1812 0 : nRow=1;
1813 0 : startRow=row;
1814 0 : endRow=row;
1815 0 : vb.modelVisCube().xyPlane(row)=Complex(0.0,0.0);
1816 : }
1817 :
1818 0 : findConvFunction(*image, vb);
1819 :
1820 0 : Nant_p = vb.msColumns().antenna().nrow();
1821 : // NAnt set but not used
1822 : //Int NAnt=0;
1823 : //if (doPointing) NAnt = findPointingOffsets(vb,l_offsets,m_offsets,true);
1824 :
1825 : // Get the uvws in a form that Fortran can use
1826 0 : Matrix<Double> uvw(3, vb.uvw().nelements());
1827 0 : uvw=0.0;
1828 0 : Vector<Double> dphase(vb.uvw().nelements());
1829 0 : dphase=0.0;
1830 : //NEGATING to correct for an image inversion problem
1831 0 : for (Int i=startRow;i<=endRow;i++)
1832 : {
1833 0 : for (Int idim=0;idim<2;idim++) uvw(idim,i)=-vb.uvw()(i)(idim);
1834 0 : uvw(2,i)=vb.uvw()(i)(2);
1835 : }
1836 :
1837 0 : doUVWRotation_p=true;
1838 : //rotateUVW(uvw, dphase, vb);
1839 0 : girarUVW(uvw, dphase, vb);
1840 0 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
1841 :
1842 : // This is the convention for dphase
1843 : // dphase*=-1.0;
1844 :
1845 : //Check if ms has changed then cache new spw and chan selection
1846 0 : if(vb.newMS())
1847 0 : matchAllSpwChans(vb);
1848 :
1849 : //Here we redo the match or use previous match
1850 : //
1851 : //Channel matching for the actual spectral window of buffer
1852 : //
1853 0 : if(doConversion_p[vb.spectralWindow()])
1854 0 : matchChannel(vb.spectralWindow(), vb);
1855 : else
1856 : {
1857 0 : chanMap.resize();
1858 0 : chanMap=multiChanMap_p[vb.spectralWindow()];
1859 : }
1860 : //No point in reading data if its not matching in frequency
1861 0 : if(max(chanMap)==-1) return;
1862 : //-----COMMONG TILL HERE
1863 :
1864 : // Cube<Int> flags(vb.flagCube().shape());
1865 : // flags=0;
1866 : // flags(vb.flagCube())=true;
1867 :
1868 0 : Cube<Complex> data;
1869 0 : Cube<Int> flags;
1870 0 : getInterpolateArrays(vb, data, flags);
1871 :
1872 0 : VBStore vbs;
1873 : // setupVBStore(vbs,vb, vb.imagingWeight(),vb.modelVisCube(),uvw,flags, dphase, dopsf);
1874 0 : Bool tmpDoPSF=false;
1875 :
1876 0 : setupVBStore(vbs,vb, vb.imagingWeight(),data,uvw,flags, dphase,tmpDoPSF,griddedData.shape().asVector());
1877 :
1878 : // visResampler_p->setParams(uvScale,uvOffset,dphase);
1879 : // visResampler_p->setMaps(chanMap, polMap);
1880 0 : resampleGridToData(vbs, griddedData, vb);//, uvw, flags, dphase);
1881 :
1882 : // De-gridding
1883 : // visResampler_p->GridToData(vbs, griddedData);
1884 :
1885 :
1886 0 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1887 0 : }
1888 : //
1889 : //-------------------------------------------------------------------------
1890 : // De-gridding
1891 0 : void AWProjectFT::resampleGridToData(VBStore& vbs, Array<Complex>& griddedData_l,
1892 : const VisBuffer& /*vb*/)
1893 : {
1894 0 : LogIO log_l(LogOrigin("AWProjectFT", "resampleGridToData[R&D]"));
1895 : // Timer tim;
1896 : // tim.mark();
1897 0 : visResampler_p->GridToData(vbs, griddedData_l);
1898 : // runTime+=tim.user();
1899 0 : }
1900 : //
1901 : //---------------------------------------------------------------
1902 : //
1903 : // Finalize the FFT to the Sky. Here we actually do the FFT and
1904 : // return the resulting image
1905 0 : ImageInterface<Complex>& AWProjectFT::getImage(Matrix<Float>& weights,
1906 : Bool fftNormalization)
1907 : {
1908 0 : LogIO log_l(LogOrigin("AWProjectFT", "getImage[R&D]"));
1909 : //
1910 : // There are three objects held by the FTMachine objects: (1)
1911 : // *image, (2) *lattice, and (3) griddedData.
1912 : //
1913 : // (1) is the Dirty Image (ImageInterface<Float>)
1914 : //
1915 : // (2) is the Lattice of the Dirty Image (Lattice<Float>)
1916 : //
1917 : // (3) appears to be a reference to (2). Since FFT of *lattice is
1918 : // done in place, griddedData (and *lattice) have the gridded data
1919 : // in them before, and the Dirty Image data after the FFT.
1920 : //
1921 : // Question: Why three objects for the same precise information?
1922 : // --SB (Dec. 2010).
1923 0 : AlwaysAssert(image, AipsError);
1924 :
1925 0 : weights.resize(sumWeight.shape());
1926 0 : convertArray(weights, sumWeight);
1927 :
1928 : //
1929 : // If the weights are all zero then we cannot normalize otherwise
1930 : // we don't care.
1931 : //
1932 0 : if(max(weights)==0.0)
1933 : log_l // UUU //<< LogIO::SEVERE
1934 0 : << "No useful data in " << name() << ". Weights all zero"
1935 0 : << LogIO::POST;
1936 : // UUU else
1937 : {
1938 : log_l << LogIO::DEBUGGING
1939 0 : << "Starting FFT and scaling of image" << LogIO::POST;
1940 : //
1941 : // x and y transforms (lattice has the gridded vis. Make the
1942 : // dirty images)
1943 : //
1944 0 : if (useDoubleGrid_p)
1945 : {
1946 0 : ArrayLattice<DComplex> darrayLattice(griddedData2);
1947 0 : LatticeFFT::cfft2d(darrayLattice,false);
1948 0 : griddedData.resize(griddedData2.shape());
1949 0 : convertArray(griddedData, griddedData2);
1950 0 : SynthesisUtilMethods::getResource("mem peak in getImage");
1951 :
1952 : //Don't need the double-prec grid anymore...
1953 0 : griddedData2.resize();
1954 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
1955 0 : lattice=arrayLattice;
1956 0 : }
1957 : else
1958 : {
1959 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
1960 : //cerr << "##### " << griddedData2.shape() << endl;
1961 : //if (!makingPSF) storeArrayAsImage(String("cgrid_awp.im"), image->coordinates(), griddedData);
1962 0 : lattice=arrayLattice;
1963 0 : LatticeFFT::cfft2d(*lattice,false);
1964 : }
1965 0 : const IPosition latticeShape = lattice->shape();
1966 : //
1967 : // Now normalize the dirty image.
1968 : //
1969 : // Since *lattice is not copied to *image till the end of this
1970 : // method, normalizeImage also needs to work with Lattices
1971 : // (rather than ImageInterface).
1972 : //
1973 : // //normalizeImage(*lattice,sumWeight,*avgPB_p,fftNormalization);
1974 : // normalizeImage(*lattice,sumWeight,*avgPB_p, *avgPBSq_p, fftNormalization);
1975 : //if (!makingPSF) storeImg(String("uvgrid0.im"), *image);
1976 :
1977 : // nx ny normalization from GridFT...
1978 : {
1979 0 : Int inx = lattice->shape()(0);
1980 0 : Int iny = lattice->shape()(1);
1981 0 : Vector<Complex> correction(inx);
1982 0 : correction=Complex(1.0, 0.0);
1983 : // Do the Grid-correction
1984 0 : IPosition cursorShape(4, inx, 1, 1, 1);
1985 0 : IPosition axisPath(4, 0, 1, 2, 3);
1986 0 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
1987 0 : LatticeIterator<Complex> lix(*lattice, lsx);
1988 0 : for(lix.reset();!lix.atEnd();lix++)
1989 : {
1990 0 : Int pol=lix.position()(2);
1991 0 : Int chan=lix.position()(3);
1992 : {
1993 0 : if(fftNormalization)
1994 : {
1995 0 : if(weights(pol,chan)!=0.0)
1996 : {
1997 0 : Complex rnorm(Float(inx)*Float(iny)/( weights(pol,chan) ));
1998 0 : lix.rwCursor()*=rnorm;
1999 : }
2000 : else
2001 0 : lix.woCursor()=0.0;
2002 : }
2003 : else
2004 : {
2005 0 : Complex rnorm(Float(inx)*Float(iny));
2006 0 : lix.rwCursor()*=rnorm;
2007 : }
2008 : }
2009 : }
2010 0 : }
2011 : //if (!makingPSF) storeImg(String("uvgrid.im"), *image);
2012 0 : if(!isTiled)
2013 : {
2014 : //
2015 : // Check the section from the image BEFORE converting to a lattice
2016 : //
2017 0 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
2018 0 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
2019 0 : IPosition stride(4, 1);
2020 0 : IPosition trc(blc+image->shape()-stride);
2021 : //
2022 : // Do the copy
2023 : //
2024 0 : image->put(griddedData(blc, trc));
2025 :
2026 :
2027 0 : if(!arrayLattice.null()) arrayLattice=0;
2028 0 : if(!lattice.null()) lattice=0;
2029 0 : griddedData.resize(IPosition(1,0));
2030 0 : }
2031 0 : }
2032 :
2033 : // {
2034 : // // TempImage<Complex> tt(lattice->shape(), image->coordinates());
2035 : // // tt.put(lattice->get());
2036 : // if (!makingPSF) throw(AipsError("Exiting from AWP.cc: 1765"));
2037 :
2038 : // }
2039 :
2040 0 : return *image;
2041 0 : }
2042 : //
2043 : //---------------------------------------------------------------
2044 : //
2045 : // Get weight image
2046 0 : void AWProjectFT::getWeightImage(ImageInterface<Float>& weightImage,
2047 : Matrix<Float>& weights)
2048 : {
2049 0 : LogIO log_l(LogOrigin("AWProjectFT", "getWeightImage[R&D]"));
2050 :
2051 0 : weights.resize(sumWeight.shape());
2052 0 : convertArray(weights,sumWeight);
2053 :
2054 0 : const IPosition latticeShape = weightImage.shape();
2055 0 : const IPosition avgpbShape = avgPB_p->shape();
2056 : // cout << "AWP::getWeightImage : weightimage shape : " << latticeShape << " and avgpb shape : " << avgpbShape << " " << sumWeight << endl;
2057 :
2058 0 : Int nx=latticeShape(0);
2059 0 : Int ny=latticeShape(1);
2060 :
2061 0 : IPosition loc(2, 0);
2062 0 : IPosition cursorShape(4, nx, ny, 1, 1);
2063 0 : IPosition axisPath(4, 0, 1, 2, 3);
2064 0 : LatticeStepper lsx(latticeShape, cursorShape, axisPath);
2065 0 : LatticeIterator<Float> lix(weightImage, lsx);
2066 0 : LatticeIterator<Float> liy(*avgPB_p,lsx);
2067 0 : for(lix.reset();!lix.atEnd();lix++)
2068 : {
2069 : //Int pol=lix.position()(2);
2070 : //Int chan=lix.position()(3);
2071 : //lix.rwCursor()=weights(pol,chan);
2072 0 : lix.rwCursor()=liy.rwCursor();
2073 : }
2074 0 : }
2075 : //
2076 : //---------------------------------------------------------------
2077 : //
2078 0 : Bool AWProjectFT::toRecord(RecordInterface& outRec, Bool withImage)
2079 : {
2080 0 : LogIO log_l(LogOrigin("AWProjectFT", "toRecord[R&D]"));
2081 :
2082 : // Save the current AWProjectFT object to an output state record
2083 0 : Bool retval = true;
2084 0 : Double cacheVal=(Double) cachesize;
2085 0 : outRec.define("cache", cacheVal);
2086 0 : outRec.define("tile", tilesize);
2087 :
2088 0 : Vector<Double> phaseValue(2);
2089 0 : String phaseUnit;
2090 0 : phaseValue=mTangent_p.getAngle().getValue();
2091 0 : phaseUnit= mTangent_p.getAngle().getUnit();
2092 0 : outRec.define("phasevalue", phaseValue);
2093 0 : outRec.define("phaseunit", phaseUnit);
2094 :
2095 0 : Vector<Double> dirValue(3);
2096 0 : String dirUnit;
2097 0 : dirValue=mLocation_p.get("m").getValue();
2098 0 : dirUnit=mLocation_p.get("m").getUnit();
2099 0 : outRec.define("dirvalue", dirValue);
2100 0 : outRec.define("dirunit", dirUnit);
2101 :
2102 0 : outRec.define("padding", padding_p);
2103 0 : outRec.define("maxdataval", maxAbsData);
2104 :
2105 0 : Vector<Int> center_loc(4), offset_loc(4);
2106 0 : for (Int k=0; k<4 ; k++)
2107 : {
2108 0 : center_loc(k)=centerLoc(k);
2109 0 : offset_loc(k)=offsetLoc(k);
2110 : }
2111 0 : outRec.define("centerloc", center_loc);
2112 0 : outRec.define("offsetloc", offset_loc);
2113 0 : outRec.define("sumofweights", sumWeight);
2114 0 : outRec.define("sumofcfweights", sumCFWeight);
2115 0 : if(withImage && image)
2116 : {
2117 0 : ImageInterface<Complex>& tempimage(*image);
2118 0 : Record imageContainer;
2119 0 : String error;
2120 0 : retval = (retval || tempimage.toRecord(error, imageContainer));
2121 0 : outRec.defineRecord("image", imageContainer);
2122 0 : }
2123 0 : return retval;
2124 0 : }
2125 : //
2126 : //---------------------------------------------------------------
2127 : //
2128 0 : Bool AWProjectFT::fromRecord(const RecordInterface& inRec)
2129 : {
2130 0 : LogIO log_l(LogOrigin("AWProjectFT", "fromRecord[R&D]"));
2131 :
2132 0 : Bool retval = true;
2133 0 : imageCache=0; lattice=0; arrayLattice=0;
2134 : Double cacheVal;
2135 0 : inRec.get("cache", cacheVal);
2136 0 : cachesize=(Long) cacheVal;
2137 0 : inRec.get("tile", tilesize);
2138 :
2139 0 : Vector<Double> phaseValue(2);
2140 0 : inRec.get("phasevalue",phaseValue);
2141 0 : String phaseUnit;
2142 0 : inRec.get("phaseunit",phaseUnit);
2143 0 : Quantity val1(phaseValue(0), phaseUnit);
2144 0 : Quantity val2(phaseValue(1), phaseUnit);
2145 0 : MDirection phasecenter(val1, val2);
2146 :
2147 0 : mTangent_p=phasecenter;
2148 : // This should be passed down too but the tangent plane is
2149 : // expected to be specified in all meaningful cases.
2150 0 : tangentSpecified_p=true;
2151 0 : Vector<Double> dirValue(3);
2152 0 : String dirUnit;
2153 0 : inRec.get("dirvalue", dirValue);
2154 0 : inRec.get("dirunit", dirUnit);
2155 0 : MVPosition dummyMVPos(dirValue(0), dirValue(1), dirValue(2));
2156 0 : MPosition mLocation(dummyMVPos, MPosition::ITRF);
2157 0 : mLocation_p=mLocation;
2158 :
2159 0 : inRec.get("padding", padding_p);
2160 0 : inRec.get("maxdataval", maxAbsData);
2161 :
2162 0 : Vector<Int> center_loc(4), offset_loc(4);
2163 0 : inRec.get("centerloc", center_loc);
2164 0 : inRec.get("offsetloc", offset_loc);
2165 0 : uInt ndim4 = 4;
2166 0 : centerLoc=IPosition(ndim4, center_loc(0), center_loc(1), center_loc(2),
2167 0 : center_loc(3));
2168 0 : offsetLoc=IPosition(ndim4, offset_loc(0), offset_loc(1), offset_loc(2),
2169 0 : offset_loc(3));
2170 0 : inRec.get("sumofweights", sumWeight);
2171 0 : inRec.get("sumofcfweights", sumCFWeight);
2172 0 : if(inRec.nfields() > 12 )
2173 : {
2174 0 : Record imageAsRec=inRec.asRecord("image");
2175 0 : if(!image) image= new TempImage<Complex>();
2176 :
2177 0 : String error;
2178 0 : retval = (retval || image->fromRecord(error, imageAsRec));
2179 :
2180 : // Might be changing the shape of sumWeight
2181 0 : init();
2182 :
2183 0 : if(isTiled)
2184 0 : lattice=CountedPtr<Lattice<Complex> > (image, false);
2185 : else
2186 : {
2187 : //
2188 : // Make the grid the correct shape and turn it into an
2189 : // array lattice Check the section from the image BEFORE
2190 : // converting to a lattice
2191 : //
2192 0 : IPosition gridShape(4, nx, ny, npol, nchan);
2193 0 : griddedData.resize(gridShape);
2194 0 : griddedData=Complex(0.0);
2195 0 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
2196 0 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
2197 0 : IPosition start(4, 0);
2198 0 : IPosition stride(4, 1);
2199 0 : IPosition trc(blc+image->shape()-stride);
2200 0 : griddedData(blc, trc) = image->getSlice(start, image->shape());
2201 :
2202 0 : arrayLattice = new ArrayLattice<Complex>(griddedData);
2203 0 : lattice=arrayLattice;
2204 0 : }
2205 :
2206 0 : AlwaysAssert(image, AipsError);
2207 0 : };
2208 0 : return retval;
2209 0 : }
2210 : //
2211 : //---------------------------------------------------------------
2212 : //
2213 0 : void AWProjectFT::ok()
2214 : {
2215 0 : AlwaysAssert(image, AipsError);
2216 0 : }
2217 : //----------------------------------------------------------------------
2218 : //
2219 : // Make a straightforward 2D interferometry honest-to-god
2220 : // image. This returns a complex image, without conversion to
2221 : // Stokes. The polarization representation is that required for the
2222 : // visibilities.
2223 : //
2224 0 : void AWProjectFT::makeImage(FTMachine::Type type,
2225 : VisSet& vs,
2226 : ImageInterface<Complex>& theImage,
2227 : Matrix<Float>& weight)
2228 : {
2229 0 : LogIO log_l(LogOrigin("AWProjectFT", "makeImage[R&D]"));
2230 :
2231 0 : if(type==FTMachine::COVERAGE)
2232 : log_l << "Type COVERAGE not defined for Fourier transforms"
2233 0 : << LogIO::EXCEPTION;
2234 :
2235 :
2236 : // Initialize the gradients
2237 0 : ROVisIter& vi(vs.iter());
2238 :
2239 : // Loop over all visibilities and pixels
2240 0 : VisBuffer vb(vi);
2241 :
2242 : // Initialize put (i.e. transform to Sky) for this model
2243 0 : vi.origin();
2244 :
2245 0 : if(vb.polFrame()==MSIter::Linear)
2246 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::LINEAR);
2247 : else
2248 0 : StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::CIRCULAR);
2249 :
2250 0 : initializeToSky(theImage,weight,vb);
2251 :
2252 : //
2253 : // Loop over the visibilities, putting VisBuffers
2254 : //
2255 0 : paChangeDetector.reset();
2256 :
2257 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk())
2258 : {
2259 0 : for (vi.origin(); vi.more(); vi++)
2260 : {
2261 0 : if (type==FTMachine::PSF) makingPSF=true;
2262 0 : findConvFunction(theImage,vb);
2263 :
2264 0 : switch(type)
2265 : {
2266 0 : case FTMachine::RESIDUAL:
2267 0 : vb.visCube()=vb.correctedVisCube();
2268 0 : vb.visCube()-=vb.modelVisCube();
2269 0 : put(vb, -1, false);
2270 0 : break;
2271 0 : case FTMachine::MODEL:
2272 0 : vb.visCube()=vb.modelVisCube();
2273 0 : put(vb, -1, false);
2274 0 : break;
2275 0 : case FTMachine::CORRECTED:
2276 0 : vb.visCube()=vb.correctedVisCube();
2277 0 : put(vb, -1, false);
2278 0 : break;
2279 0 : case FTMachine::PSF:
2280 0 : vb.visCube()=Complex(1.0,0.0);
2281 0 : makingPSF = true;
2282 0 : put(vb, -1, true);
2283 0 : break;
2284 0 : case FTMachine::OBSERVED:
2285 : default:
2286 0 : put(vb, -1, false);
2287 0 : break;
2288 : }
2289 : }
2290 : }
2291 0 : finalizeToSky();
2292 : // Normalize by dividing out weights, etc.
2293 0 : getImage(weight, true);
2294 0 : }
2295 : //
2296 : //-------------------------------------------------------------------------
2297 : //
2298 0 : void AWProjectFT::setPAIncrement(const Quantity& computePAIncrement,
2299 : const Quantity& rotateOTFPAIncrement)
2300 : {
2301 0 : LogIO log_l(LogOrigin("AWProjectFT", "setPAIncrement[R&D]"));
2302 :
2303 : // Quantity tmp(abs(paIncrement.getValue("deg")), "deg");
2304 : // if (paIncrement.getValue("rad") < 0)
2305 : // {
2306 : // rotateApertureOTF_p = false;
2307 : // convFuncCtor_p->setRotateCF(tmp.getValue("rad"), rotateApertureOTF_p);
2308 : // }
2309 : // else
2310 : // {
2311 : // rotateApertureOTF_p = true;
2312 : // convFuncCtor_p->setRotateCF(360.0*M_PI/180.0, rotateApertureOTF_p);
2313 : // }
2314 :
2315 :
2316 0 : rotateOTFPAIncr_p = rotateOTFPAIncrement.getValue("rad");
2317 0 : computePAIncr_p = computePAIncrement.getValue("rad");
2318 0 : convFuncCtor_p->setRotateCF(computePAIncr_p, rotateOTFPAIncr_p);
2319 :
2320 0 : paChangeDetector.setTolerance(computePAIncrement);
2321 0 : visResampler_p->setPATolerance(computePAIncrement.getValue("rad"));
2322 : log_l << LogIO::NORMAL <<"Setting PA increment to "
2323 0 : << computePAIncrement.getValue("deg") << " deg" << endl;
2324 0 : cfCache_p->setPAChangeDetector(paChangeDetector);
2325 0 : }
2326 : //
2327 : //-------------------------------------------------------------------------
2328 : //
2329 0 : Bool AWProjectFT::verifyShapes(IPosition pbShape, IPosition skyShape)
2330 : {
2331 0 : LogIO log_l(LogOrigin("AWProjectFT", "verifyShapes[R&D]"));
2332 :
2333 0 : if ((pbShape(0) != skyShape(0)) && // X-axis
2334 0 : (pbShape(1) != skyShape(1)) && // Y-axis
2335 0 : (pbShape(2) != skyShape(2))) // Poln-axis
2336 : {
2337 : log_l << "Sky and/or polarization shape of the avgPB and the sky model do not match."
2338 0 : << LogIO::EXCEPTION;
2339 0 : return false;
2340 : }
2341 0 : return true;
2342 :
2343 0 : }
2344 : //
2345 : //-------------------------------------------------------------------------
2346 : //
2347 0 : void AWProjectFT::makeThGridCoords(VBStore& vbs, const Vector<Int>& gridShape)
2348 : {
2349 0 : Float dNx=(gridShape(0)/8.0), dNy=(gridShape(1)/8.0);
2350 0 : Int GNx=SynthesisUtils::nint(gridShape(0)/dNx),
2351 0 : GNy=SynthesisUtils::nint(gridShape(1)/dNy);
2352 :
2353 0 : vbs.BLCXi.resize(GNx,GNy);
2354 0 : vbs.BLCYi.resize(GNx,GNy);
2355 0 : vbs.TRCXi.resize(GNx,GNy);
2356 0 : vbs.TRCYi.resize(GNx,GNy);
2357 0 : vbs.BLCXi=vbs.BLCYi=vbs.TRCXi=vbs.TRCYi=-1;
2358 :
2359 0 : for (Int i=0;i<GNx;i++)
2360 0 : for (Int j=0;j<GNy;j++)
2361 : {
2362 0 : vbs.BLCXi(i,j)=max(0,SynthesisUtils::nint(i*dNx));
2363 0 : vbs.BLCYi(i,j)=max(0,SynthesisUtils::nint(j*dNy));
2364 0 : vbs.TRCXi(i,j)=min(gridShape(0), SynthesisUtils::nint((i+1)*dNx-1.0));
2365 0 : vbs.TRCYi(i,j)=min(gridShape(1), SynthesisUtils::nint((j+1)*dNy-1.0));
2366 : }
2367 : // for (Int i=0;i<GNx;i++)
2368 : // for (Int j=0;j<GNy;j++)
2369 : // cerr << i << " " << j << ": " << vbs.BLCXi(i,j) << " " << vbs.BLCYi(i,j) << " " << vbs.TRCXi(i,j) << " " << vbs.TRCYi(i,j) << endl;
2370 0 : }
2371 : //
2372 : //-------------------------------------------------------------------------
2373 : //
2374 0 : void AWProjectFT::setupVBStore(VBStore& vbs,
2375 : const VisBuffer& vb,
2376 : const Matrix<Float>& imagingweight,
2377 : const Cube<Complex>& visData,
2378 : const Matrix<Double>& uvw,
2379 : const Cube<Int>& flagCube,
2380 : const Vector<Double>& dphase,
2381 : const Bool& dopsf,
2382 : const Vector<Int>& /*gridShape*/)
2383 : {
2384 0 : LogIO log_l(LogOrigin("AWProjectFT", "setupVBStore[R&D]"));
2385 :
2386 : // Vector<Int> ConjCFMap, CFMap;
2387 :
2388 0 : vbs.vb_p = &vb;
2389 0 : makeCFPolMap(vb,cfStokes_p,CFMap_p);
2390 0 : makeConjPolMap(vb,CFMap_p,ConjCFMap_p);
2391 :
2392 0 : visResampler_p->setParams(uvScale,uvOffset,dphase);
2393 0 : visResampler_p->setMaps(chanMap, polMap);
2394 0 : visResampler_p->setCFMaps(CFMap_p, ConjCFMap_p);
2395 0 : visResampler_p->setFreqMaps(expandedSpwFreqSel_p,expandedSpwConjFreqSel_p);
2396 : //
2397 : // Set up VBStore object to point to the relavent info. of the VB.
2398 : //
2399 0 : vbs.imRefFreq_p = imRefFreq_p;
2400 0 : vbs.nRow_p = vb.nRow();
2401 0 : vbs.beginRow_p = 0;
2402 0 : vbs.endRow_p = vbs.nRow_p;
2403 0 : vbs.spwID_p = vb.spectralWindow();
2404 0 : vbs.nDataPol_p = flagCube.shape()[0];
2405 0 : vbs.nDataChan_p = flagCube.shape()[1];
2406 :
2407 0 : vbs.antenna1_p.reference(vb.antenna1());
2408 0 : vbs.antenna2_p.reference(vb.antenna2());
2409 0 : vbs.paQuant_p = Quantity(getPA(vb),"rad");
2410 0 : vbs.corrType_p.reference(vb.corrType());
2411 0 : vbs.uvw_p.reference(uvw);
2412 0 : vbs.imagingWeight_p.reference(imagingweight);
2413 0 : vbs.visCube_p.reference(visData);
2414 : // vbs.freq_p.reference(interpVisFreq_p);
2415 0 : vbs.freq_p.reference(vb.frequency());
2416 : // vbs.rowFlag_p.assign(vb.flagRow());
2417 0 : vbs.rowFlag_p.reference(vb.flagRow());
2418 0 : if(!usezero_p)
2419 0 : for (Int rownr=0; rownr<vbs.nRow_p; rownr++)
2420 0 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) vbs.rowFlag_p(rownr)=true;
2421 :
2422 : // Really nice way of converting a Cube<Int> to Cube<Bool>.
2423 : // However the VBS objects should ultimately be references
2424 : // directly to bool cubes.
2425 : // vbs.rowFlag.resize(rowFlags.shape()); vbs.rowFlag = false; vbs.rowFlag(rowFlags) = true;
2426 0 : vbs.flagCube_p.resize(flagCube.shape()); vbs.flagCube_p = false; vbs.flagCube_p(flagCube!=0) = true;
2427 0 : vbs.conjBeams_p=conjBeams_p;
2428 :
2429 0 : timer_p.mark();
2430 :
2431 0 : Vector<Double> pointingOffset(convFuncCtor_p->findPointingOffset(*image, vb));
2432 0 : if (makingPSF)
2433 : {
2434 0 : cfwts2_p->invokeGC(vbs.spwID_p);
2435 : //SynthesisUtilMethods::getResource(String("AfterWTS_CFBClearing"),String("memusage"));
2436 :
2437 0 : visResampler_p->makeVBRow2CFMap(*cfwts2_p,*convFuncCtor_p, vb,
2438 0 : paChangeDetector.getParAngleTolerance(),
2439 0 : chanMap,polMap,pointingOffset);
2440 : //SynthesisUtilMethods::getResource(String("PSFGridding"),String("memusage"));
2441 : }
2442 : else
2443 : {
2444 : // If the Wt. CFs are still in the memory, clear them. They
2445 : // won't be required again (though with the silly check below,
2446 : // if the in-memory Wt. CFs are less than 1KB, they will be
2447 : // left in memory).
2448 0 : if (cfwts2_p->memUsage() > 1000)
2449 : {
2450 : //SynthesisUtilMethods::getResource(String("BeforeWTSClearing"),String("memusage"));
2451 : //cerr << "CFWTS memusage before clearing: " << cfwts2_p->memUsage() << endl;
2452 :
2453 0 : cfwts2_p->clear();
2454 :
2455 : //cerr << "CFWTS memusage after clearing: " << cfwts2_p->memUsage() << endl;
2456 : //SynthesisUtilMethods::getResource(String("AfterWTSClearing"),String("memusage"));
2457 : }
2458 :
2459 0 : cfs2_p->invokeGC(vbs.spwID_p);
2460 : //SynthesisUtilMethods::getResource(String("After_CFBClearing"),String("memusage"));
2461 0 : visResampler_p->makeVBRow2CFMap(*cfs2_p,*convFuncCtor_p, vb,
2462 0 : paChangeDetector.getParAngleTolerance(),
2463 0 : chanMap,polMap,pointingOffset);
2464 : //SynthesisUtilMethods::getResource(String("ImGridding"),String("memusage"));
2465 : }
2466 :
2467 :
2468 : // VBRow2CFMapType theMap(visResampler_p->getVBRow2CFMap());
2469 0 : VBRow2CFBMapType& theMap=visResampler_p->getVBRow2CFBMap();
2470 : //
2471 : // For AzElApertures, this rotates the CFs.
2472 : //
2473 0 : convFuncCtor_p->prepareConvFunction(vb,theMap);
2474 :
2475 :
2476 : // CFBStruct cfbst_pub;
2477 : //UUU theMap[0]->getAsStruct(cfbst_pub);
2478 : //UUU vbs.cfBSt_p=cfbst_pub;
2479 0 : vbs.accumCFs_p=((vbs.uvw_p.nelements() == 0) && dopsf);
2480 :
2481 :
2482 : // for(int ii=0;ii<cfbst_pub.shape[0];ii++)
2483 : // for(int jj=0;jj<cfbst_pub.shape[1];jj++)
2484 : // for(int kk=0;kk<cfbst_pub.shape[2];kk++)
2485 : // {
2486 : // CFCStruct cfcst=cfbst_pub.getCFB(ii,jj,kk);
2487 : // cerr << "[" << ii << "," << jj << "," << kk << "]:"
2488 : // << cfcst.sampling << " "
2489 : // << cfcst.xSupport << " "
2490 : // << cfcst.ySupport
2491 : // << endl;
2492 : // }
2493 :
2494 :
2495 : // The following code is required only for GPU or multi-threaded
2496 : //gridder. Currently does not work without the rest of the
2497 : //GPU/multi-threaded infrastructure (though, I (SB) thought this
2498 : //was designed to be benign for normal gridding).
2499 : //
2500 : //makeThGridCoords(vbs,gridShape);
2501 :
2502 0 : runTime1_p += timer_p.real();
2503 0 : visResampler_p->initializeDataBuffers(vbs);
2504 : // visResampler_p->setConvFunc(cfs_p);
2505 0 : }
2506 :
2507 : //
2508 : //---------------------------------------------------------------
2509 : //
2510 : // Predict the coherences as well as their derivatives w.r.t. the
2511 : // pointing offsets.
2512 : //
2513 0 : void AWProjectFT::nget(VisBuffer& vb,
2514 : // These offsets should be appropriate for the VB
2515 : Array<Float>& l_off, Array<Float>& m_off,
2516 : Cube<Complex>& Mout,
2517 : Cube<Complex>& dMout1,
2518 : Cube<Complex>& dMout2,
2519 : Int Conj, Int doGrad)
2520 : {
2521 0 : LogIO log_l(LogOrigin("AWProjectFT", "nget[R&D]"));
2522 : Int startRow, endRow, nRow;
2523 0 : nRow=vb.nRow();
2524 0 : startRow=0;
2525 0 : endRow=nRow-1;
2526 :
2527 0 : Mout = dMout1 = dMout2 = Complex(0,0);
2528 :
2529 0 : findConvFunction(*image, vb);
2530 0 : Int NAnt=0;
2531 0 : Nant_p = vb.msColumns().antenna().nrow();
2532 0 : if (doPointing)
2533 0 : NAnt = findPointingOffsets(vb,l_offsets,m_offsets,false);
2534 :
2535 0 : l_offsets=l_off;
2536 0 : m_offsets=m_off;
2537 0 : Matrix<Double> uvw(3, vb.uvw().nelements());
2538 0 : uvw=0.0;
2539 0 : Vector<Double> dphase(vb.uvw().nelements());
2540 0 : dphase=0.0;
2541 : //NEGATING to correct for an image inversion problem
2542 0 : for (Int i=startRow;i<=endRow;i++)
2543 : {
2544 0 : for (Int idim=0;idim<2;idim++) uvw(idim,i)=-vb.uvw()(i)(idim);
2545 0 : uvw(2,i)=vb.uvw()(i)(2);
2546 : }
2547 :
2548 0 : doUVWRotation_p=true;
2549 : //rotateUVW(uvw, dphase, vb);
2550 0 : girarUVW(uvw, dphase, vb);
2551 0 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
2552 :
2553 : // This is the convention for dphase
2554 : // dphase*=-1.0;
2555 :
2556 0 : Cube<Int> flags(vb.flagCube().shape());
2557 0 : flags=0;
2558 0 : flags(vb.flagCube())=true;
2559 :
2560 : //Check if ms has changed then cache new spw and chan selection
2561 0 : if(vb.newMS())
2562 0 : matchAllSpwChans(vb);
2563 :
2564 : //Here we redo the match or use previous match
2565 : //
2566 : //Channel matching for the actual spectral window of buffer
2567 : //
2568 0 : if(doConversion_p[vb.spectralWindow()])
2569 0 : matchChannel(vb.spectralWindow(), vb);
2570 : else
2571 : {
2572 0 : chanMap.resize();
2573 0 : chanMap=multiChanMap_p[vb.spectralWindow()];
2574 : }
2575 :
2576 0 : Vector<Int> rowFlags(vb.nRow());
2577 0 : rowFlags=0;
2578 0 : rowFlags(vb.flagRow())=true;
2579 0 : if(!usezero_p)
2580 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
2581 0 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
2582 :
2583 0 : IPosition s,gradS;
2584 0 : Cube<Complex> visdata,gradVisAzData,gradVisElData;
2585 : //
2586 : // visdata now references the Mout data structure rather than to the internal VB storeage.
2587 : //
2588 0 : visdata.reference(Mout);
2589 :
2590 0 : if (doGrad)
2591 : {
2592 : // The following should reference some slice of dMout?
2593 0 : gradVisAzData.reference(dMout1);
2594 0 : gradVisElData.reference(dMout2);
2595 : }
2596 0 : visResampler_p->setParams(uvScale,uvOffset,dphase);
2597 0 : visResampler_p->setMaps(chanMap, polMap);
2598 0 : visResampler_p->setFreqMaps(expandedSpwFreqSel_p,expandedSpwConjFreqSel_p);
2599 : // Vector<Int> ConjCFMap, CFMap;
2600 0 : makeCFPolMap(vb,cfStokes_p,CFMap_p);
2601 0 : makeConjPolMap(vb,CFMap_p,ConjCFMap_p);
2602 0 : visResampler_p->setCFMaps(CFMap_p, ConjCFMap_p);
2603 : //
2604 : // Begin the actual de-gridding.
2605 : //
2606 0 : if(isTiled)
2607 : {
2608 0 : log_l << "The sky model is tiled" << LogIO::NORMAL << LogIO::POST;
2609 0 : Double invLambdaC=vb.frequency()(0)/C::c;
2610 0 : Vector<Double> uvLambda(2);
2611 0 : Vector<Int> centerLoc2D(2);
2612 0 : centerLoc2D=0;
2613 :
2614 : // Loop over all rows
2615 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
2616 : {
2617 :
2618 : // Calculate uvw for this row at the center frequency
2619 0 : uvLambda(0)=uvw(0, rownr)*invLambdaC;
2620 0 : uvLambda(1)=uvw(1, rownr)*invLambdaC;
2621 0 : centerLoc2D=gridder->location(centerLoc2D, uvLambda);
2622 :
2623 : // Is this point on the grid?
2624 0 : if(gridder->onGrid(centerLoc2D))
2625 : {
2626 :
2627 : // Get the tile
2628 0 : Array<Complex>* dataPtr=getDataPointer(centerLoc2D, true);
2629 0 : gridder->setOffset(IPosition(2, offsetLoc(0), offsetLoc(1)));
2630 0 : Int aNx=dataPtr->shape()(0);
2631 0 : Int aNy=dataPtr->shape()(1);
2632 :
2633 : // Now use FORTRAN to do the gridding. Remember to
2634 : // ensure that the shape and offsets of the tile are
2635 : // accounted for.
2636 :
2637 0 : Vector<Double> actualOffset(3);
2638 0 : for (Int i=0;i<2;i++)
2639 0 : actualOffset(i)=uvOffset(i)-Double(offsetLoc(i));
2640 :
2641 0 : actualOffset(2)=uvOffset(2);
2642 0 : IPosition s(vb.modelVisCube().shape());
2643 :
2644 0 : Int ScanNo=0, tmpPAI;
2645 0 : Double area=1.0;
2646 0 : tmpPAI = 1;
2647 0 : runFortranGetGrad(uvw,dphase,visdata,s,
2648 : gradVisAzData,gradVisElData,
2649 : Conj,flags,rowFlags,rownr,
2650 0 : actualOffset,dataPtr,aNx,aNy,npol,nchan,vb,NAnt,ScanNo,sigma,
2651 0 : l_offsets,m_offsets,area,doGrad,tmpPAI);
2652 0 : }
2653 : }
2654 0 : }
2655 : else
2656 : {
2657 0 : IPosition s(vb.modelVisCube().shape());
2658 0 : Int ScanNo=0, tmpPAI, trow=-1;
2659 0 : Double area=1.0;
2660 0 : tmpPAI = 1;
2661 0 : runFortranGetGrad(uvw,dphase,visdata/*vb.modelVisCube()*/,s,
2662 : gradVisAzData, gradVisElData,
2663 : Conj,flags,rowFlags,trow,
2664 0 : uvOffset,&griddedData,nx,ny,npol,nchan,vb,NAnt,ScanNo,sigma,
2665 0 : l_offsets,m_offsets,area,doGrad,tmpPAI);
2666 0 : }
2667 :
2668 0 : }
2669 0 : void AWProjectFT::get(VisBuffer& vb,
2670 : VisBuffer& gradVBAz,
2671 : VisBuffer& gradVBEl,
2672 : Cube<Float>& pointingOffsets,
2673 : Int row, // default row=-1
2674 : Type whichVBColumn, // default whichVBColumn = FTMachine::MODEL
2675 : Type whichGradVBColumn,// default whichGradVBColumn = FTMachine::MODEL
2676 : Int Conj, Int doGrad) // default Conj=0, doGrad=1
2677 : {
2678 : // If row is -1 then we pass through all rows
2679 : Int startRow, endRow, nRow;
2680 0 : if (row==-1)
2681 : {
2682 0 : nRow=vb.nRow();
2683 0 : startRow=0;
2684 0 : endRow=nRow-1;
2685 0 : initVisBuffer(vb,whichVBColumn);
2686 0 : if (doGrad)
2687 : {
2688 0 : initVisBuffer(gradVBAz, whichGradVBColumn);
2689 0 : initVisBuffer(gradVBEl, whichGradVBColumn);
2690 : }
2691 : }
2692 : else
2693 : {
2694 0 : nRow=1;
2695 0 : startRow=row;
2696 0 : endRow=row;
2697 0 : initVisBuffer(vb,whichVBColumn,row);
2698 0 : if (doGrad)
2699 : {
2700 0 : initVisBuffer(gradVBAz, whichGradVBColumn,row);
2701 0 : initVisBuffer(gradVBEl, whichGradVBColumn,row);
2702 : }
2703 : }
2704 :
2705 0 : findConvFunction(*image, vb);
2706 :
2707 0 : Nant_p = vb.msColumns().antenna().nrow();
2708 0 : Int NAnt=0;
2709 0 : if (doPointing)
2710 0 : NAnt = findPointingOffsets(vb,pointingOffsets,l_offsets,m_offsets,false);
2711 :
2712 0 : Matrix<Double> uvw(3, vb.uvw().nelements());
2713 0 : uvw=0.0;
2714 0 : Vector<Double> dphase(vb.uvw().nelements());
2715 0 : dphase=0.0;
2716 : //NEGATING to correct for an image inversion problem
2717 0 : for (Int i=startRow;i<=endRow;i++)
2718 : {
2719 0 : for (Int idim=0;idim<2;idim++) uvw(idim,i)=-vb.uvw()(i)(idim);
2720 0 : uvw(2,i)=vb.uvw()(i)(2);
2721 : }
2722 :
2723 0 : doUVWRotation_p=true;
2724 : //rotateUVW(uvw, dphase, vb);
2725 0 : girarUVW(uvw, dphase, vb);
2726 0 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
2727 :
2728 : // This is the convention for dphase
2729 : //dphase*=-1.0;
2730 :
2731 :
2732 0 : Cube<Int> flags(vb.flagCube().shape());
2733 0 : flags=0;
2734 0 : flags(vb.flagCube())=true;
2735 : //
2736 : //Check if ms has changed then cache new spw and chan selection
2737 : //
2738 0 : if(vb.newMS()) matchAllSpwChans(vb);
2739 :
2740 : //Here we redo the match or use previous match
2741 : //
2742 : //Channel matching for the actual spectral window of buffer
2743 : //
2744 0 : if(doConversion_p[vb.spectralWindow()])
2745 0 : matchChannel(vb.spectralWindow(), vb);
2746 : else
2747 : {
2748 0 : chanMap.resize();
2749 0 : chanMap=multiChanMap_p[vb.spectralWindow()];
2750 : }
2751 :
2752 0 : Vector<Int> rowFlags(vb.nRow());
2753 0 : rowFlags=0;
2754 0 : rowFlags(vb.flagRow())=true;
2755 0 : if(!usezero_p)
2756 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
2757 0 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
2758 :
2759 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
2760 0 : if (vb.antenna1()(rownr) != vb.antenna2()(rownr))
2761 0 : rowFlags(rownr) = (vb.flagRow()(rownr)==true);
2762 :
2763 0 : IPosition s,gradS;
2764 0 : Cube<Complex> visdata,gradVisAzData,gradVisElData;
2765 0 : if (whichVBColumn == FTMachine::MODEL)
2766 : {
2767 0 : s = vb.modelVisCube().shape();
2768 0 : visdata.reference(vb.modelVisCube());
2769 : }
2770 0 : else if (whichVBColumn == FTMachine::OBSERVED)
2771 : {
2772 0 : s = vb.visCube().shape();
2773 0 : visdata.reference(vb.visCube());
2774 : }
2775 :
2776 0 : if (doGrad)
2777 : {
2778 0 : if (whichGradVBColumn == FTMachine::MODEL)
2779 : {
2780 : // gradS = gradVBAz.modelVisCube().shape();
2781 0 : gradVisAzData.reference(gradVBAz.modelVisCube());
2782 0 : gradVisElData.reference(gradVBEl.modelVisCube());
2783 : }
2784 0 : else if (whichGradVBColumn == FTMachine::OBSERVED)
2785 : {
2786 : // gradS = gradVBAz.visCube().shape();
2787 0 : gradVisAzData.reference(gradVBAz.visCube());
2788 0 : gradVisElData.reference(gradVBEl.visCube());
2789 : }
2790 : }
2791 0 : visResampler_p->setParams(uvScale,uvOffset,dphase);
2792 0 : visResampler_p->setFreqMaps(expandedSpwFreqSel_p,expandedSpwConjFreqSel_p);
2793 0 : visResampler_p->setMaps(chanMap, polMap);
2794 : // Vector<Int> ConjCFMap, CFMap;
2795 0 : makeCFPolMap(vb,cfStokes_p,CFMap_p);
2796 0 : makeConjPolMap(vb,CFMap_p,ConjCFMap_p);
2797 0 : visResampler_p->setCFMaps(CFMap_p, ConjCFMap_p);
2798 :
2799 0 : if(isTiled)
2800 : {
2801 0 : Double invLambdaC=vb.frequency()(0)/C::c;
2802 0 : Vector<Double> uvLambda(2);
2803 0 : Vector<Int> centerLoc2D(2);
2804 0 : centerLoc2D=0;
2805 :
2806 : // Loop over all rows
2807 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
2808 : {
2809 :
2810 : // Calculate uvw for this row at the center frequency
2811 0 : uvLambda(0)=uvw(0, rownr)*invLambdaC;
2812 0 : uvLambda(1)=uvw(1, rownr)*invLambdaC;
2813 0 : centerLoc2D=gridder->location(centerLoc2D, uvLambda);
2814 :
2815 : // Is this point on the grid?
2816 0 : if(gridder->onGrid(centerLoc2D))
2817 : {
2818 :
2819 : // Get the tile
2820 0 : Array<Complex>* dataPtr=getDataPointer(centerLoc2D, true);
2821 0 : gridder->setOffset(IPosition(2, offsetLoc(0), offsetLoc(1)));
2822 0 : Int aNx=dataPtr->shape()(0);
2823 0 : Int aNy=dataPtr->shape()(1);
2824 :
2825 : // Now use FORTRAN to do the gridding. Remember to
2826 : // ensure that the shape and offsets of the tile are
2827 : // accounted for.
2828 :
2829 0 : Vector<Double> actualOffset(3);
2830 0 : for (Int i=0;i<2;i++)
2831 0 : actualOffset(i)=uvOffset(i)-Double(offsetLoc(i));
2832 :
2833 0 : actualOffset(2)=uvOffset(2);
2834 0 : IPosition s(vb.modelVisCube().shape());
2835 :
2836 0 : Int ScanNo=0, tmpPAI;
2837 0 : Double area=1.0;
2838 0 : tmpPAI = 1;
2839 0 : runFortranGetGrad(uvw,dphase,visdata,s,
2840 : gradVisAzData,gradVisElData,
2841 : Conj,flags,rowFlags,rownr,
2842 0 : actualOffset,dataPtr,aNx,aNy,npol,nchan,vb,NAnt,ScanNo,sigma,
2843 0 : l_offsets,m_offsets,area,doGrad,tmpPAI);
2844 0 : }
2845 : }
2846 0 : }
2847 : else
2848 : {
2849 :
2850 0 : IPosition s(vb.modelVisCube().shape());
2851 0 : Int ScanNo=0, tmpPAI;
2852 0 : Double area=1.0;
2853 :
2854 0 : tmpPAI = 1;
2855 :
2856 0 : runFortranGetGrad(uvw,dphase,visdata/*vb.modelVisCube()*/,s,
2857 : gradVisAzData, gradVisElData,
2858 : Conj,flags,rowFlags,row,
2859 0 : uvOffset,&griddedData,nx,ny,npol,nchan,vb,NAnt,ScanNo,sigma,
2860 0 : l_offsets,m_offsets,area,doGrad,tmpPAI);
2861 : // runFortranGet(uvw,dphase,vb.modelVisCube(),s,Conj,flags,rowFlags,row,
2862 : // uvOffset,&griddedData,nx,ny,npol,nchan,vb,NAnt,ScanNo,sigma,
2863 : // l_offsets,m_offsets,area,doGrad,tmpPAI);
2864 0 : }
2865 0 : }
2866 : //
2867 : //---------------------------------------------------------------
2868 : //
2869 0 : void AWProjectFT::get(VisBuffer& vb, Cube<Complex>& modelVis,
2870 : Array<Complex>& griddedVis, Vector<Double>& scale,
2871 : Int row)
2872 : {
2873 :
2874 : (void)scale; //Suppress the warning
2875 :
2876 0 : Int nX=griddedVis.shape()(0);
2877 0 : Int nY=griddedVis.shape()(1);
2878 0 : Vector<Double> offset(2);
2879 0 : offset(0)=Double(nX)/2.0;
2880 0 : offset(1)=Double(nY)/2.0;
2881 : // If row is -1 then we pass through all rows
2882 : Int startRow, endRow, nRow;
2883 0 : if (row==-1)
2884 : {
2885 0 : nRow=vb.nRow();
2886 0 : startRow=0;
2887 0 : endRow=nRow-1;
2888 0 : modelVis.set(Complex(0.0,0.0));
2889 : }
2890 : else
2891 : {
2892 0 : nRow=1;
2893 0 : startRow=row;
2894 0 : endRow=row;
2895 0 : modelVis.xyPlane(row)=Complex(0.0,0.0);
2896 : }
2897 :
2898 0 : Int NAnt=0;
2899 :
2900 0 : if (doPointing)
2901 0 : NAnt = findPointingOffsets(vb,l_offsets, m_offsets,true);
2902 :
2903 :
2904 : //
2905 : // Get the uvws in a form that Fortran can use
2906 : //
2907 0 : Matrix<Double> uvw(3, vb.uvw().nelements());
2908 0 : uvw=0.0;
2909 0 : Vector<Double> dphase(vb.uvw().nelements());
2910 0 : dphase=0.0;
2911 : //
2912 : //NEGATING to correct for an image inversion problem
2913 : //
2914 0 : for (Int i=startRow;i<=endRow;i++)
2915 : {
2916 0 : for (Int idim=0;idim<2;idim++) uvw(idim,i)=-vb.uvw()(i)(idim);
2917 0 : uvw(2,i)=vb.uvw()(i)(2);
2918 : }
2919 :
2920 0 : doUVWRotation_p=true;
2921 : //rotateUVW(uvw, dphase, vb);
2922 0 : girarUVW(uvw, dphase, vb);
2923 0 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
2924 :
2925 : // This is the convention for dphase
2926 : //dphase*=-1.0;
2927 :
2928 0 : Cube<Int> flags(vb.flagCube().shape());
2929 0 : flags=0;
2930 0 : flags(vb.flagCube())=true;
2931 :
2932 : //Check if ms has changed then cache new spw and chan selection
2933 0 : if(vb.newMS())
2934 0 : matchAllSpwChans(vb);
2935 :
2936 : //Channel matching for the actual spectral window of buffer
2937 0 : if(doConversion_p[vb.spectralWindow()])
2938 0 : matchChannel(vb.spectralWindow(), vb);
2939 : else
2940 : {
2941 0 : chanMap.resize();
2942 0 : chanMap=multiChanMap_p[vb.spectralWindow()];
2943 : }
2944 :
2945 0 : Vector<Int> rowFlags(vb.nRow());
2946 0 : rowFlags=0;
2947 0 : rowFlags(vb.flagRow())=true;
2948 0 : if(!usezero_p)
2949 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
2950 0 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
2951 :
2952 0 : visResampler_p->setParams(uvScale,uvOffset,dphase);
2953 0 : visResampler_p->setMaps(chanMap, polMap);
2954 0 : visResampler_p->setFreqMaps(expandedSpwFreqSel_p,expandedSpwConjFreqSel_p);
2955 :
2956 0 : IPosition s(modelVis.shape());
2957 0 : Int Conj=0,doGrad=0,ScanNo=0;
2958 0 : Double area=1.0;
2959 0 : Int tmpPAI=1;
2960 0 : runFortranGet(uvw,dphase,vb.modelVisCube(),s,Conj,flags,rowFlags,row,
2961 0 : offset,&griddedVis,nx,ny,npol,nchan,vb,NAnt,ScanNo,sigma,
2962 0 : l_offsets,m_offsets,area,doGrad,tmpPAI);
2963 0 : }
2964 : //
2965 : //----------------------------------------------------------------------
2966 : //
2967 0 : void AWProjectFT::initVisBuffer(VisBuffer& vb, Type whichVBColumn)
2968 : {
2969 0 : if (whichVBColumn == FTMachine::MODEL) vb.modelVisCube()=Complex(0.0,0.0);
2970 0 : else if (whichVBColumn == FTMachine::OBSERVED) vb.visCube()=Complex(0.0,0.0);
2971 0 : }
2972 : //
2973 : //----------------------------------------------------------------------
2974 : //
2975 0 : void AWProjectFT::initVisBuffer(VisBuffer& vb, Type whichVBColumn, Int row)
2976 : {
2977 0 : if (whichVBColumn == FTMachine::MODEL)
2978 0 : vb.modelVisCube().xyPlane(row)=Complex(0.0,0.0);
2979 0 : else if (whichVBColumn == FTMachine::OBSERVED)
2980 0 : vb.visCube().xyPlane(row)=Complex(0.0,0.0);
2981 0 : }
2982 :
2983 0 : void AWProjectFT::ComputeResiduals(VisBuffer&vb, Bool useCorrected)
2984 : {
2985 0 : VBStore vbs;
2986 :
2987 0 : vbs.nRow_p = vb.nRow();
2988 0 : vbs.beginRow_p = 0;
2989 0 : vbs.endRow_p = vbs.nRow_p;
2990 :
2991 0 : vbs.modelCube_p.reference(vb.modelVisCube());
2992 0 : if (useCorrected) vbs.correctedCube_p.reference(vb.correctedVisCube());
2993 0 : else vbs.visCube_p.reference(vb.visCube());
2994 0 : vbs.useCorrected_p = useCorrected;
2995 0 : visResampler_p->ComputeResiduals(vbs);
2996 0 : }
2997 :
2998 : // void AWProjectFT::ComputeResiduals(VBStore &vbs)
2999 : // {
3000 : // vbs.nRow_p = vb.nRow();
3001 : // vbs.modelCube_p.reference(vb.modelVisCube());
3002 : // if (useCorrected) vbs.correctedCube_p.reference(vb.correctedVisCube());
3003 : // else vbs.visCube_p.reference(vb.visCube());
3004 : // vbs.useCorrected_p = useCorrected;
3005 : // visResampler_p->ComputeResiduals(vbs);
3006 : // }
3007 :
3008 : } //# NAMESPACE CASA - END
3009 :
3010 :
3011 :
3012 :
3013 : // //
3014 : // //---------------------------------------------------------------
3015 : // //
3016 : // void AWProjectFT::normalizeImage(Lattice<Complex>& skyImage,
3017 : // const Matrix<Double>& sumOfWts,
3018 : // Lattice<Float>& sensitivityImage,
3019 : // Lattice<Complex>& sensitivitySqImage,
3020 : // Bool fftNorm)
3021 : // {
3022 : // //
3023 : // // Apply the gridding correction
3024 : // //
3025 : // Int inx = skyImage.shape()(0);
3026 : // Int iny = skyImage.shape()(1);
3027 : // // Vector<Complex> correction(inx);
3028 :
3029 : // Vector<Float> sincConv(nx);
3030 : // Float centerX=nx/2;
3031 : // for (Int ix=0;ix<nx;ix++)
3032 : // {
3033 : // Float x=C::pi*Float(ix-centerX)/(Float(nx)*Float(convSampling));
3034 : // if(ix==centerX) sincConv(ix)=1.0;
3035 : // else sincConv(ix)=sin(x)/x;
3036 : // }
3037 :
3038 : // IPosition cursorShape(4, inx, 1, 1, 1);
3039 : // IPosition axisPath(4, 0, 1, 2, 3);
3040 : // LatticeStepper lsx(skyImage.shape(), cursorShape, axisPath);
3041 : // LatticeIterator<Complex> lix(skyImage, lsx);
3042 :
3043 : // LatticeStepper lavgpb(sensitivityImage.shape(),cursorShape,axisPath);
3044 : // LatticeIterator<Float> liavgpb(sensitivityImage, lavgpb);
3045 : // LatticeStepper lavgpbSq(sensitivitySqImage.shape(),cursorShape,axisPath);
3046 : // LatticeIterator<Complex> liavgpbSq(sensitivitySqImage, lavgpbSq);
3047 :
3048 : // for(lix.reset(),liavgpb.reset(),liavgpbSq.reset();
3049 : // !lix.atEnd();
3050 : // lix++,liavgpb++,liavgpbSq++)
3051 : // {
3052 : // Int pol=lix.position()(2);
3053 : // Int chan=lix.position()(3);
3054 : // // Int iy=lix.position()(1);
3055 :
3056 : // // gridder->correctX1D(correction,iy);
3057 :
3058 : // Vector<Complex> PBCorrection(liavgpb.rwVectorCursor().shape());
3059 : // Vector<Float> avgPBVec(liavgpb.rwVectorCursor().shape());
3060 : // Vector<Complex> avgPBSqVec(liavgpbSq.rwVectorCursor().shape());
3061 :
3062 : // avgPBSqVec= liavgpbSq.rwVectorCursor();
3063 : // avgPBVec = liavgpb.rwVectorCursor();
3064 :
3065 : // for(int i=0;i<PBCorrection.shape();i++)
3066 : // {
3067 : // PBCorrection(i)=(avgPBSqVec(i)/avgPBVec(i));///(sincConv(i)*sincConv(iy));
3068 : // if ((abs(avgPBVec(i))) >= pbLimit_p)
3069 : // lix.rwVectorCursor()(i) /= PBCorrection(i);
3070 : // }
3071 :
3072 : // // if(sumOfWts(pol, chan)>0.0) // UUU Changing this again....
3073 : // {
3074 : // if(fftNorm)
3075 : // {
3076 : // if(sumOfWts(pol, chan)>0.0) // UUU Moving this inside the fftNorm check.
3077 : // {
3078 : // Complex rnorm(Float(inx)*Float(iny)/sumOfWts(pol,chan));
3079 : // lix.rwCursor()*=rnorm;
3080 : // }
3081 : // else
3082 : // lix.woCursor()=0.0;
3083 : // }
3084 : // else
3085 : // {
3086 : // Complex rnorm(Float(inx)*Float(iny));
3087 : // lix.rwCursor()*=rnorm;
3088 : // }
3089 : // }
3090 : // //else
3091 : // //lix.woCursor()=0.0;
3092 : // }
3093 : // }
3094 : // //
3095 : // //---------------------------------------------------------------
3096 : // //
3097 : // void AWProjectFT::normalizeImage(Lattice<Complex>& skyImage,
3098 : // const Matrix<Double>& sumOfWts,
3099 : // Lattice<Float>& sensitivityImage,
3100 : // Bool fftNorm)
3101 : // {
3102 : // //
3103 : // // Apply the gridding correction
3104 : // //
3105 : // Int inx = skyImage.shape()(0);
3106 : // Int iny = skyImage.shape()(1);
3107 : // Vector<Complex> correction(inx);
3108 :
3109 : // Vector<Float> sincConv(nx);
3110 : // Float centerX=nx/2;
3111 : // for (Int ix=0;ix<nx;ix++)
3112 : // {
3113 : // Float x=C::pi*Float(ix-centerX)/(Float(nx)*Float(convSampling));
3114 : // if(ix==centerX) sincConv(ix)=1.0;
3115 : // else sincConv(ix)=sin(x)/x;
3116 : // }
3117 :
3118 : // IPosition cursorShape(4, inx, 1, 1, 1);
3119 : // IPosition axisPath(4, 0, 1, 2, 3);
3120 : // LatticeStepper lsx(skyImage.shape(), cursorShape, axisPath);
3121 : // // LatticeIterator<Complex> lix(skyImage, lsx);
3122 : // LatticeIterator<Complex> lix(skyImage, lsx);
3123 :
3124 : // LatticeStepper lavgpb(sensitivityImage.shape(),cursorShape,axisPath);
3125 : // // Array<Float> senArray;sensitivityImage.get(senArray,true);
3126 : // // ArrayLattice<Float> senLat(senArray,true);
3127 : // // LatticeIterator<Float> liavgpb(senLat, lavgpb);
3128 : // LatticeIterator<Float> liavgpb(sensitivityImage, lavgpb);
3129 :
3130 : // for(lix.reset(),liavgpb.reset();
3131 : // !lix.atEnd();
3132 : // lix++,liavgpb++)
3133 : // {
3134 : // Int pol=lix.position()(2);
3135 : // Int chan=lix.position()(3);
3136 :
3137 : // // if(sumOfWts(pol, chan)>0.0) // UUU Changing this again.....
3138 : // {
3139 : // Int iy=lix.position()(1);
3140 : // gridder->correctX1D(correction,iy);
3141 :
3142 : // Vector<Float> avgPBVec(liavgpb.rwVectorCursor().shape());
3143 :
3144 : // avgPBVec = liavgpb.rwVectorCursor();
3145 :
3146 : // for(int i=0;i<avgPBVec.shape();i++)
3147 : // {
3148 : // //
3149 : // // This with the PS functions
3150 : // //
3151 : // // PBCorrection(i)=FUNC(avgPBVec(i))*sincConv(i)*sincConv(iy);
3152 : // // if ((abs(PBCorrection(i)*correction(i))) >= pbLimit_p)
3153 : // // lix.rwVectorCursor()(i) /= PBCorrection(i)*correction(i);
3154 : // // else if (!makingPSF)
3155 : // // lix.rwVectorCursor()(i) /= correction(i)*sincConv(i)*sincConv(iy);
3156 : // //
3157 : // // This without the PS functions
3158 : // //
3159 : // Float tt = pbFunc(avgPBVec(i),pbLimit_p)*sincConv(i)*sincConv(iy);
3160 : // // PBCorrection(i)=pbFunc(avgPBVec(i),pbLimit_p)*sincConv(i)*sincConv(iy);
3161 : // // lix.rwVectorCursor()(i) /= PBCorrection(i);
3162 : // // lix.rwVectorCursor()(i) *= tt;
3163 : // // if (!makingPSF)
3164 :
3165 : // lix.rwVectorCursor()(i) /= tt;
3166 : // }
3167 :
3168 : // if(fftNorm)
3169 : // {
3170 : // if(sumOfWts(pol, chan)>0.0) // UUU Moved this inside the fftNorm check.
3171 : // {
3172 : // Complex rnorm(Float(inx)*Float(iny)/sumOfWts(pol,chan));
3173 : // lix.rwCursor()*=rnorm;
3174 : // }
3175 : // else
3176 : // lix.woCursor()=0.0;
3177 : // }
3178 : // else
3179 : // {
3180 : // Complex rnorm(Float(inx)*Float(iny));
3181 : // lix.rwCursor()*=rnorm;
3182 : // }
3183 : // }
3184 : // //else
3185 : // //lix.woCursor()=0.0;
3186 : // }
3187 : // }
|
