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/TransformMachines2/AWProjectFT.h>
44 : #include <synthesis/TransformMachines2/CFStore2.h>
45 : #include <synthesis/MeasurementComponents/ExpCache.h>
46 : #include <synthesis/MeasurementComponents/CExp.h>
47 : #include <synthesis/TransformMachines2/AWVisResampler.h>
48 : #include <synthesis/TransformMachines2/VBStore.h>
49 :
50 : #include <casacore/scimath/Mathematics/FFTServer.h>
51 : #include <casacore/scimath/Mathematics/MathFunc.h>
52 : #include <casacore/measures/Measures/MeasTable.h>
53 : #include <iostream>
54 : #include <casacore/casa/OS/Timer.h>
55 :
56 : #include <synthesis/TransformMachines2/ATerm.h>
57 : #include <synthesis/TransformMachines2/NoOpATerm.h>
58 : #include <synthesis/TransformMachines2/AWConvFunc.h>
59 : #include <synthesis/TransformMachines2/EVLAAperture.h>
60 : #include <iomanip>
61 : //#define CONVSIZE (1024*2)
62 : // #define OVERSAMPLING 2
63 : #define USETABLES 0 // If equal to 1, use tabulated exp() and
64 : // complex exp() functions.
65 : #define MAXPOINTINGERROR 250.0 // Max. pointing error in arcsec used to
66 : // determine the resolution of the
67 : // tabulated exp() function.
68 : #define DORES true
69 :
70 :
71 : using namespace casacore;
72 : namespace casa { //# NAMESPACE CASA - BEGIN
73 : using namespace vi;
74 : #define NEED_UNDERSCORES
75 : namespace refim{
76 : extern "C"
77 : {
78 : //
79 : // The Gridding Convolution Function (GCF) used by the underlying
80 : // gridder written in FORTRAN.
81 : //
82 : // The arguments must all be pointers and the value of the GCF at
83 : // the given (u,v) point is returned in the weight variable. Making
84 : // this a function which returns a complex value (namely the weight)
85 : // has problems when called in FORTRAN - I (SB) don't understand
86 : // why.
87 : //
88 : #if defined(NEED_UNDERSCORES)
89 : #define nwcppeij nwcppeij_
90 : #endif
91 : //
92 : //---------------------------------------------------------------
93 : //
94 : IlluminationConvFunc awEij2;
95 0 : void awcppeij2(Double *griduvw, Double *area,
96 : Double *raoff1, Double *decoff1,
97 : Double *raoff2, Double *decoff2,
98 : Int *doGrad,
99 : Complex *weight,
100 : Complex *dweight1,
101 : Complex *dweight2,
102 : Double *currentCFPA)
103 : {
104 0 : Complex w,d1,d2;
105 0 : awEij2.getValue(griduvw, raoff1, raoff2, decoff1, decoff2,
106 : area,doGrad,w,d1,d2,*currentCFPA);
107 0 : *weight = w;
108 0 : *dweight1 = d1;
109 0 : *dweight2 = d2;
110 0 : }
111 : }
112 :
113 135 : ATerm* AWProjectFT::createTelescopeATerm(const String& telescopeName, const Bool& isATermOn)
114 : {
115 :
116 135 : if (!isATermOn) return new NoOpATerm();
117 :
118 : // MSObservationColumns msoc(ms.observation());
119 : // String ObsName=msoc.telescopeName()(0);
120 134 : if ((telescopeName == "EVLA") || (telescopeName == "VLA"))
121 134 : return new EVLAAperture();
122 : else
123 : {
124 0 : LogIO os(LogOrigin("AWProjectFT2", "createTelescopeATerm",WHERE));
125 0 : os << "Telescope name ('"+
126 0 : telescopeName+"') in the MS not recognized to create the telescope specific ATerm"
127 0 : << LogIO::EXCEPTION;
128 0 : }
129 :
130 0 : return NULL;
131 : }
132 :
133 135 : CountedPtr<ConvolutionFunction> AWProjectFT::makeCFObject(const String& telescopeName,
134 : const Bool aTermOn,
135 : const Bool psTermOn,
136 : const Bool wTermOn,
137 : const Bool,// mTermOn,
138 : const Bool wbAWP,
139 : const Bool conjBeams)
140 : {
141 : (void)wTermOn;//Unused
142 135 : CountedPtr<ATerm> apertureFunction = AWProjectFT::createTelescopeATerm(telescopeName, aTermOn);
143 135 : CountedPtr<PSTerm> psTerm = new PSTerm();
144 135 : CountedPtr<WTerm> wTerm = new WTerm();
145 : //if (cfBufferSize > 0) apertureFunction->setConvSize(cfBufferSize);
146 : //
147 : // Selectively switch off CFTerms.
148 : //
149 135 : if (aTermOn == false) {apertureFunction->setOpCode(CFTerms::NOOP);}
150 135 : if (psTermOn == false) psTerm->setOpCode(CFTerms::NOOP);
151 135 : if (wTermOn == False) wTerm->setOpCode(CFTerms::NOOP);
152 : //
153 : // Construct the CF object with appropriate CFTerms.
154 : //
155 135 : CountedPtr<ConvolutionFunction> awConvFunc;
156 : // awConvFunc = new AWConvFunc(apertureFunction,psTerm,wTerm, !wbAWP);
157 : //if ((ftmName=="mawprojectft") || (mTermOn))
158 :
159 135 : awConvFunc = new AWConvFunc(apertureFunction,psTerm,wTerm,wbAWP, conjBeams);
160 :
161 270 : return awConvFunc;
162 135 : }
163 : //
164 : //---------------------------------------------------------------
165 : //
166 0 : AWProjectFT::AWProjectFT()
167 0 : : FTMachine(), padding_p(1.0), nWPlanes_p(1),
168 0 : imageCache(0), cachesize(0), tilesize(16),
169 0 : gridder(0), isTiled(false), lattice( ),
170 0 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
171 0 : pointingToImage(0), usezero_p(false), avgPB_p(nullptr),
172 0 : epJ_p(nullptr),
173 0 : doPBCorrection(true), conjBeams_p(true),/*cfCache_p(cfcache),*/ paChangeDetector(),
174 0 : rotateOTFPAIncr_p(0.1),
175 0 : Second("s"),Radian("rad"),Day("d"), pbNormalized_p(false), paNdxProcessed_p(),
176 0 : visResampler_p(nullptr), sensitivityPatternQualifier_p(-1),sensitivityPatternQualifierStr_p(""),
177 0 : rotatedConvFunc_p(),
178 0 : runTime1_p(0.0), previousSPWID_p(-1), self_p(nullptr), vb2CFBMap_p(nullptr), po_p(nullptr),wbAWP_p(true), timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0)
179 : {
180 : // convSize=0;
181 0 : tangentSpecified_p=false;
182 0 : lastIndex_p=0;
183 0 : paChangeDetector.reset();
184 0 : pbLimit_p=5e-2;
185 : //
186 : // Get various parameters from the visibilities.
187 : //
188 0 : doPointing=1;
189 :
190 0 : maxConvSupport=-1;
191 : //
192 : // Set up the Conv. Func. disk cache manager object.
193 : //
194 : // if (!cfCache_p.null()) delete &cfCache_p;
195 : // cfCache_p=cfcache;
196 0 : convSampling=-1;
197 : //convSize=CONVSIZE;
198 0 : Long hostRAM = (HostInfo::memoryTotal(true)*1024); // In bytes
199 0 : hostRAM = hostRAM/(sizeof(Float)*2); // In complex pixels
200 0 : if (cachesize > hostRAM) cachesize=hostRAM;
201 0 : sigma=1.0;
202 0 : canComputeResiduals_p=DORES;
203 : // cfs2_p = &cfCache_p->memCache2_p[0];//new CFStore2;
204 : // cfwts2_p = &cfCache_p->memCacheWt2_p[0];//new CFStore2;
205 0 : pop_p->init();
206 0 : CFBuffer::initCFBStruct(cfbst_pub);
207 : // rotatedConvFunc_p.data=new Array<Complex>();
208 : // self_p.reset(this);
209 0 : vb2CFBMap_p = new VB2CFBMap();
210 0 : po_p = new PointingOffsets();
211 0 : po_p->setOverSampling(convSampling);
212 0 : }
213 : //
214 : //---------------------------------------------------------------
215 : //
216 93 : AWProjectFT::AWProjectFT(Int nWPlanes, Long icachesize,
217 : CountedPtr<CFCache>& cfcache,
218 : CountedPtr<ConvolutionFunction>& cf,
219 : CountedPtr<VisibilityResamplerBase>& visResampler,
220 : Bool applyPointingOffset,
221 : vector<float> pointingOffsetSigDev,
222 : Bool doPBCorr,
223 : Int itilesize,
224 : Float pbLimit,
225 : Bool usezero,
226 : Bool conjBeams,
227 : Bool doublePrecGrid,
228 93 : PolOuterProduct::MuellerType muellerType)
229 93 : : FTMachine(cfcache,cf), padding_p(1.0), nWPlanes_p(nWPlanes),
230 93 : imageCache(0), cachesize(icachesize), tilesize(itilesize),
231 93 : gridder(0), isTiled(false), lattice( ),
232 93 : maxAbsData(0.0), centerLoc(IPosition(4,0)), offsetLoc(IPosition(4,0)),
233 93 : pointingToImage(0), usezero_p(usezero), avgPB_p(nullptr),
234 : // convFunc_p(), convWeights_p(),
235 93 : epJ_p(),
236 93 : doPBCorrection(doPBCorr), conjBeams_p(conjBeams),
237 93 : /*cfCache_p(cfcache),*/ paChangeDetector(),
238 93 : rotateOTFPAIncr_p(0.1),
239 93 : Second("s"),Radian("rad"),Day("d"), pbNormalized_p(false),
240 93 : visResampler_p(visResampler), sensitivityPatternQualifier_p(-1),sensitivityPatternQualifierStr_p(""),
241 279 : rotatedConvFunc_p(), runTime1_p(0.0), previousSPWID_p(-1),self_p(nullptr), vb2CFBMap_p(nullptr), po_p(nullptr),wbAWP_p(true), timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0)
242 : {
243 : //convSize=0;
244 93 : tangentSpecified_p=false;
245 93 : lastIndex_p=0;
246 93 : paChangeDetector.reset();
247 93 : pbLimit_p=pbLimit;
248 : //
249 : // Get various parameters from the visibilities.
250 : //
251 93 : if (applyPointingOffset) doPointing=1; else doPointing=0;
252 :
253 93 : maxConvSupport=-1;
254 : //
255 : // Set up the Conv. Func. disk cache manager object.
256 : //
257 : // if (!cfCache_p.null()) delete &cfCache_p;
258 : // cfCache_p=cfcache;
259 93 : convSampling=convFuncCtor_p->getOversampling();
260 : //convSize=CONVSIZE;
261 93 : Long hostRAM = (HostInfo::memoryTotal(true)*1024); // In bytes
262 93 : hostRAM = hostRAM/(sizeof(Float)*2); // In complex pixels
263 93 : if (cachesize > hostRAM) cachesize=hostRAM;
264 93 : sigma=1.0;
265 93 : canComputeResiduals_p=DORES;
266 93 : if (!cfCache_p.null())
267 : {
268 0 : cfs2_p = CountedPtr<CFStore2>(&(cfCache_p->memCache2_p)[0],false);//new CFStore2;
269 0 : cfwts2_p = CountedPtr<CFStore2>(&cfCache_p->memCacheWt2_p[0],false);//new CFStore2;
270 : }
271 93 : pop_p->init();
272 93 : useDoubleGrid_p=doublePrecGrid;
273 : // rotatedConvFunc_p.data=new Array<Complex>();
274 93 : CFBuffer::initCFBStruct(cfbst_pub);
275 93 : muellerType_p = muellerType;
276 : // self_p.reset(this);
277 93 : vb2CFBMap_p = new VB2CFBMap();
278 93 : po_p = new PointingOffsets();
279 93 : po_p->setOverSampling(convSampling);
280 93 : vb2CFBMap_p->setPOSigmaDev(pointingOffsetSigDev);
281 93 : wbAWP_p=convFuncCtor_p->isWBAWP();
282 :
283 93 : }
284 : //
285 : //---------------------------------------------------------------
286 : //
287 0 : AWProjectFT::AWProjectFT(const RecordInterface& stateRec)
288 0 : : FTMachine(),Second("s"),Radian("rad"),Day("d"),visResampler_p(nullptr), self_p(nullptr), vb2CFBMap_p(nullptr), po_p(nullptr),wbAWP_p(true), timemass_p(0.0), timegrid_p(0.0), timedegrid_p(0.0)
289 : {
290 : //
291 : // Construct from the input state record
292 : //
293 0 : String error;
294 :
295 0 : if (!fromRecord(stateRec)) {
296 0 : LogIO log_l(LogOrigin("AWProjectFT2", "AWProjectFT[R&D]"));
297 0 : log_l << "Failed to create " << name() << " object." << LogIO::EXCEPTION;
298 0 : };
299 0 : maxConvSupport=-1;
300 : //convSampling=OVERSAMPLING;
301 0 : convSampling=-1;
302 0 : visResampler_p->init(useDoubleGrid_p);
303 : //convSize=CONVSIZE;
304 0 : canComputeResiduals_p=DORES;
305 0 : if (!cfCache_p.null())
306 : {
307 0 : cfs2_p = CountedPtr<CFStore2>(&cfCache_p->memCache2_p[0],false);//new CFStore2;
308 0 : cfwts2_p = CountedPtr<CFStore2>(&cfCache_p->memCacheWt2_p[0],false);//new CFStore2;
309 : }
310 0 : pop_p->init();
311 : // self_p.reset(this);
312 0 : vb2CFBMap_p = new VB2CFBMap();
313 0 : po_p = new PointingOffsets();
314 0 : po_p->setOverSampling(convSampling);
315 0 : }
316 : //
317 : //----------------------------------------------------------------------
318 : //
319 159 : AWProjectFT::AWProjectFT(const AWProjectFT& other):FTMachine()
320 : {
321 159 : operator=(other);
322 159 : }
323 : //
324 : //---------------------------------------------------------------
325 : //
326 : // This is nasty, we should use CountedPointers here.
327 252 : AWProjectFT::~AWProjectFT()
328 : {
329 252 : if(imageCache) delete imageCache;
330 252 : imageCache=0;
331 252 : if(gridder) delete gridder;
332 252 : gridder=0;
333 252 : }
334 : //
335 : //---------------------------------------------------------------
336 : //
337 318 : AWProjectFT& AWProjectFT::operator=(const AWProjectFT& other)
338 : {
339 318 : if(this!=&other)
340 : {
341 : //Do the base parameters
342 318 : FTMachine::operator=(other);
343 :
344 :
345 318 : padding_p=other.padding_p;
346 :
347 318 : nWPlanes_p=other.nWPlanes_p;
348 318 : imageCache=other.imageCache;
349 318 : cachesize=other.cachesize;
350 318 : tilesize=other.tilesize;
351 318 : cfRefFreq_p = other.cfRefFreq_p;
352 318 : if(other.gridder==0)
353 318 : gridder=0;
354 : else
355 : {
356 0 : uvScale.resize();
357 0 : uvOffset.resize();
358 0 : uvScale=other.uvScale;
359 0 : uvOffset=other.uvOffset;
360 0 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
361 0 : uvScale, uvOffset,
362 0 : "SF");
363 : }
364 :
365 318 : isTiled=other.isTiled;
366 318 : lattice=0;
367 :
368 318 : maxAbsData=other.maxAbsData;
369 318 : centerLoc=other.centerLoc;
370 318 : offsetLoc=other.offsetLoc;
371 318 : pointingToImage=other.pointingToImage;
372 318 : usezero_p=other.usezero_p;
373 :
374 :
375 318 : padding_p=other.padding_p;
376 318 : imageCache=other.imageCache;
377 318 : cachesize=other.cachesize;
378 318 : tilesize=other.tilesize;
379 318 : isTiled=other.isTiled;
380 318 : maxAbsData=other.maxAbsData;
381 318 : centerLoc=other.centerLoc;
382 318 : offsetLoc=other.offsetLoc;
383 318 : pointingToImage=other.pointingToImage;
384 318 : useDoubleGrid_p=other.useDoubleGrid_p;
385 318 : doPBCorrection = other.doPBCorrection;
386 318 : maxConvSupport= other.maxConvSupport;
387 :
388 318 : epJ_p=other.epJ_p;
389 : //convSize=other.convSize;
390 318 : lastIndex_p=other.lastIndex_p;
391 318 : paChangeDetector=other.paChangeDetector;
392 318 : pbLimit_p=other.pbLimit_p;
393 : //
394 : // Get various parameters from the visibilities.
395 : //
396 318 : doPointing=other.doPointing;
397 :
398 318 : maxConvSupport=other.maxConvSupport;
399 : //
400 : // Set up the Conv. Func. disk cache manager object.
401 : //
402 318 : cfCache_p=other.cfCache_p;
403 318 : convSampling=other.convSampling;
404 : //convSize=other.convSize;
405 318 : cachesize=other.cachesize;
406 :
407 318 : currentCFPA=other.currentCFPA;
408 318 : lastPAUsedForWtImg = other.lastPAUsedForWtImg;
409 318 : avgPB_p = other.avgPB_p;
410 :
411 318 : convFuncCtor_p = other.convFuncCtor_p;
412 318 : pbNormalized_p = other.pbNormalized_p;
413 318 : sensitivityPatternQualifier_p = other.sensitivityPatternQualifier_p;
414 318 : sensitivityPatternQualifierStr_p = other.sensitivityPatternQualifierStr_p;
415 318 : visResampler_p=other.visResampler_p; // Copy the counted pointer
416 : // visResampler_p=other.visResampler_p->clone();
417 : // *visResampler_p = *other.visResampler_p; // Call the appropriate operator=()
418 :
419 318 : rotatedConvFunc_p = other.rotatedConvFunc_p;
420 318 : cfs2_p = other.cfs2_p;
421 318 : cfwts2_p = other.cfwts2_p;
422 318 : paNdxProcessed_p = other.paNdxProcessed_p;
423 318 : imRefFreq_p = other.imRefFreq_p;
424 318 : conjBeams_p = other.conjBeams_p;
425 318 : rotateOTFPAIncr_p=other.rotateOTFPAIncr_p;
426 318 : computePAIncr_p=other.computePAIncr_p;
427 318 : runTime1_p = other.runTime1_p;
428 318 : muellerType_p = other.muellerType_p;
429 318 : previousSPWID_p = other.previousSPWID_p;
430 318 : vb2CFBMap_p = other.vb2CFBMap_p;
431 318 : po_p = other.po_p;
432 : // self_p = other.self_p;
433 318 : wbAWP_p=other.wbAWP_p;
434 318 : timemass_p=0.0;
435 318 : timegrid_p=0.0;
436 318 : timedegrid_p=0.0;
437 : };
438 318 : return *this;
439 : };
440 : //
441 : //----------------------------------------------------------------------
442 : //
443 211 : void AWProjectFT::init()
444 : {
445 422 : LogIO log_l(LogOrigin("AWProjectFT2", "init[R&D]"));
446 :
447 211 : nx = image->shape()(0);
448 211 : ny = image->shape()(1);
449 211 : npol = image->shape()(2);
450 211 : nchan = image->shape()(3);
451 :
452 :
453 211 : sumWeight.resize(npol, nchan);
454 211 : sumCFWeight.resize(npol, nchan);
455 :
456 211 : wConvSize=max(1, nWPlanes_p);
457 :
458 211 : CoordinateSystem cs=image->coordinates();
459 211 : uvScale.resize(3);
460 211 : uvScale=0.0;
461 211 : uvScale(0)=Float(nx)*cs.increment()(0);
462 211 : uvScale(1)=Float(ny)*cs.increment()(1);
463 211 : uvScale(2)=Float(wConvSize)*abs(cs.increment()(0));
464 :
465 211 : Int index= cs.findCoordinate(Coordinate::SPECTRAL);
466 211 : SpectralCoordinate spCS = cs.spectralCoordinate(index);
467 211 : imRefFreq_p = spCS.referenceValue()(0);
468 :
469 211 : uvOffset.resize(3);
470 211 : uvOffset(0)=nx/2;
471 211 : uvOffset(1)=ny/2;
472 211 : uvOffset(2)=0;
473 :
474 211 : if(gridder) delete gridder;
475 211 : gridder=0;
476 422 : gridder = new ConvolveGridder<Double, Complex>(IPosition(2, nx, ny),
477 211 : uvScale, uvOffset,
478 211 : "SF");
479 :
480 : // Set up image cache needed for gridding.
481 211 : if(imageCache) delete imageCache;
482 211 : imageCache=0;
483 :
484 : // The tile size should be large enough that the
485 : // extended convolution function can fit easily
486 211 : if(isTiled)
487 : {
488 0 : Float tileOverlap=0.5;
489 0 : tilesize=min(256,tilesize);
490 0 : IPosition tileShape=IPosition(4,tilesize,tilesize,npol,nchan);
491 0 : Vector<Float> tileOverlapVec(4);
492 0 : tileOverlapVec=0.0;
493 0 : tileOverlapVec(0)=tileOverlap;
494 0 : tileOverlapVec(1)=tileOverlap;
495 : if (sizeof(long) < 4) // 32-bit machine
496 : {
497 : Int tmpCacheVal=static_cast<Int>(cachesize);
498 : imageCache=new LatticeCache <Complex> (*image, tmpCacheVal, tileShape,
499 : tileOverlapVec,
500 : (tileOverlap>0.0));
501 : }
502 : else // 64-bit machine
503 : {
504 0 : Long tmpCacheVal=cachesize;
505 0 : imageCache=new LatticeCache <Complex> (*image, tmpCacheVal, tileShape,
506 : tileOverlapVec,
507 0 : (tileOverlap>0.0));
508 : }
509 0 : }
510 :
511 : #if(USETABLES)
512 : Double StepSize;
513 : Int N=500000;
514 : StepSize = abs((((2*nx)/uvScale(0))/(sigma) +
515 : MAXPOINTINGERROR*1.745329E-02*(sigma)/3600.0))/N;
516 : if (!awEij2.isReady())
517 : {
518 : log_l << "Making lookup table for exp function with a resolution of "
519 : << StepSize << " radians. "
520 : << "Memory used: " << sizeof(Float)*N/(1024.0*1024.0)<< " MB."
521 : << LogIO::NORMAL
522 : <<LogIO::POST;
523 :
524 : awEij2.setSigma(sigma);
525 : awEij2.initExpTable(N,StepSize);
526 : // ExpTab.build(N,StepSize);
527 :
528 : log_l << "Making lookup table for complex exp function with a resolution of "
529 : << 2*M_PI/N << " radians. "
530 : << "Memory used: " << 2*sizeof(Float)*N/(1024.0*1024.0) << " MB."
531 : << LogIO::NORMAL
532 : << LogIO::POST;
533 : awEij2.initCExpTable(N);
534 : // CExpTab.build(N);
535 : }
536 : #endif
537 : // vpSJ->reset();
538 211 : paChangeDetector.reset();
539 211 : makingPSF = false;
540 :
541 211 : log_l << "Using " << visResampler_p->name() << " visibility resampler (a.k.a. \"gridder/degridder\")" << LogIO::POST;
542 211 : vb2CFBMap_p->computePhaseScreen_p = visResampler_p->needCFPhaseScreen();
543 : //cerr << "Current runTime = " << runTime << endl;
544 211 : }
545 : //
546 : //---------------------------------------------------------------
547 : //
548 0 : MDirection::Convert AWProjectFT::makeCoordinateMachine(const VisBuffer2& vb,
549 : const MDirection::Types& From,
550 : const MDirection::Types& To,
551 : MEpoch& last)
552 : {
553 0 : LogIO log_l(LogOrigin("AWProjectFT2","makeCoordinateMachine[R&D]"));
554 0 : Double time = getCurrentTimeStamp(vb);
555 :
556 0 : MEpoch epoch(Quantity(time,Second),MEpoch::TAI);
557 : // epoch = MEpoch(Quantity(time,Second),MEpoch::TAI);
558 : //
559 : // ...now make an object to hold the observatory position info...
560 : //
561 0 : MPosition pos;
562 0 : String ObsName=(vb.subtableColumns()).observation().telescopeName()(vb.arrayId()(0));
563 :
564 0 : if (!MeasTable::Observatory(pos,ObsName))
565 0 : log_l << "Observatory position for "+ ObsName + " not found"
566 0 : << LogIO::EXCEPTION;
567 : //
568 : // ...now make a Frame object out of the observatory position and
569 : // time objects...
570 : //
571 0 : MeasFrame frame(epoch,pos);
572 : //
573 : // ...finally make the convert machine.
574 : //
575 0 : MDirection::Convert mac(MDirection::Ref(From,frame),
576 0 : MDirection::Ref(To,frame));
577 :
578 0 : MEpoch::Convert toLAST = MEpoch::Convert(MEpoch::Ref(MEpoch::TAI,frame),
579 0 : MEpoch::Ref(MEpoch::LAST,frame));
580 0 : last = toLAST(epoch);
581 :
582 0 : return mac;
583 0 : }
584 : //
585 : //---------------------------------------------------------------
586 : //
587 0 : int AWProjectFT::findPointingOffsets(const VisBuffer2& vb,
588 : Array<Float> &l_off,
589 : Array<Float> &m_off,
590 : Bool Evaluate)
591 : {
592 : // LogIO log_l(LogOrigin("AWProjectFT2", "findPointingOffsets[R&D]"));
593 0 : Int NAnt = 0;
594 0 : MEpoch LAST;
595 0 : Double thisTime = getCurrentTimeStamp(vb);
596 : // Array<Float> pointingOffsets = epJ->nearest(thisTime);
597 0 : if (epJ_p.null()) return 0;
598 0 : Array<Float> pointingOffsets; epJ_p->nearest(thisTime,pointingOffsets);
599 0 : NAnt=pointingOffsets.shape()(2);
600 0 : l_off.resize(IPosition(3,1,1,NAnt)); // Poln x NChan x NAnt
601 0 : m_off.resize(IPosition(3,1,1,NAnt)); // Poln x NChan x NAnt
602 : // Can't figure out how to do the damn slicing of [Pol,NChan,NAnt,1] array
603 : // into [Pol,NChan,NAnt] array
604 : //
605 0 : IPosition tndx(3,0,0,0), sndx(4,0,0,0,0);
606 0 : for(tndx(2)=0;tndx(2)<NAnt; tndx(2)++,sndx(2)++)
607 : {
608 0 : sndx(0)=0; l_off(tndx) = pointingOffsets(sndx);
609 0 : sndx(0)=2; m_off(tndx) = pointingOffsets(sndx);
610 : }
611 0 : return NAnt;
612 : if (!Evaluate) return NAnt;
613 :
614 : //
615 : // Make a Coordinate Conversion Machine to go from (Az,El) to
616 : // (HA,Dec).
617 : //
618 : MDirection::Convert toAzEl = makeCoordinateMachine(vb,MDirection::HADEC,
619 : MDirection::AZEL,
620 : LAST);
621 : MDirection::Convert toHADec = makeCoordinateMachine(vb,MDirection::AZEL,
622 : MDirection::HADEC,
623 : LAST);
624 : //
625 : // ...and now hope that it all works and works correctly!!!
626 : //
627 : Quantity dAz(0,Radian),dEl(0,Radian);
628 : //
629 : // An array of shape [2,1,1]!
630 : //
631 : Array<Double> phaseDir = vb.subtableColumns().field().phaseDir().getColumn();
632 : Double RA0 = phaseDir(IPosition(3,0,0,0));
633 : Double Dec0 = phaseDir(IPosition(3,1,0,0));
634 : //
635 : // Compute reference (HA,Dec)
636 : //
637 : Double LST = LAST.get(Day).getValue();
638 : Double SDec0 = sin(Dec0), CDec0=cos(Dec0);
639 : LST -= floor(LST); // Extract the fractional day
640 : LST *= 2*C::pi;// Convert to Raidan
641 :
642 : Double HA0;
643 : HA0 = LST - RA0;
644 : Quantity QHA0(HA0,Radian), QDEC0(Dec0,Radian);
645 : //
646 : // Convert reference (HA,Dec) to reference (Az,El)
647 : //
648 : MDirection PhaseCenter(QHA0, QDEC0,MDirection::Ref(MDirection::HADEC));
649 : MDirection AzEl0 = toAzEl(PhaseCenter);
650 :
651 : MDirection tmpHADec = toHADec(AzEl0);
652 :
653 : Double Az0_Rad = AzEl0.getAngle(Radian).getValue()(0);
654 : Double El0_Rad = AzEl0.getAngle(Radian).getValue()(1);
655 :
656 : //
657 : // Convert the antenna pointing offsets from (Az,El)-->(RA,Dec)-->(l,m)
658 : //
659 :
660 : for(IPosition n(3,0,0,0);n(2)<=NAnt;n(2)++)
661 : {
662 : //
663 : // From (Az,El) -> (HA,Dec)
664 : //
665 : // Add (Az,El) offsets to the reference (Az,El)
666 : //
667 : dAz.setValue(l_off(n)+Az0_Rad); dEl.setValue(m_off(n)+El0_Rad);
668 : // dAz.setValue(0.0+Az0_Rad); dEl.setValue(0.0+El0_Rad);
669 : MDirection AzEl(dAz,dEl,MDirection::Ref(MDirection::AZEL));
670 : //
671 : // Convert offsetted (Az,El) to (HA,Dec) and then to (RA,Dec)
672 : //
673 : MDirection HADec = toHADec(AzEl);
674 : Double HA,Dec,RA, dRA;
675 : HA = HADec.getAngle(Radian).getValue()(0);
676 : Dec = HADec.getAngle(Radian).getValue()(1);
677 : RA = LST - HA;
678 : dRA = RA - RA0;
679 : //
680 : // Convert offsetted (RA,Dec) -> (l,m)
681 : //
682 : l_off(n) = sin(dRA)*cos(Dec);
683 : m_off(n) = sin(Dec)*CDec0-cos(Dec)*SDec0*cos(dRA);
684 : }
685 :
686 : return NAnt+1;
687 0 : }
688 : //
689 : //---------------------------------------------------------------
690 : //
691 0 : int AWProjectFT::findPointingOffsets(const VisBuffer2& vb,
692 : Cube<Float>& pointingOffsets,
693 : Array<Float> &l_off,
694 : Array<Float> &m_off,
695 : Bool Evaluate)
696 : {
697 : // LogIO log_l(LogOrigin("AWProjectFT2", "findPointingOffsets[R&D]"));
698 0 : Int NAnt = 0;
699 : //Float tmp;
700 : // TBD: adapt the following to VisCal mechanism:
701 0 : MEpoch LAST;
702 :
703 0 : NAnt=pointingOffsets.shape()(2);
704 0 : l_off.resize(IPosition(3,2,1,NAnt));
705 0 : m_off.resize(IPosition(3,2,1,NAnt));
706 0 : IPosition ndx(3,0,0,0),ndx1(3,0,0,0);
707 0 : for(ndx(2)=0;ndx(2)<NAnt;ndx(2)++)
708 : {
709 0 : ndx1=ndx;
710 0 : ndx(0)=0;ndx1(0)=0; //tmp=l_off(ndx) = pointingOffsets(ndx1);//Axis_0,Pol_0,Ant_i
711 0 : ndx(0)=1;ndx1(0)=1; //tmp=l_off(ndx) = pointingOffsets(ndx1);//Axis_0,Pol_1,Ant_i
712 0 : ndx(0)=0;ndx1(0)=2; //tmp=m_off(ndx) = pointingOffsets(ndx1);//Axis_1,Pol_0,Ant_i
713 0 : ndx(0)=1;ndx1(0)=3; //tmp=m_off(ndx) = pointingOffsets(ndx1);//Axis_1,Pol_1,Ant_i
714 : }
715 :
716 0 : return NAnt;
717 : if (!Evaluate) return NAnt;
718 :
719 : //
720 : // Make a Coordinate Conversion Machine to go from (Az,El) to
721 : // (HA,Dec).
722 : //
723 : MDirection::Convert toAzEl = makeCoordinateMachine(vb,MDirection::HADEC,
724 : MDirection::AZEL,
725 : LAST);
726 : MDirection::Convert toHADec = makeCoordinateMachine(vb,MDirection::AZEL,
727 : MDirection::HADEC,
728 : LAST);
729 : //
730 : // ...and now hope that it all works and works correctly!!!
731 : //
732 : Quantity dAz(0,Radian),dEl(0,Radian);
733 : //
734 : // An array of shape [2,1,1]!
735 : //
736 : Array<Double> phaseDir = vb.subtableColumns().field().phaseDir().getColumn();
737 : Double RA0 = phaseDir(IPosition(3,0,0,0));
738 : Double Dec0 = phaseDir(IPosition(3,1,0,0));
739 : //
740 : // Compute reference (HA,Dec)
741 : //
742 : Double LST = LAST.get(Day).getValue();
743 : Double SDec0 = sin(Dec0), CDec0=cos(Dec0);
744 : LST -= floor(LST); // Extract the fractional day
745 : LST *= 2*C::pi;// Convert to Raidan
746 :
747 : Double HA0;
748 : HA0 = LST - RA0;
749 : Quantity QHA0(HA0,Radian), QDEC0(Dec0,Radian);
750 : //
751 : // Convert reference (HA,Dec) to reference (Az,El)
752 : //
753 : MDirection PhaseCenter(QHA0, QDEC0,MDirection::Ref(MDirection::HADEC));
754 : MDirection AzEl0 = toAzEl(PhaseCenter);
755 :
756 : MDirection tmpHADec = toHADec(AzEl0);
757 :
758 : Double Az0_Rad = AzEl0.getAngle(Radian).getValue()(0);
759 : Double El0_Rad = AzEl0.getAngle(Radian).getValue()(1);
760 :
761 : //
762 : // Convert the antenna pointing offsets from (Az,El)-->(RA,Dec)-->(l,m)
763 : //
764 :
765 : for(IPosition n(3,0,0,0);n(2)<=NAnt;n(2)++)
766 : {
767 : //
768 : // From (Az,El) -> (HA,Dec)
769 : //
770 : // Add (Az,El) offsets to the reference (Az,El)
771 : //
772 : dAz.setValue(l_off(n)+Az0_Rad); dEl.setValue(m_off(n)+El0_Rad);
773 : // dAz.setValue(0.0+Az0_Rad); dEl.setValue(0.0+El0_Rad);
774 : MDirection AzEl(dAz,dEl,MDirection::Ref(MDirection::AZEL));
775 : //
776 : // Convert offsetted (Az,El) to (HA,Dec) and then to (RA,Dec)
777 : //
778 : MDirection HADec = toHADec(AzEl);
779 : Double HA,Dec,RA, dRA;
780 : HA = HADec.getAngle(Radian).getValue()(0);
781 : Dec = HADec.getAngle(Radian).getValue()(1);
782 : RA = LST - HA;
783 : dRA = RA - RA0;
784 : //
785 : // Convert offsetted (RA,Dec) -> (l,m)
786 : //
787 : l_off(n) = sin(dRA)*cos(Dec);
788 : m_off(n) = sin(Dec)*CDec0-cos(Dec)*SDec0*cos(dRA);
789 : }
790 :
791 : return NAnt+1;
792 0 : }
793 : //
794 : //---------------------------------------------------------------
795 : // The method below is for making the sensitivty image directly from
796 : // image-plane models of the PB. This is efficiently done via
797 : // ConvolutionFunction classes, at least for the ray traced models.
798 : //
799 : // For Wide-band AWP, this is computed in a wide-band sense by
800 : // accumulating WT*CF and FT'ing the resulting complex grid. This
801 : // is relatively more expensive but is a one-time cost. This method
802 : // (defined in C++ as the name plus the function signature) is
803 : // therefore overloaded in AWProjectWBFT which, in case the
804 : // sensitivity pattern (.pb image) is not found, only issues a
805 : // message informing that the first gridding cycle will do this
806 : // accumulation and therefore will be slower than the subsequent
807 : // ones. The steps necessary to compute the WB sensitivity pattern are:
808 : //
809 : // 1. Accumulate CFs (done via AWProjectWBFT::resampleCFToGrid())
810 : // 2. FT the grid (done via AWProjectWBFT::ftWeightImage()).
811 : // This gives PBs in the feed pol basis.
812 : // 3. Compute the Stokes PBs from feed basis
813 : //
814 : // The accumulation is done via the
815 : // AWProjectWBFT::resampleCFToGrid() method (which, for accumulating
816 : // the CFs, is the equivalent of ::resampledDataToGrid()). FT and
817 : // averaging across feed-pol planes is done in the following method:
818 : //
819 : // AWProjectWBFT::makeWBSensitivityImage(Lattice<T>&,
820 : // ImageInterface<Float>&,
821 : // const Matrix<Float>&,
822 : // const Bool&)
823 : //
824 : // (which uses the AWProjectWBFT::ftWeightImage() method for FFT).
825 : // This method is expected to be a NoOp if avgPBReady_p state
826 : // variable is True. Since these operations can be done only at the
827 : // _end_ of the (first) gridding cycle, they are naturally triggered
828 : // in the ImageInterface<Complex>& AWProjectWBFT::getImage() method,
829 : // which is also overloaded to first call
830 : // AWProjectWBFT::makeWBSensitivityImage() (which makes the PB and
831 : // points the class variables avgPB_p to it). And then call
832 : // AWProjectFT::getImage() (which normalizes the sky images).
833 : //
834 0 : void AWProjectFT::makeSensitivityImage(const VisBuffer2& vb,
835 : const ImageInterface<Complex>& imageTemplate,
836 : ImageInterface<Float>& sensitivityImage)
837 : {
838 0 : if (convFuncCtor_p->makeAverageResponse(vb, imageTemplate, sensitivityImage))
839 0 : cfCache_p->flush(sensitivityImage,sensitivityPatternQualifierStr_p);
840 0 : }
841 : //
842 : //---------------------------------------------------------------
843 : //
844 27035 : void AWProjectFT::makeCFPolMap(const VisBuffer2& vb, const Vector<Int>& locCfStokes,
845 : Vector<Int>& polM)
846 : {
847 27035 : Vector<Int> msStokes = vb.correlationTypes();
848 27035 : Int nPol = msStokes.nelements();
849 27035 : polM.resize(polMap.shape());
850 27035 : polM = -1;
851 :
852 81105 : for(Int i=0;i<nPol;i++)
853 81105 : for(uInt j=0;j<locCfStokes.nelements();j++)
854 81105 : if (locCfStokes(j) == msStokes(i))
855 54070 : {polM(i) = j;break;}
856 27035 : }
857 : //
858 : //---------------------------------------------------------------
859 : //
860 : // Given a polMap (mapping of which Visibility polarization is
861 : // gridded onto which grid plane), make a map of the conjugate
862 : // planes of the grid E.g, for Stokes-I and -V imaging, the two
863 : // planes of the uv-grid are [LL,RR]. For input VisBuffer2
864 : // visibilites in order [RR,RL,LR,LL], polMap = [1,-1,-1,0]. The
865 : // conjugate map will be [0,-1,-1,1].
866 : //
867 27035 : void AWProjectFT::makeConjPolMap(const VisBuffer2& vb,
868 : const Vector<Int> cfPolMap,
869 : Vector<Int>& conjPolMap)
870 : {
871 27035 : if (conjPolMap.nelements() > 0) return;
872 :
873 282 : LogIO log_l(LogOrigin("AWProjectFT2", "makConjPolMap[R&D]"));
874 :
875 : //
876 : // All the Natak (Drama) below with slicers etc. is to extract the
877 : // Poln. info. for the first IF only (not much "information
878 : // hiding" for the code to slice arrays in a general fashion).
879 : //
880 : // Extract the shape of the array to be sliced.
881 : //
882 : // Array<Int> stokesForAllIFs = vb.subtableColumns().polarization().corrType().getColumn();
883 : log_l << "############....temp code!!!!!!!!!! "
884 282 : << SynthesisUtils::mapSpwIDToPolID(vb, vb.spectralWindows()[0])
885 282 : << LogIO::DEBUG1;
886 141 : Vector<Int> polIDs=SynthesisUtils::mapSpwIDToPolID(vb,vb.spectralWindows()[0]);
887 141 : if (polIDs.nelements()==0)
888 0 : log_l << "Internal Error: Selected SPW did not map to any pol!" << LogIO::EXCEPTION;
889 : //
890 : // Use the first selected Spw to determine the pol. mapping in the
891 : // VB. This implies that the code assumes that all selected SPWs
892 : // have the same pol. setup.
893 : //
894 141 : Int polID=polIDs(0);
895 141 : Array<Int> stokesForAllIFs = (vb.subtableColumns()).polarization().corrType().get(polID);
896 141 : IPosition stokesShape(stokesForAllIFs.shape());
897 141 : IPosition firstIFStart(stokesShape),firstIFLength(stokesShape);
898 : //
899 : // Set up the start and length IPositions to extract only the
900 : // first column of the array. The following is required since the
901 : // array could have only one column as well.
902 : //
903 141 : firstIFStart(0)=0;firstIFLength(0)=stokesShape(0);
904 141 : for(uInt i=1;i<stokesShape.nelements();i++) {firstIFStart(i)=0;firstIFLength(i)=1;}
905 : //
906 : // Construct the slicer and produce the slice. .nonDegenerate
907 : // required to ensure the result of slice is a pure vector.
908 : //
909 282 : Vector<Int> visStokes = stokesForAllIFs(Slicer(firstIFStart,firstIFLength)).nonDegenerate();
910 :
911 141 : conjPolMap = cfPolMap;
912 :
913 141 : Int i,j,N = cfPolMap.nelements();
914 423 : for(i=0;i<N;i++)
915 282 : if (cfPolMap[i] > -1)
916 : {
917 282 : if (visStokes[i] == Stokes::RR)
918 : {
919 141 : conjPolMap[i]=-1;
920 282 : for(j=0;j<N;j++) if (visStokes[j] == Stokes::LL) break;
921 141 : conjPolMap[i]=cfPolMap[j];
922 : }
923 141 : else if (visStokes[i] == Stokes::LL)
924 : {
925 141 : conjPolMap[i]=-1;
926 141 : for(j=0;j<N;j++) if (visStokes[j] == Stokes::RR) break;
927 141 : conjPolMap[i]=cfPolMap[j];
928 : }
929 0 : else if (visStokes[i] == Stokes::LR)
930 : {
931 0 : conjPolMap[i]=-1;
932 0 : for(j=0;j<N;j++) if (visStokes[j] == Stokes::RL) break;
933 0 : conjPolMap[i]=cfPolMap[j];
934 : }
935 0 : else if (visStokes[i] == Stokes::RL)
936 : {
937 0 : conjPolMap[i]=-1;
938 0 : for(j=0;j<N;j++) if (visStokes[j] == Stokes::LR) break;
939 0 : conjPolMap[i]=cfPolMap[j];
940 : }
941 : }
942 141 : }
943 : //
944 : //---------------------------------------------------------------
945 : // Convert the CFs in the supplied CFStore to
946 : // A(nu)<Convolution>A(nu_*)
947 : //
948 0 : void AWProjectFT::makeWBCFWt(CFStore2& cfs, const Double imRefFreq)
949 : {
950 0 : LogIO log_l(LogOrigin("AWProjectFT2", "makeWBCFWt[R&D]"));
951 0 : log_l << "Converting WTCFs to wide-band versions" << LogIO::POST;
952 :
953 0 : Vector<Int> ant1List, ant2List;
954 0 : Vector<Quantity> paList;
955 0 : ant1List = cfs.getAnt1List();
956 0 : ant2List = cfs.getAnt2List();
957 0 : paList = cfs.getPAList();
958 :
959 0 : if (paNdxProcessed_p.nelements() == 0) {paNdxProcessed_p.resize(1); paNdxProcessed_p[0]=false;}
960 0 : CountedPtr<CFBuffer> cfb_l, cfb_clone;
961 0 : Quantity dPA(360.0,"deg");
962 0 : for (uInt pa=0;pa<paList.nelements();pa++)
963 0 : for (uint a1=0;a1<ant1List.nelements(); a1++)
964 0 : for (uint a2=0;a2<ant2List.nelements(); a2++)
965 : {
966 0 : if (paNdxProcessed_p.nelements() < pa) {paNdxProcessed_p.resize(pa+1,true); paNdxProcessed_p[pa]=false;}
967 0 : if (!paNdxProcessed_p[pa])
968 : {
969 0 : paNdxProcessed_p[pa]=true;
970 0 : Vector<Double> wVals, fVals; PolMapType mVals, mNdx, conjMVals, conjMNdx;
971 : Double fIncr, wIncr;
972 0 : cfb_l = cfs.getCFBuffer(paList[pa], dPA, ant1List(a1), ant2List(a2));
973 0 : cfb_l->getCoordList(fVals,wVals,mNdx, mVals, conjMNdx, conjMVals, fIncr, wIncr);
974 :
975 0 : cfb_clone=cfb_l->clone();
976 0 : cfb_clone->getCoordList(fVals,wVals,mNdx, mVals, conjMNdx, conjMVals, fIncr, wIncr);
977 :
978 0 : CountedPtr<Array<Complex> > cfc_l, conjCFC_l, result_l;
979 : //
980 : // Damn! Convolver does not work for complex convolution!
981 : // Convolver<Complex> convolver;
982 :
983 0 : for (Int nw=0;nw<(Int)wVals.nelements(); nw++)
984 0 : for (Int nf=0;nf<(Int)fVals.nelements(); nf++)
985 0 : for (Int ipol=0;ipol<(Int)conjMVals.nelements();ipol++)
986 0 : for (Int mRow=0;mRow<(Int)conjMVals[ipol].nelements(); mRow++)
987 : {
988 : Double f, cf;
989 0 : f=fVals[nf];
990 0 : cf=sqrt(2*imRefFreq*imRefFreq - f*f);
991 0 : Int conjNF = cfb_l->nearestFreqNdx(cf);
992 :
993 0 : result_l = cfb_l->getCFCellPtr(nf,nw,conjMNdx[ipol][mRow])->getStorage();
994 0 : cfc_l = cfb_clone->getCFCellPtr(nf,nw,conjMNdx[ipol][mRow])->getStorage();
995 0 : conjCFC_l = cfb_clone->getCFCellPtr(conjNF,nw,conjMNdx[ipol][mRow])->getStorage();
996 :
997 0 : SynthesisUtils::libreConvolver(*result_l,*conjCFC_l);
998 : }
999 0 : }
1000 : }
1001 0 : }
1002 : //
1003 : // Locate a convlution function. It will be either in the cache
1004 : // (mem. or disk cache) or will be computed and cached for possible
1005 : // later use.
1006 : //
1007 27778 : void AWProjectFT::findConvFunction(const ImageInterface<Complex>& image,
1008 : const VisBuffer2& vb)
1009 : {
1010 27778 : if (!paChangeDetector.changed(vb,0)) return;
1011 211 : Int cfSource=CFDefs::NOTCACHED;
1012 211 : CoordinateSystem ftcoords;
1013 : // Think of a generic call to get the key-values. And a
1014 : // overloadable method (or an externally supplied one?) to convert
1015 : // the values to key-ids. That will ensure that AWProjectFT
1016 : // remains the A-Projection algorithm implementation configurable
1017 : // by the behaviour of the supplied objects.
1018 211 : Float pa=getVBPA(vb);
1019 : //UUU// ok();
1020 211 : visResampler_p->setMaps(chanMap, polMap); //UUU Added here.
1021 211 : visResampler_p->setFreqMaps(expandedSpwFreqSel_p,expandedSpwConjFreqSel_p);
1022 :
1023 211 : lastPAUsedForWtImg = currentCFPA = pa;
1024 :
1025 : //Vector<Vector<Double> > pointingOffset(convFuncCtor_p->findPointingOffset(image,vb, doPointing));
1026 :
1027 : // PO::setOverSampling needs to be called here since
1028 : // convFuncCtor_p gets that value from a combination of (1)
1029 : // ATerm_OVERSAMPLING env. variable, (2) ATERM.OVERSAMPLING in
1030 : // ~/.casa and (3) from existing CFCache. This setting in the AWP
1031 : // constructor will only get the default value from ATerm.h
1032 211 : convSampling=convFuncCtor_p->getOversampling();
1033 211 : po_p->setOverSampling(convFuncCtor_p->getOversampling());
1034 : // PO::fetchPointingOffset() only updates the internal cache in PO
1035 : // class. PO::pullPointingOffset() is required to extract in the
1036 : // calling class.
1037 211 : po_p->fetchPointingOffset(image,vb, doPointing);
1038 :
1039 211 : Float dPA = paChangeDetector.getParAngleTolerance().getValue("rad");
1040 211 : Quantity dPAQuant = Quantity(paChangeDetector.getParAngleTolerance());
1041 :
1042 211 : vb2CFBMap_p->setDoPointing(doPointing);
1043 422 : cfSource = vb2CFBMap_p->makeVBRow2CFBMap(*cfs2_p,
1044 : vb,
1045 : dPAQuant,
1046 211 : chanMap,polMap,po_p);
1047 :
1048 211 : if (cfSource == CFDefs::NOTCACHED)
1049 : {
1050 14 : PolMapType polMat, polIndexMat, conjPolMat, conjPolIndexMat;
1051 14 : Vector<Int> visPolMap(vb.correlationTypes());
1052 14 : polMat = pop_p->makePolMat(visPolMap,polMap);
1053 14 : polIndexMat = pop_p->makePol2CFMat(visPolMap,polMap);
1054 :
1055 14 : conjPolMat = pop_p->makeConjPolMat(visPolMap,polMap);
1056 14 : conjPolIndexMat = pop_p->makeConjPol2CFMat(visPolMap,polMap);
1057 :
1058 14 : convFuncCtor_p->setPolMap(polMap);
1059 14 : convFuncCtor_p->setSpwSelection(spwChanSelFlag_p);
1060 14 : convFuncCtor_p->setSpwFreqSelection(spwFreqSel_p);
1061 :
1062 : // USEFUL DEBUG MESSAGE
1063 : //cerr << "Freq. selection: " << expandedSpwFreqSel_p << endl << expandedSpwConjFreqSel_p << endl;
1064 14 : Bool pleaseDoAlsoFillTheCF=!dryRun();
1065 14 : convFuncCtor_p->makeConvFunction(image,vb,wConvSize,
1066 14 : pop_p, pa, dPA, uvScale, uvOffset,spwFreqSel_p,
1067 : *cfs2_p, *cfwts2_p, pleaseDoAlsoFillTheCF);
1068 14 : }
1069 : //
1070 : // If one-time-operations in the CFCache not yet done, set the
1071 : // pol. index maps in the CFCache.
1072 : //
1073 211 : if (!cfCache_p->OTODone())
1074 : {
1075 : //Vector<Int> visPolMap(vb.corrType());
1076 79 : Vector<Int> visPolMap(vb.correlationTypes());
1077 :
1078 79 : PolMapType polMat, conjPolMat;
1079 79 : polMat = pop_p->makePolMat(visPolMap,polMap);
1080 79 : conjPolMat = pop_p->makeConjPolMat(visPolMap,polMap);
1081 :
1082 79 : PolMapType pNdx, cpNdx;
1083 79 : pNdx = pop_p->makePol2CFMat(visPolMap,polMap);
1084 79 : cpNdx = pop_p->makeConjPol2CFMat(visPolMap,polMap);
1085 :
1086 79 : cfCache_p->initPolMaps(pNdx,cpNdx);
1087 :
1088 79 : cfs2_p->initMaps(vb,spwFreqSel_p,imRefFreq_p);
1089 79 : cfwts2_p->initMaps(vb,spwFreqSel_p,imRefFreq_p);
1090 79 : }
1091 : //
1092 : // Load the average PB (sensitivity pattern) from the cache. If
1093 : // not found in the cache, make one and cache it.
1094 : //
1095 211 : std::tuple<int, double>cubeinfo(1,-1.0);
1096 : double freqofBegChan;
1097 211 : spectralCoord_p.toWorld(freqofBegChan, 0.0);
1098 :
1099 211 : cubeinfo=std::make_tuple(image.shape()(3),freqofBegChan);
1100 :
1101 211 : if(!avgPBReady_p)
1102 77 : avgPBReady_p = (cfCache_p->loadAvgPB(avgPB_p,sensitivityPatternQualifierStr_p, cubeinfo) != CFDefs::NOTCACHED);
1103 :
1104 211 : if(avgPBReady_p){
1105 165 : LatticeExprNode le( max( *avgPB_p ) );
1106 165 : Float avgPB_max=le.getFloat();
1107 :
1108 165 : if(avgPB_max <= 0.0) avgPBReady_p = false;
1109 165 : }
1110 :
1111 211 : if(!avgPBReady_p) makeSensitivityImage(vb,image,*avgPB_p);
1112 :
1113 :
1114 :
1115 211 : verifyShapes(avgPB_p->shape(), image.shape());
1116 :
1117 211 : if (paChangeDetector.changed(vb,0)) paChangeDetector.update(vb,0);
1118 : //
1119 : // Write some useful info. to the logger.
1120 : //
1121 211 : if (cfSource != CFDefs::MEMCACHE)
1122 : {
1123 : // If dry run, write the uvgrid as an image for later use in
1124 : // filling the empty CFCache. Only the co-ordinate system of
1125 : // the uvgrid is required later.
1126 14 : if (dryRun())
1127 : {
1128 7 : PagedImage<Complex> thisGrid(image.shape(),image.coordinates(),
1129 21 : cfCache_p->getCacheDir()+"/uvgrid.im");
1130 7 : }
1131 : // Save only the CF Cube for the current value of PA (not the
1132 : // entire CFStore -- CFs for PA values encountered earlier
1133 : // than current value have already need made persistent).
1134 14 : cfs2_p->makePersistent(cfCache_p->getCacheDir().c_str(),"","", Quantity(pa,"rad"),dPAQuant,0,0);
1135 14 : cfwts2_p->makePersistent(cfCache_p->getCacheDir().c_str(),"","WT",Quantity(pa,"rad"),dPAQuant,0,0);
1136 14 : Double memUsed=cfs2_p->memUsage();
1137 14 : String unit(" KB");
1138 14 : memUsed = (Int)(memUsed/1024.0+0.5);
1139 14 : if (memUsed > 1024) {memUsed /=1024; unit=" MB";}
1140 :
1141 28 : LogIO log_l(LogOrigin("AWProjectFT2", "findConvFunction[R&D]"));
1142 : log_l << "Convolution function memory footprint:"
1143 : << (Int)(memUsed) << unit << " out of a maximum of "
1144 14 : << HostInfo::memoryTotal(true)/1024 << " MB" << LogIO::POST;
1145 :
1146 : //
1147 : // Initialize any internal maps that may be used later for
1148 : // efficient access.
1149 : //
1150 14 : cfs2_p->initMaps(vb,spwFreqSel_p,imRefFreq_p);
1151 14 : cfwts2_p->initMaps(vb,spwFreqSel_p,imRefFreq_p);
1152 14 : }
1153 211 : }
1154 : //
1155 : //------------------------------------------------------------------------------
1156 : // Vectorized initializeToVis. See design related comments in the
1157 : // .h file.
1158 : //
1159 : //
1160 : // The functional goals here are:
1161 : //
1162 : // 1. If doPBCorrection==true (i.e. the input sky-image is flat-sky
1163 : // image), divide the image by the PB
1164 : //
1165 : // 2. Convert from Stokes to Feed (correlation) frame
1166 : //
1167 : // 3. FFT the image to produce gridded vis.
1168 : //
1169 : // 4. And since the same image buffer is used to accumulate the
1170 : // model, multiply the sky-image back with PB
1171 : //
1172 : // The call to non-vectorized version of initializeToVis() which is
1173 : // pure virtual and hence the local version is called, is only to do
1174 : // the FFT. FFT is *NOT* in place.
1175 : //
1176 : // These operations effecitvely maintains the model image as PB*Sky
1177 : // and supplies flat-sky model for prediction. This can probably be
1178 : // achieve with fewer operations and same memory buffers....but that
1179 : // for later (SB)
1180 0 : void AWProjectFT::initializeToVis(Block<CountedPtr<ImageInterface<Complex> > > & compImageVec,
1181 : PtrBlock<SubImage<Float> *> & modelImageVec,
1182 : PtrBlock<SubImage<Float> *>& weightImageVec,
1183 : PtrBlock<SubImage<Float> *>& /*fluxScaleVec*/,
1184 : Block<Matrix<Float> >& weightsVec,
1185 : const VisBuffer2& vb)
1186 : {
1187 0 : LogIO log_p(LogOrigin("AWProjectFT2","initToVis[V][R&D]"));
1188 : //
1189 : // Setting the image below is crucial since init() and
1190 : // initMaps(vb) below expect this to be set.
1191 : //
1192 0 : image=&(*compImageVec[0]);
1193 :
1194 : log_p << "Total flux in model image (before avgPB normalization): "
1195 0 : << sum((*(modelImageVec[0])).get())
1196 0 : << LogIO::POST;
1197 0 : if(doPBCorrection)
1198 : {
1199 : // Make the sensitivity Image if applicable
1200 0 : init();
1201 0 : initMaps(vb);
1202 0 : findConvFunction(*(compImageVec[0]), vb); // Pure virtual -- call local version
1203 :
1204 0 : if (isDryRun) return;
1205 :
1206 : // Get the sensitivity Image
1207 0 : Matrix<Float> tempWts;
1208 0 : tempWts.resize();
1209 0 : getWeightImage(*(weightImageVec[0]), tempWts); // Pure virtual -- call local version
1210 :
1211 : // Normalize the model image by the sensitivity image only.
1212 : // No local implementation -- call FTMachine version
1213 :
1214 : // Divide by avgPB ///// PBWeight
1215 : //
1216 :
1217 : // Divide by sqrt(avgPB) ////// PBSQWeight
1218 : //
1219 0 : normalizeImage( *(modelImageVec[0]) , weightsVec[0], *(weightImageVec[0]),
1220 0 : false, (Float)pbLimit_p, (Int)4);
1221 0 : }
1222 : log_p << "Total flux in model image (after avgPB normalization): "
1223 0 : << sum((*(modelImageVec[0])).get()) << LogIO::POST;
1224 :
1225 : // Convert from Stokes planes to Correlation planes
1226 : // No local implementation -- call FTMachine version
1227 0 : stokesToCorrelation(*(modelImageVec[0]), *(compImageVec[0]));
1228 :
1229 : // Call initializeToVis
1230 0 : initializeToVis(*(compImageVec[0]), vb); // Pure virtual
1231 :
1232 : // Multiply the flat-sky model by the PB.
1233 : // No local implementation -- call FTMachine version
1234 0 : if(doPBCorrection)
1235 : // Multiply by avgPB /////// PBWeight
1236 : //normalizeImage( *(modelImageVec[0]) , weightsVec[0], *(weightImageVec[0]) , false, (Float)pbLimit_p, (Int)3);
1237 :
1238 : // Multiply by sqrt(avgPB) ///// PBSQWeight
1239 0 : normalizeImage( *(modelImageVec[0]) , weightsVec[0], *(weightImageVec[0]) , false, (Float)pbLimit_p, (Int)5);
1240 0 : }
1241 : //
1242 : //------------------------------------------------------------------------------
1243 : //
1244 35 : void AWProjectFT::initializeToVis(ImageInterface<Complex>& iimage,
1245 : const VisBuffer2& vb)
1246 : {
1247 70 : LogIO log_l(LogOrigin("AWProjectFT2", "initializeToVis[R&D]"));
1248 35 : image=&iimage;
1249 :
1250 35 : ok();
1251 :
1252 35 : init();
1253 35 : makingPSF = false;
1254 35 : initMaps(vb);
1255 :
1256 35 : findConvFunction(*image, vb);
1257 35 : if (isDryRun) return;
1258 : //
1259 : // Initialize the maps for polarization and channel. These maps
1260 : // translate visibility indices into image indices
1261 : //
1262 :
1263 35 : nx = image->shape()(0);
1264 35 : ny = image->shape()(1);
1265 35 : npol = image->shape()(2);
1266 35 : nchan = image->shape()(3);
1267 :
1268 : //
1269 : // If we are memory-based then read the image in and create an
1270 : // ArrayLattice otherwise just use the PagedImage
1271 : //
1272 :
1273 35 : isTiled=false;
1274 :
1275 : {
1276 35 : IPosition gridShape(4, nx, ny, npol, nchan);
1277 35 : griddedData.resize(gridShape);
1278 35 : griddedData=Complex(0.0);
1279 :
1280 35 : IPosition stride(4, 1);
1281 70 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
1282 70 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
1283 70 : IPosition trc(blc+image->shape()-stride);
1284 :
1285 35 : IPosition start(4, 0);
1286 35 : griddedData(blc, trc) = image->getSlice(start, image->shape());
1287 :
1288 35 : lattice=new ArrayLattice<Complex>(griddedData);
1289 35 : }
1290 :
1291 : //AlwaysAssert(lattice, AipsError);
1292 :
1293 35 : log_l << LogIO::DEBUGGING << "Starting FFT of image" << LogIO::POST;
1294 :
1295 35 : Vector<Float> sincConv(nx);
1296 35 : Float centerX=nx/2;
1297 33315 : for (Int ix=0;ix<nx;ix++)
1298 : {
1299 33280 : Float x=C::pi*Float(ix-centerX)/(Float(nx)*Float(convSampling));
1300 33280 : if(ix==centerX) sincConv(ix)=1.0;
1301 33245 : else sincConv(ix)=sin(x)/x;
1302 33280 : sincConv(ix) = 1.0;
1303 : }
1304 :
1305 35 : if(cfCache_p->avgPBReady()) //SB
1306 : {
1307 : // normalizeAvgPB();
1308 :
1309 20 : IPosition cursorShape(4, nx, 1, 1, 1);
1310 20 : IPosition axisPath(4, 0, 1, 2, 3);
1311 20 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
1312 20 : LatticeIterator<Complex> lix(*lattice, lsx);
1313 :
1314 20 : verifyShapes(avgPB_p->shape(), image->shape());
1315 20 : Array<Float> avgBuf; avgPB_p->get(avgBuf);
1316 : // If the total-power sensitivity pattern peak is too low, warn
1317 : // the user. This usually is indicative of a rat somewhere in
1318 : // the pipes upstream...
1319 20 : if ((sensitivityPatternQualifier_p==0) && (max(avgBuf) < 1e-04))
1320 : log_l << "Normalization by PB requested but either PB was not"
1321 : <<" found in the cache or is ill-formed. Peak = "
1322 0 : << max(avgBuf)// << " " << sensitivityPatternQualifier_p
1323 0 : << LogIO::WARN << LogIO::POST;
1324 :
1325 :
1326 20 : LatticeStepper lpb(avgPB_p->shape(),cursorShape,axisPath);
1327 20 : LatticeIterator<Float> lipb(*avgPB_p, lpb);
1328 :
1329 20 : Vector<Complex> griddedVis;
1330 : //
1331 : // Grid correct in anticipation of the convolution by the
1332 : // convFunc. Each polarization plane is corrected by the
1333 : // appropraite primary beam.
1334 : //
1335 22548 : for(lix.reset(),lipb.reset();!lix.atEnd();lix++,lipb++)
1336 : {
1337 22528 : Int iy=lix.position()(1);
1338 22528 : griddedVis = lix.rwVectorCursor();
1339 :
1340 22528 : Vector<Float> PBCorrection(lipb.rwVectorCursor().shape());
1341 22528 : PBCorrection = lipb.rwVectorCursor();
1342 30431232 : for(int ix=0;ix<nx;ix++)
1343 : {
1344 30408704 : if (doPBCorrection)
1345 : {
1346 30408704 : PBCorrection(ix) = pbFunc(PBCorrection(ix),pbLimit_p)*(sincConv(ix)*sincConv(iy));
1347 30408704 : lix.rwVectorCursor()(ix) /= (PBCorrection(ix));
1348 : }
1349 : else
1350 0 : lix.rwVectorCursor()(ix) /= (1.0/(sincConv(ix)*sincConv(iy)));
1351 : }
1352 22528 : }
1353 20 : }
1354 : //
1355 : // Now do the FFT2D in place
1356 : //
1357 :
1358 35 : LatticeFFT::cfft2d(*lattice);
1359 35 : log_l << LogIO::DEBUGGING << "Finished FFT" << LogIO::POST;
1360 35 : }
1361 : //
1362 : //------------------------------------------------------------------------------
1363 : //
1364 0 : void AWProjectFT::initializeToVis(ImageInterface<Complex>& iimage,
1365 : const VisBuffer2& vb,
1366 : Array<Complex>& griddedVis,
1367 : Vector<Double>& uvscale)
1368 : {
1369 0 : initializeToVis(iimage, vb);
1370 0 : griddedVis.assign(griddedData); //using the copy for storage
1371 0 : uvscale.assign(uvScale);
1372 :
1373 0 : }
1374 : //
1375 : //---------------------------------------------------------------
1376 : //
1377 35 : void AWProjectFT::finalizeToVis()
1378 : {
1379 35 : visResampler_p->runTimeDG_p=0.0;
1380 35 : logIO()<< LogIO::NORMAL2 << "Time degrid " << timedegrid_p << LogIO::POST;
1381 35 : timedegrid_p=0.0;
1382 :
1383 35 : if(!lattice.null()) lattice=0;
1384 35 : griddedData.resize();
1385 :
1386 35 : if(isTiled)
1387 : {
1388 0 : AlwaysAssert(imageCache, AipsError);
1389 0 : AlwaysAssert(image, AipsError);
1390 0 : ostringstream o;
1391 0 : imageCache->flush();
1392 0 : imageCache->showCacheStatistics(o);
1393 :
1394 0 : LogIO log_l(LogOrigin("AWProjectFT2", "finalizeToVis[R&D]"));
1395 0 : log_l << o.str() << LogIO::POST;
1396 0 : }
1397 35 : if(pointingToImage) delete pointingToImage;
1398 35 : pointingToImage=0;
1399 35 : }
1400 : //
1401 : //---------------------------------------------------------------
1402 : //
1403 : // Initialize the FFT to the Sky. Here we have to setup and
1404 : // initialize the grid.
1405 : //
1406 176 : void AWProjectFT::initializeToSky(ImageInterface<Complex>& iimage,
1407 : Matrix<Float>& weight,
1408 : const VisBuffer2& vb)
1409 : {
1410 352 : LogIO log_l(LogOrigin("AWProjectFT2", "initializeToSky[R&D]"));
1411 :
1412 : // image always points to the image
1413 176 : image=&iimage;
1414 :
1415 176 : init();
1416 176 : initMaps(vb);
1417 176 : log_l << "Computed maps using FTMachine::initMaps. " << "polMap = " << polMap << LogIO::POST;
1418 176 : visResampler_p->setMaps(chanMap, polMap);
1419 176 : visResampler_p->setFreqMaps(expandedSpwFreqSel_p,expandedSpwConjFreqSel_p);
1420 :
1421 : // Initialize the maps for polarization and channel. These maps
1422 : // translate visibility indices into image indices
1423 176 : nx = image->shape()(0);
1424 176 : ny = image->shape()(1);
1425 176 : npol = image->shape()(2);
1426 176 : nchan = image->shape()(3);
1427 :
1428 176 : sumWeight=0.0;
1429 176 : sumCFWeight = 0.0;
1430 176 : weight.resize(sumWeight.shape());
1431 176 : weight=0.0;
1432 :
1433 : //
1434 : // Construct the HPG. This is now donw in AWVRHPG. Not the right
1435 : // place, I think, but good enough for testing. (21Dec2020).
1436 : //
1437 : // const std::array<unsigned, 4> gridSize{(uInt)nx,(uInt)ny,(uInt)npol,(uInt)nchan};
1438 : // const std::array<float, 2> gridScale{(float)uvScale(0), (float)uvScale(1)};
1439 : // unsigned max_added_tasks=1; // This allocates 2 CUDA streams in AWVRHPG.
1440 : // visResampler_p->initGridder(max_added_tasks,gridSize,gridScale);
1441 :
1442 :
1443 : //
1444 : // Initialize for in memory or to disk gridding. lattice will
1445 : // point to the appropriate Lattice, either the ArrayLattice for
1446 : // in memory gridding or to the image for to disk gridding.
1447 : //
1448 176 : if(isTiled)
1449 : {
1450 0 : imageCache->flush();
1451 0 : image->set(Complex(0.0));
1452 0 : lattice=CountedPtr<Lattice<Complex> > (image, false);
1453 : }
1454 : else
1455 : {
1456 176 : IPosition gridShape(4, nx, ny, npol, nchan);
1457 176 : if(!useDoubleGrid_p){
1458 0 : griddedData.resize(gridShape);
1459 0 : griddedData=Complex(0.0);
1460 : }
1461 : else {
1462 176 : griddedData2.resize(gridShape);
1463 176 : griddedData2=DComplex(0.0);
1464 : }
1465 176 : }
1466 :
1467 : //cerr << "initializeToSky for grid" << endl;
1468 176 : if(useDoubleGrid_p)
1469 176 : visResampler_p->initializeToSky(griddedData2, sumWeight);
1470 : else
1471 0 : visResampler_p->initializeToSky(griddedData, sumWeight);
1472 176 : }
1473 : //
1474 : //---------------------------------------------------------------
1475 : //
1476 169 : void AWProjectFT::finalizeToSky()
1477 : {
1478 :
1479 169 : logIO() << LogIO::NORMAL2 << "time to massage data " << timemass_p << LogIO::POST;
1480 169 : logIO() << LogIO::NORMAL2 << "time gridding " << timegrid_p << LogIO::POST;
1481 169 : timemass_p=0.0;
1482 169 : timegrid_p=0.0;
1483 :
1484 : //
1485 : // Now we flush the cache and report statistics For memory based,
1486 : // we don't write anything out yet.
1487 : //
1488 : // LogIO log_l(LogOrigin("AWProjectFT2", "finalizeToSky[R&D]"));
1489 :
1490 169 : if(isTiled)
1491 : {
1492 0 : AlwaysAssert(image, AipsError);
1493 0 : AlwaysAssert(imageCache, AipsError);
1494 0 : imageCache->flush();
1495 0 : ostringstream o;
1496 0 : imageCache->showCacheStatistics(o);
1497 0 : }
1498 169 : if(pointingToImage) delete pointingToImage;
1499 169 : pointingToImage=0;
1500 :
1501 169 : paChangeDetector.reset();
1502 169 : cfCache_p->flush();
1503 169 : if(useDoubleGrid_p)
1504 169 : visResampler_p->finalizeToSky(griddedData2, sumWeight);
1505 : else
1506 0 : visResampler_p->finalizeToSky(griddedData, sumWeight);
1507 169 : }
1508 : //
1509 : //---------------------------------------------------------------
1510 : //
1511 0 : Array<Complex>* AWProjectFT::getDataPointer(const IPosition& centerLoc2D,
1512 : Bool readonly)
1513 : {
1514 : Array<Complex>* result;
1515 : // Is tiled: get tiles and set up offsets
1516 0 : centerLoc(0)=centerLoc2D(0);
1517 0 : centerLoc(1)=centerLoc2D(1);
1518 0 : result=&imageCache->tile(offsetLoc, centerLoc, readonly);
1519 0 : gridder->setOffset(IPosition(2, offsetLoc(0), offsetLoc(1)));
1520 0 : return result;
1521 : }
1522 :
1523 : // The following file has the runFORTRAN* stuff. Moving it to a
1524 : // separate file to reduce clutter and ultimately delete it.
1525 : #include "AWProjectFT.FORTRANSTUFF"
1526 : //
1527 : //---------------------------------------------------------------
1528 : //
1529 21641 : void AWProjectFT::put(const VisBuffer2& vb, Int /*row*/, Bool dopsf,
1530 : FTMachine::Type type)
1531 : {
1532 : // Take care of translation of Bools to Integer
1533 21641 : makingPSF=dopsf;
1534 21641 : if(dopsf)
1535 9154 : ftmType_p=refim::FTMachine::PSF;
1536 21641 : Timer tim;
1537 21641 : tim.mark();
1538 :
1539 : try
1540 : {
1541 21641 : findConvFunction(*image, vb);
1542 : }
1543 0 : catch(AipsError& x)
1544 : {
1545 0 : LogIO log_l(LogOrigin("AWProjectFT2", "put[R&D]"));
1546 0 : log_l << x.getMesg() << LogIO::WARN;
1547 0 : return;
1548 0 : }
1549 21641 : if (isDryRun) return;
1550 :
1551 20933 : Nant_p = vb.subtableColumns().antenna().nrow();
1552 :
1553 20933 : Matrix<Float> imagingweight;
1554 20933 : getImagingWeight(imagingweight, vb);
1555 :
1556 20933 : Cube<Complex> data;
1557 : //Fortran gridder need the flag as ints
1558 20933 : Cube<Int> flags;
1559 20933 : Matrix<Float> elWeight;
1560 :
1561 20933 : interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, type);
1562 :
1563 : //
1564 : // Get the uvws in a form that Fortran can use and do that
1565 : // necessary phase rotation. On a Pentium Pro 200 MHz when null,
1566 : // this step takes about 50us per uvw point. This is just barely
1567 : // noticeable for Stokes I continuum and irrelevant for other
1568 : // cases.
1569 : //
1570 20933 : Matrix<Double> uvw(negateUV(vb));
1571 :
1572 20933 : Vector<Double> dphase(vb.nRows());
1573 20933 : dphase=0.0;
1574 20933 : doUVWRotation_p=true;
1575 20933 : girarUVW(uvw, dphase, vb);
1576 20933 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
1577 :
1578 : //Here we redo the match or use previous match
1579 :
1580 : //Channel matching for the actual spectral window of buffer
1581 :
1582 20933 : matchChannel(vb);
1583 20933 : VBStore vbs;
1584 20933 : Vector<Int> gridShape = griddedData2.shape().asVector();
1585 20933 : setupVBStore(vbs,vb, elWeight,data,uvw,flags, dphase,dopsf,gridShape);
1586 20933 : timemass_p +=tim.real();
1587 20933 : tim.mark();
1588 :
1589 20933 : if (useDoubleGrid_p)
1590 : {
1591 20933 : resampleDataToGrid(griddedData2, vbs, vb, dopsf);//, *imagingweight, *data, uvw,flags,dphase,dopsf);
1592 : }
1593 : else
1594 : {
1595 0 : resampleDataToGrid(griddedData, vbs, vb, dopsf);//, *imagingweight, *data, uvw,flags,dphase,dopsf);
1596 : }
1597 20933 : timegrid_p+=tim.real();
1598 20933 : }
1599 :
1600 0 : std::shared_ptr<std::complex<double>> AWProjectFT::getGridPtr(size_t& size) const
1601 : {
1602 0 : return visResampler_p->getGridPtr(size);
1603 : }
1604 :
1605 0 : std::shared_ptr<double> AWProjectFT::getSumWeightsPtr(size_t& size) const
1606 : {
1607 0 : return visResampler_p->getSumWeightsPtr(size);
1608 : }
1609 :
1610 : //
1611 : //-------------------------------------------------------------------------
1612 : // Gridding
1613 0 : void AWProjectFT::resampleDataToGrid(Array<Complex>& griddedData_l, VBStore& vbs,
1614 : const VisBuffer2& /*vb*/, Bool& dopsf)
1615 : {
1616 0 : visResampler_p->DataToGrid(griddedData_l, vbs, sumWeight, dopsf);
1617 0 : }
1618 : //
1619 : //-------------------------------------------------------------------------
1620 : // Gridding
1621 20933 : void AWProjectFT::resampleDataToGrid(Array<DComplex>& griddedData_l, VBStore& vbs,
1622 : const VisBuffer2& /*vb*/, Bool& dopsf)
1623 : {
1624 20933 : visResampler_p->DataToGrid(griddedData_l, vbs, sumWeight, dopsf);
1625 20933 : }
1626 : //
1627 : //---------------------------------------------------------------
1628 : //
1629 6102 : void AWProjectFT::get(VisBuffer2& vb, Int /*row*/)
1630 : {
1631 6102 : findConvFunction(*image, vb);
1632 6102 : Timer tim;
1633 6102 : tim.mark();
1634 6102 : Nant_p = vb.subtableColumns().antenna().nrow();
1635 : // Get the uvws in a form that Fortran can use
1636 6102 : Matrix<Double> uvw(negateUV(vb));
1637 :
1638 6102 : Vector<Double> dphase(vb.nRows());
1639 6102 : dphase=0.0;
1640 6102 : doUVWRotation_p=true;
1641 : //rotateUVW(uvw, dphase, vb);
1642 :
1643 6102 : girarUVW(uvw, dphase, vb);
1644 6102 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
1645 :
1646 6102 : matchChannel(vb);
1647 : //No point in reading data if its not matching in frequency
1648 6102 : if(max(chanMap)==-1) return;
1649 :
1650 6102 : Cube<Complex> data;
1651 6102 : Cube<Int> flags;
1652 6102 : getInterpolateArrays(vb, data, flags);
1653 :
1654 6102 : VBStore vbs;
1655 6102 : Bool tmpDoPSF=false;
1656 :
1657 6102 : setupVBStore(vbs,vb, vb.imagingWeight(),data,uvw,flags, dphase,tmpDoPSF,griddedData.shape().asVector());
1658 6102 : resampleGridToData(vbs, griddedData, vb);//, uvw, flags, dphase);
1659 6102 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1660 6102 : timedegrid_p+=tim.real();
1661 6102 : }
1662 : //
1663 : //-------------------------------------------------------------------------
1664 : // De-gridding
1665 6102 : void AWProjectFT::resampleGridToData(VBStore& vbs, Array<Complex>& griddedData_l,
1666 : const VisBuffer2& /*vb*/)
1667 : {
1668 6102 : visResampler_p->GridToData(vbs, griddedData_l);
1669 6102 : }
1670 : //
1671 : //---------------------------------------------------------------
1672 : //
1673 : // Finalize the FFT to the Sky. Here we actually do the FFT and
1674 : // return the resulting image
1675 169 : ImageInterface<Complex>& AWProjectFT::getImage(Matrix<Float>& weights,
1676 : Bool fftNormalization)
1677 : {
1678 338 : LogIO log_l(LogOrigin("AWProjectFT2", "getImage[R&D]"));
1679 :
1680 169 : AlwaysAssert(image, AipsError);
1681 :
1682 169 : weights.resize(sumWeight.shape());
1683 169 : convertArray(weights, sumWeight);
1684 : //
1685 : // If the weights are all zero then we cannot normalize otherwise
1686 : // we don't care.
1687 : //
1688 169 : if(max(weights)==0.0)
1689 : log_l // UUU //<< LogIO::SEVERE
1690 0 : << "No useful data in " << name() << ". Weights all zero"
1691 0 : << LogIO::POST;
1692 : else
1693 : {
1694 169 : log_l << "Sum of weights: " << weights << " " << max(griddedData2) << " " << min(griddedData2) << LogIO::POST;
1695 : //cerr << "Sum of weights: " << setprecision(20) << weights << endl;
1696 : }
1697 : // UUU else
1698 : {
1699 : log_l << LogIO::DEBUGGING
1700 169 : << "Starting FFT and scaling of image" << LogIO::POST;
1701 : //
1702 : // x and y transforms (lattice has the gridded vis. Make the
1703 : // dirty images)
1704 : //
1705 169 : if (useDoubleGrid_p)
1706 : {
1707 169 : ArrayLattice<DComplex> darrayLattice(griddedData2);
1708 169 : LatticeFFT::cfft2d(darrayLattice,false);
1709 169 : griddedData.resize(griddedData2.shape());
1710 169 : convertArray(griddedData, griddedData2);
1711 169 : SynthesisUtilMethods::getResource("mem peak in getImage");
1712 :
1713 : //Don't need the double-prec grid anymore...
1714 169 : griddedData2.resize();
1715 169 : lattice=new ArrayLattice<Complex>(griddedData);
1716 169 : }
1717 : else
1718 : {
1719 0 : lattice=new ArrayLattice<Complex>(griddedData);
1720 0 : LatticeFFT::cfft2d(*lattice,false);
1721 : }
1722 169 : const IPosition latticeShape = lattice->shape();
1723 :
1724 169 : int samp=getAWConvFunc()->getOversampling();
1725 : //cerr << "SAMP " << samp << endl;
1726 : //Do sampling size correction
1727 169 : Vector<Float> sincConvX(nx);
1728 183977 : for (Int ix=0;ix<nx;ix++) {
1729 183808 : Float x=C::pi*Float(ix-nx/2)/(Float(nx)*Float(convSampling));
1730 183808 : if(ix==nx/2) {
1731 169 : sincConvX(ix)=1.0;
1732 : }
1733 : else {
1734 183639 : sincConvX(ix)=sin(x)/x;
1735 : }
1736 : }
1737 169 : Vector<Float> sincConvY(ny);
1738 183977 : for (Int ix=0;ix<ny;ix++) {
1739 183808 : Float x=C::pi*Float(ix-ny/2)/(Float(ny)*Float(convSampling));
1740 183808 : if(ix==ny/2) {
1741 169 : sincConvY(ix)=1.0;
1742 : }
1743 : else {
1744 183639 : sincConvY(ix)=sin(x)/x;
1745 : }
1746 : }
1747 :
1748 :
1749 : //cerr << convSampling << " max min of sincs " << max(sincConvX) << " " << min(sincConvX) << max(sincConvY) << " " << min(sincConvY) << endl;
1750 : //
1751 : // Now normalize the dirty image.
1752 : //
1753 : // Since *lattice is not copied to *image till the end of this
1754 : // method, normalizeImage also needs to work with Lattices
1755 : // (rather than ImageInterface).
1756 : //
1757 : // nx ny normalization from GridFT...
1758 : {
1759 169 : Int inx = lattice->shape()(0);
1760 169 : Int iny = lattice->shape()(1);
1761 169 : Vector<Complex> correction(inx);
1762 169 : correction=Complex(1.0, 0.0);
1763 : // Do the Grid-correction
1764 169 : IPosition cursorShape(4, inx, 1, 1, 1);
1765 169 : IPosition axisPath(4, 0, 1, 2, 3);
1766 169 : LatticeStepper lsx(lattice->shape(), cursorShape, axisPath);
1767 169 : LatticeIterator<Complex> lix(*lattice, lsx);
1768 295593 : for(lix.reset();!lix.atEnd();lix++)
1769 : {
1770 295424 : Int pol=lix.position()(2);
1771 295424 : Int chan=lix.position()(3);
1772 : {
1773 :
1774 295424 : Int iy=lix.position()(1);
1775 386433536 : for (Int ix=0;ix<nx;ix++) {
1776 386138112 : correction(ix)=1.0/(sincConvX(ix)*sincConvY(iy));
1777 : }
1778 : //cerr << iy << " min max corr " << min(abs(correction)) << " " << max(abs(correction)) << endl;
1779 295424 : lix.rwVectorCursor()*=correction;
1780 295424 : if(fftNormalization)
1781 : {
1782 0 : if(weights(pol,chan)!=0.0)
1783 : {
1784 0 : Complex rnorm(Float(inx)*Float(iny)/( weights(pol,chan) ));
1785 0 : lix.rwCursor()*=rnorm;
1786 : }
1787 : else
1788 0 : lix.woCursor()=0.0;
1789 : }
1790 : else
1791 : {
1792 295424 : Complex rnorm(Float(inx)*Float(iny));
1793 295424 : lix.rwCursor()*=rnorm;
1794 : }
1795 : }
1796 : }
1797 169 : }
1798 169 : if(!isTiled)
1799 : {
1800 : //
1801 : // Check the section from the image BEFORE converting to a lattice
1802 : //
1803 169 : LatticeLocker lock1 (*(image), FileLocker::Write);
1804 338 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
1805 338 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
1806 169 : IPosition stride(4, 1);
1807 338 : IPosition trc(blc+image->shape()-stride);
1808 : //
1809 : // Do the copy
1810 : //
1811 169 : image->put(griddedData(blc, trc));
1812 :
1813 :
1814 169 : if(!lattice.null()) lattice=0;
1815 169 : griddedData.resize(IPosition(1,0));
1816 169 : }
1817 169 : }
1818 :
1819 169 : return *image;
1820 169 : }
1821 : //
1822 : //---------------------------------------------------------------
1823 : //
1824 : // Get weight image
1825 76 : void AWProjectFT::getWeightImage(ImageInterface<Float>& weightImage,
1826 : Matrix<Float>& weights)
1827 : {
1828 76 : weights.resize(sumWeight.shape());
1829 76 : convertArray(weights,sumWeight);
1830 :
1831 76 : const IPosition latticeShape = weightImage.shape();
1832 76 : const IPosition avgpbShape = avgPB_p->shape();
1833 : //cout << "AWP::getWeightImage : weightimage shape : " << latticeShape << " and avgpb shape : " << avgpbShape << " nelems " << avgpbShape.nelements()<< " " << sumWeight << endl;
1834 76 : if(avgpbShape.nelements()==0 || ( avgpbShape != latticeShape) )
1835 0 : avgPB_p->resize(weightImage.shape());
1836 :
1837 76 : Int nx=latticeShape(0);
1838 76 : Int ny=latticeShape(1);
1839 :
1840 76 : int samp=getAWConvFunc()->getOversampling();
1841 : //Do sampling size correction
1842 76 : Vector<Float> sincConvX(nx);
1843 85068 : for (Int ix=0;ix<nx;ix++) {
1844 84992 : Float x=C::pi*Float(ix-nx/2)/(Float(nx)*Float(convSampling));
1845 84992 : if(ix==nx/2) {
1846 76 : sincConvX(ix)=1.0;
1847 : }
1848 : else {
1849 84916 : sincConvX(ix)=sin(x)/x;
1850 : }
1851 : }
1852 76 : Vector<Float> sincConvY(ny);
1853 85068 : for (Int ix=0;ix<ny;ix++) {
1854 84992 : Float x=C::pi*Float(ix-ny/2)/(Float(ny)*Float(convSampling));
1855 84992 : if(ix==ny/2) {
1856 76 : sincConvY(ix)=1.0;
1857 : }
1858 : else {
1859 84916 : sincConvY(ix)=sin(x)/x;
1860 : }
1861 : }
1862 :
1863 :
1864 :
1865 :
1866 : {
1867 76 : IPosition cursorShape(4, nx, ny, latticeShape(2), latticeShape(3));
1868 76 : IPosition axisPath(4, 0, 1, 2, 3);
1869 76 : LatticeStepper lsx(latticeShape, cursorShape, axisPath);
1870 76 : LatticeIterator<Float> lix(weightImage, lsx);
1871 76 : LatticeIterator<Float> liy(*avgPB_p,lsx);
1872 152 : for(lix.reset();!lix.atEnd();lix++)
1873 : {
1874 76 : lix.rwCursor()=liy.cursor();
1875 : }
1876 76 : }
1877 : {//sampling size correction
1878 76 : Vector<Float> correction(nx);
1879 76 : correction=1.0;
1880 : // Do the Grid-correction
1881 76 : IPosition cursorShape(4, nx, 1, 1, 1);
1882 76 : IPosition axisPath(4, 0, 1, 2, 3);
1883 76 : LatticeStepper lsx(weightImage.shape(), cursorShape, axisPath);
1884 76 : LatticeIterator<Float> lix(weightImage, lsx);
1885 132172 : for(lix.reset();!lix.atEnd();lix++)
1886 : {
1887 :
1888 132096 : Int iy=lix.position()(1);
1889 177865728 : for (Int ix=0;ix<nx;ix++) {
1890 177733632 : correction(ix)=1.0/(sincConvX(ix)*sincConvY(iy));
1891 : }
1892 132096 : lix.rwVectorCursor()*=correction;
1893 : }
1894 76 : }
1895 76 : }
1896 : //---------------------------------------------------------------
1897 144 : void AWProjectFT::setWeightImage(ImageInterface<Float>& weightImage){
1898 : //cerr <<"@@@loading weightimage" << endl;
1899 144 : IPosition latticeShape = weightImage.shape();
1900 144 : CoordinateSystem cs=weightImage.coordinates();
1901 144 : avgPB_p=new TempImage<Float>(latticeShape, cs);
1902 144 : avgPB_p->copyData(weightImage);
1903 144 : avgPBReady_p=True;
1904 :
1905 144 : }
1906 :
1907 : //---------------------------------------------------------------
1908 : //
1909 0 : Bool AWProjectFT::toRecord(RecordInterface& outRec, Bool withImage)
1910 : {
1911 : // Save the current AWProjectFT object to an output state record
1912 0 : Bool retval = true;
1913 0 : String error;
1914 : //save the base class variables
1915 0 : if(!FTMachine::toRecord(error, outRec, withImage, ""))
1916 0 : return false;
1917 0 : Double cacheVal=(Double) cachesize;
1918 0 : outRec.define("cache", cacheVal);
1919 0 : outRec.define("tile", tilesize);
1920 :
1921 0 : Vector<Double> phaseValue(2);
1922 0 : String phaseUnit;
1923 0 : phaseValue=mTangent_p.getAngle().getValue();
1924 0 : phaseUnit= mTangent_p.getAngle().getUnit();
1925 0 : outRec.define("phasevalue", phaseValue);
1926 0 : outRec.define("phaseunit", phaseUnit);
1927 :
1928 0 : Vector<Double> dirValue(3);
1929 0 : String dirUnit;
1930 0 : dirValue=mLocation_p.get("m").getValue();
1931 0 : dirUnit=mLocation_p.get("m").getUnit();
1932 0 : outRec.define("dirvalue", dirValue);
1933 0 : outRec.define("dirunit", dirUnit);
1934 :
1935 0 : outRec.define("padding", padding_p);
1936 0 : outRec.define("maxdataval", maxAbsData);
1937 :
1938 0 : Vector<Int> center_loc(4), offset_loc(4);
1939 0 : for (Int k=0; k<4 ; k++)
1940 : {
1941 0 : center_loc(k)=centerLoc(k);
1942 0 : offset_loc(k)=offsetLoc(k);
1943 : }
1944 0 : outRec.define("centerloc", center_loc);
1945 0 : outRec.define("offsetloc", offset_loc);
1946 0 : outRec.define("sumofweights", sumWeight);
1947 0 : outRec.define("sumofcfweights", sumCFWeight);
1948 0 : if(withImage && image)
1949 : {
1950 0 : ImageInterface<Complex>& tempimage(*image);
1951 0 : Record imageContainer;
1952 0 : String error;
1953 0 : retval = (retval || tempimage.toRecord(error, imageContainer));
1954 0 : outRec.defineRecord("image", imageContainer);
1955 0 : }
1956 0 : return retval;
1957 0 : }
1958 : //
1959 : //---------------------------------------------------------------
1960 : //
1961 0 : Bool AWProjectFT::fromRecord(const RecordInterface& inRec)
1962 : {
1963 0 : Bool retval = true;
1964 0 : String error;
1965 0 : if(!FTMachine::fromRecord(error, inRec))
1966 0 : return false;
1967 0 : imageCache=0; lattice=0;
1968 : Double cacheVal;
1969 0 : inRec.get("cache", cacheVal);
1970 0 : cachesize=(Long) cacheVal;
1971 0 : inRec.get("tile", tilesize);
1972 :
1973 0 : Vector<Double> phaseValue(2);
1974 0 : inRec.get("phasevalue",phaseValue);
1975 0 : String phaseUnit;
1976 0 : inRec.get("phaseunit",phaseUnit);
1977 0 : Quantity val1(phaseValue(0), phaseUnit);
1978 0 : Quantity val2(phaseValue(1), phaseUnit);
1979 0 : MDirection phasecenter(val1, val2);
1980 :
1981 0 : mTangent_p=phasecenter;
1982 : // This should be passed down too but the tangent plane is
1983 : // expected to be specified in all meaningful cases.
1984 0 : tangentSpecified_p=true;
1985 0 : Vector<Double> dirValue(3);
1986 0 : String dirUnit;
1987 0 : inRec.get("dirvalue", dirValue);
1988 0 : inRec.get("dirunit", dirUnit);
1989 0 : MVPosition dummyMVPos(dirValue(0), dirValue(1), dirValue(2));
1990 0 : MPosition mLocation(dummyMVPos, MPosition::ITRF);
1991 0 : mLocation_p=mLocation;
1992 :
1993 0 : inRec.get("padding", padding_p);
1994 0 : inRec.get("maxdataval", maxAbsData);
1995 :
1996 0 : Vector<Int> center_loc(4), offset_loc(4);
1997 0 : inRec.get("centerloc", center_loc);
1998 0 : inRec.get("offsetloc", offset_loc);
1999 0 : uInt ndim4 = 4;
2000 0 : centerLoc=IPosition(ndim4, center_loc(0), center_loc(1), center_loc(2),
2001 0 : center_loc(3));
2002 0 : offsetLoc=IPosition(ndim4, offset_loc(0), offset_loc(1), offset_loc(2),
2003 0 : offset_loc(3));
2004 0 : inRec.get("sumofweights", sumWeight);
2005 0 : inRec.get("sumofcfweights", sumCFWeight);
2006 0 : if(inRec.nfields() > 12 )
2007 : {
2008 0 : Record imageAsRec=inRec.asRecord("image");
2009 0 : if(!image) image= new TempImage<Complex>();
2010 :
2011 0 : String error;
2012 0 : retval = (retval || image->fromRecord(error, imageAsRec));
2013 :
2014 : // Might be changing the shape of sumWeight
2015 0 : init();
2016 :
2017 0 : if(isTiled)
2018 0 : lattice=CountedPtr<Lattice<Complex> > (image, false);
2019 : else
2020 : {
2021 : //
2022 : // Make the grid the correct shape and turn it into an
2023 : // array lattice Check the section from the image BEFORE
2024 : // converting to a lattice
2025 : //
2026 0 : IPosition gridShape(4, nx, ny, npol, nchan);
2027 0 : griddedData.resize(gridShape);
2028 0 : griddedData=Complex(0.0);
2029 0 : IPosition blc(4, (nx-image->shape()(0)+(nx%2==0))/2,
2030 0 : (ny-image->shape()(1)+(ny%2==0))/2, 0, 0);
2031 0 : IPosition start(4, 0);
2032 0 : IPosition stride(4, 1);
2033 0 : IPosition trc(blc+image->shape()-stride);
2034 0 : griddedData(blc, trc) = image->getSlice(start, image->shape());
2035 :
2036 0 : lattice=new ArrayLattice<Complex>(griddedData);
2037 0 : }
2038 :
2039 0 : AlwaysAssert(image, AipsError);
2040 0 : };
2041 0 : return retval;
2042 0 : }
2043 : //
2044 : //---------------------------------------------------------------
2045 : //
2046 27070 : void AWProjectFT::ok()
2047 : {
2048 27070 : AlwaysAssert(image, AipsError);
2049 27070 : }
2050 : //
2051 : //-------------------------------------------------------------------------
2052 : //
2053 93 : void AWProjectFT::setPAIncrement(const Quantity& computePAIncrement,
2054 : const Quantity& rotateOTFPAIncrement)
2055 : {
2056 186 : LogIO log_l(LogOrigin("AWProjectFT2", "setPAIncrement[R&D]"));
2057 :
2058 93 : rotateOTFPAIncr_p = rotateOTFPAIncrement.getValue("rad");
2059 93 : computePAIncr_p = computePAIncrement.getValue("rad");
2060 93 : convFuncCtor_p->setRotateCF(computePAIncr_p, rotateOTFPAIncr_p);
2061 :
2062 93 : paChangeDetector.setTolerance(computePAIncrement);
2063 93 : visResampler_p->setPATolerance(computePAIncrement.getValue("rad"));
2064 : log_l << LogIO::NORMAL <<"Setting PA increment to "
2065 93 : << computePAIncrement.getValue("deg") << " deg" << endl;
2066 93 : cfCache_p->setPAChangeDetector(paChangeDetector);
2067 93 : }
2068 : //
2069 : //-------------------------------------------------------------------------
2070 : //
2071 0 : Bool AWProjectFT::verifyShapes(IPosition pbShape, IPosition skyShape)
2072 : {
2073 0 : LogIO log_l(LogOrigin("AWProjectFT2", "verifyShapes[R&D]"));
2074 :
2075 0 : if ((pbShape(0) != skyShape(0)) && // X-axis
2076 0 : (pbShape(1) != skyShape(1)) && // Y-axis
2077 0 : (pbShape(2) != skyShape(2))) // Poln-axis
2078 : {
2079 : log_l << "Sky and/or polarization shape of the avgPB and the sky model do not match."
2080 0 : << LogIO::EXCEPTION;
2081 0 : return false;
2082 : }
2083 0 : return true;
2084 :
2085 0 : }
2086 : //
2087 : //-------------------------------------------------------------------------
2088 : //
2089 27035 : void AWProjectFT::setupVBStore(VBStore& vbs,
2090 : const VisBuffer2& vb,
2091 : const Matrix<Float>& imagingweight,
2092 : const Cube<Complex>& visData,
2093 : const Matrix<Double>& uvw,
2094 : const Cube<Int>& flagCube,
2095 : const Vector<Double>& dphase,
2096 : const Bool& dopsf,
2097 : const Vector<Int>& /*gridShape*/)
2098 : {
2099 27035 : vbs.vb_p = &vb;
2100 27035 : vbs.wbAWP_p=wbAWP_p;
2101 27035 : vbs.ftmType_p=ftmType_p;
2102 27035 : vbs.nWPlanes_p = nWPlanes_p;
2103 27035 : makeCFPolMap(vb,cfStokes_p,CFMap_p);
2104 27035 : makeConjPolMap(vb,CFMap_p,ConjCFMap_p);
2105 :
2106 27035 : visResampler_p->setParams(uvScale,uvOffset,dphase);
2107 27035 : visResampler_p->setMaps(chanMap, polMap);
2108 27035 : visResampler_p->setCFMaps(CFMap_p, ConjCFMap_p);
2109 27035 : visResampler_p->setFreqMaps(expandedSpwFreqSel_p,expandedSpwConjFreqSel_p);
2110 : //
2111 : // Set up VBStore object to point to the relavent info. of the VB.
2112 : //
2113 27035 : vbs.imRefFreq_p = imRefFreq_p;
2114 27035 : vbs.nRow_p = vb.nRows();
2115 27035 : vbs.beginRow_p = 0;
2116 27035 : vbs.endRow_p = vbs.nRow_p;
2117 27035 : vbs.spwID_p = vb.spectralWindows()(0);
2118 27035 : vbs.nDataPol_p = flagCube.shape()[0];
2119 27035 : vbs.nDataChan_p = flagCube.shape()[1];
2120 :
2121 27035 : vbs.antenna1_p.reference(vb.antenna1());
2122 27035 : vbs.antenna2_p.reference(vb.antenna2());
2123 27035 : vbs.paQuant_p = Quantity(getPA(vb),"rad");
2124 :
2125 27035 : vbs.corrType_p.reference(vb.correlationTypes());
2126 :
2127 27035 : vbs.uvw_p=uvw;
2128 27035 : vbs.imagingWeight_p.reference(imagingweight);
2129 27035 : vbs.visCube_p.reference(visData);
2130 :
2131 27035 : vbs.freq_p.reference(vb.getFrequencies(0));
2132 :
2133 27035 : vbs.rowFlag_p.reference(vb.flagRow());
2134 27035 : if(!usezero_p)
2135 0 : for (Int rownr=0; rownr<vbs.nRow_p; rownr++)
2136 0 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) vbs.rowFlag_p(rownr)=true;
2137 :
2138 27035 : vbs.flagCube_p.resize(flagCube.shape()); vbs.flagCube_p = false; vbs.flagCube_p(flagCube!=0) = true;
2139 :
2140 27035 : vbs.conjBeams_p=conjBeams_p;
2141 :
2142 : //timer_p.mark();
2143 :
2144 27035 : po_p->fetchPointingOffset(*image, vb, doPointing);
2145 :
2146 27035 : if (makingPSF || (vbs.ftmType_p==casa::refim::FTMachine::WEIGHT)){
2147 8446 : cfwts2_p->invokeGC(vbs.spwID_p);
2148 8446 : vb2CFBMap_p->setDoPointing(doPointing);
2149 16892 : vb2CFBMap_p->makeVBRow2CFBMap(*cfwts2_p,
2150 : vb,
2151 8446 : paChangeDetector.getParAngleTolerance(),
2152 8446 : chanMap,polMap,po_p);
2153 : }
2154 : else
2155 : {
2156 : // If the Wt. CFs are still in the memory, clear them. They
2157 : // won't be required again (though with the silly check below,
2158 : // if the in-memory Wt. CFs are less than 1KB, they will be
2159 : // left in memory).
2160 18589 : if (cfwts2_p->memUsage() > 1000)
2161 : {
2162 21 : cfwts2_p->clear();
2163 : }
2164 :
2165 18589 : cfs2_p->invokeGC(vbs.spwID_p);
2166 18589 : vb2CFBMap_p->setDoPointing(doPointing);
2167 37178 : vb2CFBMap_p->makeVBRow2CFBMap(*cfs2_p, vb,
2168 18589 : paChangeDetector.getParAngleTolerance(),
2169 18589 : chanMap,polMap,po_p);
2170 : }
2171 :
2172 : //
2173 : // Trigger the computation of phase gradiant corresponding to the
2174 : // field offset (from the VB) w.r.t. the image phase center.
2175 : //
2176 : //
2177 : // For AzElApertures, this rotates the CFs.
2178 : //
2179 27035 : convFuncCtor_p->prepareConvFunction(vb,*vb2CFBMap_p);
2180 :
2181 27035 : vbs.accumCFs_p=((vbs.uvw_p.nelements() == 0) && dopsf);
2182 27035 : visResampler_p->setVB2CFMap(vb2CFBMap_p);
2183 :
2184 : // The following code is required only for GPU or multi-threaded
2185 : //gridder. Currently does not work without the rest of the
2186 : //GPU/multi-threaded infrastructure (though, I (SB) thought this
2187 : //was designed to be benign for normal gridding).
2188 : //
2189 27035 : visResampler_p->initializeDataBuffers(vbs);
2190 27035 : }
2191 : //
2192 : //---------------------------------------------------------------
2193 : //
2194 0 : void AWProjectFT::get(VisBuffer2& vb, Cube<Complex>& modelVis,
2195 : Array<Complex>& griddedVis, Vector<Double>& scale,
2196 : Int row)
2197 : {
2198 :
2199 : (void)scale; //Suppress the warning
2200 :
2201 0 : Int nX=griddedVis.shape()(0);
2202 0 : Int nY=griddedVis.shape()(1);
2203 0 : Vector<Double> offset(2);
2204 0 : offset(0)=Double(nX)/2.0;
2205 0 : offset(1)=Double(nY)/2.0;
2206 : // If row is -1 then we pass through all rows
2207 : Int startRow, endRow, nRow;
2208 0 : if (row==-1)
2209 : {
2210 0 : nRow=vb.nRows();
2211 0 : startRow=0;
2212 0 : endRow=nRow-1;
2213 0 : modelVis.set(Complex(0.0,0.0));
2214 : }
2215 : else
2216 : {
2217 0 : nRow=1;
2218 0 : startRow=row;
2219 0 : endRow=row;
2220 : }
2221 :
2222 0 : Int NAnt=0;
2223 :
2224 0 : if (doPointing)
2225 0 : NAnt = findPointingOffsets(vb,l_offsets, m_offsets,true);
2226 :
2227 :
2228 : //
2229 : // Get the uvws in a form that Fortran can use
2230 : //
2231 0 : Matrix<Double> uvw(negateUV(vb));
2232 0 : Vector<Double> dphase(vb.nRows());
2233 0 : dphase=0.0;
2234 0 : doUVWRotation_p=true;
2235 : //rotateUVW(uvw, dphase, vb);
2236 0 : girarUVW(uvw, dphase, vb);
2237 0 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
2238 :
2239 : // This is the convention for dphase
2240 : //dphase*=-1.0;
2241 :
2242 0 : Cube<Int> flags(vb.flagCube().shape());
2243 0 : flags=0;
2244 0 : flags(vb.flagCube())=true;
2245 :
2246 :
2247 0 : matchChannel(vb);
2248 0 : Vector<Int> rowFlags(vb.nRows());
2249 0 : rowFlags=0;
2250 0 : rowFlags(vb.flagRow())=true;
2251 0 : if(!usezero_p)
2252 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
2253 0 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
2254 :
2255 0 : visResampler_p->setParams(uvScale,uvOffset,dphase);
2256 0 : visResampler_p->setMaps(chanMap, polMap);
2257 0 : visResampler_p->setFreqMaps(expandedSpwFreqSel_p,expandedSpwConjFreqSel_p);
2258 :
2259 0 : IPosition s(modelVis.shape());
2260 0 : Int Conj=0,doGrad=0,ScanNo=0;
2261 0 : Double area=1.0;
2262 0 : Int tmpPAI=1;
2263 0 : Cube<Complex> visCubeModel(vb.visCubeModel());
2264 0 : runFortranGet(uvw,dphase,visCubeModel,s,Conj,flags,rowFlags,row,
2265 0 : offset,&griddedVis,nx,ny,npol,nchan,vb,NAnt,ScanNo,sigma,
2266 0 : l_offsets,m_offsets,area,doGrad,tmpPAI);
2267 0 : vb.setVisCubeModel(visCubeModel);
2268 0 : }
2269 : //
2270 : //----------------------------------------------------------------------
2271 : //
2272 0 : void AWProjectFT::initVisBuffer(VisBuffer2& vb, Type whichVBColumn)
2273 : {
2274 0 : if (whichVBColumn == FTMachine::MODEL) vb.setVisCubeModel(Complex(0.0,0.0));
2275 0 : else if (whichVBColumn == FTMachine::OBSERVED) vb.setVisCube(Complex(0.0,0.0));
2276 0 : }
2277 : //
2278 : //----------------------------------------------------------------------
2279 : //
2280 0 : void AWProjectFT::initVisBuffer(VisBuffer2&, // vb
2281 : Type, // whichVBColumn
2282 : Int // row
2283 : )
2284 : {
2285 0 : cerr << "AWProjectFT::initVisBuffer(VisBuffer2& vb, Type whichVBColumn, Int row) disabled" << endl;
2286 0 : }
2287 :
2288 0 : void AWProjectFT::ComputeResiduals(VisBuffer2&vb, Bool useCorrected)
2289 : {
2290 0 : VBStore vbs;
2291 :
2292 0 : vbs.nRow_p = vb.nRows();
2293 0 : vbs.beginRow_p = 0;
2294 0 : vbs.endRow_p = vbs.nRow_p;
2295 :
2296 0 : vbs.modelCube_p.reference(vb.visCubeModel());
2297 0 : if (useCorrected) vbs.correctedCube_p.reference(vb.visCubeCorrected());
2298 0 : else vbs.visCube_p.reference(vb.visCube());
2299 0 : vbs.useCorrected_p = useCorrected;
2300 0 : visResampler_p->ComputeResiduals(vbs);
2301 0 : }
2302 :
2303 : //
2304 : //---------------------------------------------------------------
2305 : //---------------------------------------------------------------
2306 : //---------------------------------------------------------------
2307 : //---------------------------------------------------------------
2308 : // THIS IS FOR NON-PRODUCTION WORK. HERE SINCE IT DOES LINK TO POINTING SELFCAL CODE
2309 : //
2310 : // Predict the coherences as well as their derivatives w.r.t. the
2311 : // pointing offsets.
2312 : //
2313 0 : void AWProjectFT::nget(VisBuffer2& vb,
2314 : // These offsets should be appropriate for the VB
2315 : Array<Float>& l_off, Array<Float>& m_off,
2316 : Cube<Complex>& Mout,
2317 : Cube<Complex>& dMout1,
2318 : Cube<Complex>& dMout2,
2319 : Int Conj, Int doGrad)
2320 : {
2321 0 : LogIO log_l(LogOrigin("AWProjectFT2", "nget[R&D]"));
2322 : Int startRow, endRow, nRow;
2323 0 : nRow=vb.nRows();
2324 0 : startRow=0;
2325 0 : endRow=nRow-1;
2326 :
2327 0 : Mout = dMout1 = dMout2 = Complex(0,0);
2328 :
2329 0 : findConvFunction(*image, vb);
2330 0 : Int NAnt=0;
2331 0 : Nant_p = vb.subtableColumns().antenna().nrow();
2332 0 : if (doPointing)
2333 0 : NAnt = findPointingOffsets(vb,l_offsets,m_offsets,false);
2334 :
2335 0 : l_offsets=l_off;
2336 0 : m_offsets=m_off;
2337 0 : Matrix<Double> uvw(negateUV(vb));
2338 :
2339 0 : Vector<Double> dphase(vb.nRows());
2340 0 : dphase=0.0;
2341 0 : doUVWRotation_p=true;
2342 : //rotateUVW(uvw, dphase, vb);
2343 0 : girarUVW(uvw, dphase, vb);
2344 0 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
2345 : // This is the convention for dphase
2346 : // dphase*=-1.0;
2347 :
2348 0 : Cube<Int> flags(vb.flagCube().shape());
2349 0 : flags=0;
2350 0 : flags(vb.flagCube())=true;
2351 :
2352 0 : Vector<Int> rowFlags(vb.nRows());
2353 0 : rowFlags=0;
2354 0 : rowFlags(vb.flagRow())=true;
2355 0 : if(!usezero_p)
2356 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
2357 0 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
2358 :
2359 0 : IPosition s,gradS;
2360 0 : Cube<Complex> visdata,gradVisAzData,gradVisElData;
2361 : //
2362 : // visdata now references the Mout data structure rather than to the internal VB storeage.
2363 : //
2364 0 : visdata.reference(Mout);
2365 :
2366 0 : if (doGrad)
2367 : {
2368 : // The following should reference some slice of dMout?
2369 0 : gradVisAzData.reference(dMout1);
2370 0 : gradVisElData.reference(dMout2);
2371 : }
2372 0 : visResampler_p->setParams(uvScale,uvOffset,dphase);
2373 0 : visResampler_p->setMaps(chanMap, polMap);
2374 0 : visResampler_p->setFreqMaps(expandedSpwFreqSel_p,expandedSpwConjFreqSel_p);
2375 :
2376 0 : makeCFPolMap(vb,cfStokes_p,CFMap_p);
2377 0 : makeConjPolMap(vb,CFMap_p,ConjCFMap_p);
2378 0 : visResampler_p->setCFMaps(CFMap_p, ConjCFMap_p);
2379 : //
2380 : // Begin the actual de-gridding.
2381 : //
2382 0 : if(isTiled)
2383 : {
2384 0 : log_l << "The sky model is tiled" << LogIO::NORMAL << LogIO::POST;
2385 0 : Double invLambdaC=vb.getFrequencies(0)(0)/C::c;
2386 0 : Vector<Double> uvLambda(2);
2387 0 : Vector<Int> centerLoc2D(2);
2388 0 : centerLoc2D=0;
2389 :
2390 : // Loop over all rows
2391 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
2392 : {
2393 :
2394 : // Calculate uvw for this row at the center frequency
2395 0 : uvLambda(0)=uvw(0, rownr)*invLambdaC;
2396 0 : uvLambda(1)=uvw(1, rownr)*invLambdaC;
2397 0 : centerLoc2D=gridder->location(centerLoc2D, uvLambda);
2398 :
2399 : // Is this point on the grid?
2400 0 : if(gridder->onGrid(centerLoc2D))
2401 : {
2402 :
2403 : // Get the tile
2404 0 : Array<Complex>* dataPtr=getDataPointer(centerLoc2D, true);
2405 0 : gridder->setOffset(IPosition(2, offsetLoc(0), offsetLoc(1)));
2406 0 : Int aNx=dataPtr->shape()(0);
2407 0 : Int aNy=dataPtr->shape()(1);
2408 :
2409 : // Now use FORTRAN to do the gridding. Remember to
2410 : // ensure that the shape and offsets of the tile are
2411 : // accounted for.
2412 :
2413 0 : Vector<Double> actualOffset(3);
2414 0 : for (Int i=0;i<2;i++)
2415 0 : actualOffset(i)=uvOffset(i)-Double(offsetLoc(i));
2416 :
2417 0 : actualOffset(2)=uvOffset(2);
2418 0 : IPosition s(vb.visCubeModel().shape());
2419 :
2420 0 : Int ScanNo=0, tmpPAI;
2421 0 : Double area=1.0;
2422 0 : tmpPAI = 1;
2423 0 : runFortranGetGrad(uvw,dphase,visdata,s,
2424 : gradVisAzData,gradVisElData,
2425 : Conj,flags,rowFlags,rownr,
2426 0 : actualOffset,dataPtr,aNx,aNy,npol,nchan,vb,NAnt,ScanNo,sigma,
2427 0 : l_offsets,m_offsets,area,doGrad,tmpPAI);
2428 0 : }
2429 : }
2430 0 : }
2431 : else
2432 : {
2433 0 : IPosition s(vb.visCubeModel().shape());
2434 0 : Int ScanNo=0, tmpPAI, trow=-1;
2435 0 : Double area=1.0;
2436 0 : tmpPAI = 1;
2437 0 : runFortranGetGrad(uvw,dphase,visdata/*vb.modelVisCube()*/,s,
2438 : gradVisAzData, gradVisElData,
2439 : Conj,flags,rowFlags,trow,
2440 0 : uvOffset,&griddedData,nx,ny,npol,nchan,vb,NAnt,ScanNo,sigma,
2441 0 : l_offsets,m_offsets,area,doGrad,tmpPAI);
2442 0 : }
2443 :
2444 0 : }
2445 0 : void AWProjectFT::get(VisBuffer2& vb,
2446 : VisBuffer2& gradVBAz,
2447 : VisBuffer2& gradVBEl,
2448 : Cube<Float>& pointingOffsets,
2449 : Int row, // default row=-1
2450 : Type whichVBColumn, // default whichVBColumn = FTMachine::MODEL
2451 : Type whichGradVBColumn,// default whichGradVBColumn = FTMachine::MODEL
2452 : Int Conj, Int doGrad) // default Conj=0, doGrad=1
2453 : {
2454 : // If row is -1 then we pass through all rows
2455 : Int startRow, endRow, nRow;
2456 0 : if (row==-1)
2457 : {
2458 0 : nRow=vb.nRows();
2459 0 : startRow=0;
2460 0 : endRow=nRow-1;
2461 0 : initVisBuffer(vb,whichVBColumn);
2462 0 : if (doGrad)
2463 : {
2464 0 : initVisBuffer(gradVBAz, whichGradVBColumn);
2465 0 : initVisBuffer(gradVBEl, whichGradVBColumn);
2466 : }
2467 : }
2468 : else
2469 : {
2470 0 : nRow=1;
2471 0 : startRow=row;
2472 0 : endRow=row;
2473 0 : initVisBuffer(vb,whichVBColumn,row);
2474 0 : if (doGrad)
2475 : {
2476 0 : initVisBuffer(gradVBAz, whichGradVBColumn,row);
2477 0 : initVisBuffer(gradVBEl, whichGradVBColumn,row);
2478 : }
2479 : }
2480 :
2481 0 : findConvFunction(*image, vb);
2482 :
2483 0 : Nant_p = vb.subtableColumns().antenna().nrow();
2484 0 : Int NAnt=0;
2485 0 : if (doPointing)
2486 0 : NAnt = findPointingOffsets(vb,pointingOffsets,l_offsets,m_offsets,false);
2487 :
2488 0 : Matrix<Double> uvw(negateUV(vb));
2489 :
2490 0 : Vector<Double> dphase(vb.nRows());
2491 0 : dphase=0.0;
2492 0 : doUVWRotation_p=true;
2493 : // rotateUVW(uvw, dphase, vb);
2494 0 : girarUVW(uvw, dphase, vb);
2495 0 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
2496 :
2497 : // This is the convention for dphase
2498 : //dphase*=-1.0;
2499 :
2500 :
2501 0 : Cube<Int> flags(vb.flagCube().shape());
2502 0 : flags=0;
2503 0 : flags(vb.flagCube())=true;
2504 : //
2505 : //Check if ms has changed then cache new spw and chan selection
2506 : //
2507 : // if(vb.newMS()) matchAllSpwChans(vb);
2508 :
2509 : //Here we redo the match or use previous match
2510 : //
2511 : //Channel matching for the actual spectral window of buffer
2512 : //
2513 : // if(doConversion_p[vb.spectralWindow()])
2514 : // matchChannel(vb.spectralWindow(), vb);
2515 : // else
2516 : // {
2517 : // chanMap.resize();
2518 : // chanMap=multiChanMap_p[vb.spectralWindow()];
2519 : // }
2520 :
2521 0 : matchChannel(vb);
2522 0 : Vector<Int> rowFlags(vb.nRows());
2523 0 : rowFlags=0;
2524 0 : rowFlags(vb.flagRow())=true;
2525 0 : if(!usezero_p)
2526 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
2527 0 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
2528 :
2529 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
2530 0 : if (vb.antenna1()(rownr) != vb.antenna2()(rownr))
2531 0 : rowFlags(rownr) = (vb.flagRow()(rownr)==true);
2532 :
2533 0 : IPosition s,gradS;
2534 0 : Cube<Complex> visdata,gradVisAzData,gradVisElData;
2535 0 : if (whichVBColumn == FTMachine::MODEL)
2536 : {
2537 0 : s = vb.visCubeModel().shape();
2538 0 : visdata.reference(vb.visCubeModel());
2539 : }
2540 0 : else if (whichVBColumn == FTMachine::OBSERVED)
2541 : {
2542 0 : s = vb.visCube().shape();
2543 0 : visdata.reference(vb.visCube());
2544 : }
2545 :
2546 0 : if (doGrad)
2547 : {
2548 0 : if (whichGradVBColumn == FTMachine::MODEL)
2549 : {
2550 : // gradS = gradVBAz.modelVisCube().shape();
2551 0 : gradVisAzData.reference(gradVBAz.visCubeModel());
2552 0 : gradVisElData.reference(gradVBEl.visCubeModel());
2553 : }
2554 0 : else if (whichGradVBColumn == FTMachine::OBSERVED)
2555 : {
2556 : // gradS = gradVBAz.visCube().shape();
2557 0 : gradVisAzData.reference(gradVBAz.visCube());
2558 0 : gradVisElData.reference(gradVBEl.visCube());
2559 : }
2560 : }
2561 0 : visResampler_p->setParams(uvScale,uvOffset,dphase);
2562 0 : visResampler_p->setFreqMaps(expandedSpwFreqSel_p,expandedSpwConjFreqSel_p);
2563 0 : visResampler_p->setMaps(chanMap, polMap);
2564 : // Vector<Int> ConjCFMap, CFMap;
2565 0 : makeCFPolMap(vb,cfStokes_p,CFMap_p);
2566 0 : makeConjPolMap(vb,CFMap_p,ConjCFMap_p);
2567 0 : visResampler_p->setCFMaps(CFMap_p, ConjCFMap_p);
2568 :
2569 0 : if(isTiled)
2570 : {
2571 0 : Double invLambdaC=vb.getFrequencies(0)(0)/C::c;
2572 0 : Vector<Double> uvLambda(2);
2573 0 : Vector<Int> centerLoc2D(2);
2574 0 : centerLoc2D=0;
2575 :
2576 : // Loop over all rows
2577 0 : for (Int rownr=startRow; rownr<=endRow; rownr++)
2578 : {
2579 :
2580 : // Calculate uvw for this row at the center frequency
2581 0 : uvLambda(0)=uvw(0, rownr)*invLambdaC;
2582 0 : uvLambda(1)=uvw(1, rownr)*invLambdaC;
2583 0 : centerLoc2D=gridder->location(centerLoc2D, uvLambda);
2584 :
2585 : // Is this point on the grid?
2586 0 : if(gridder->onGrid(centerLoc2D))
2587 : {
2588 :
2589 : // Get the tile
2590 0 : Array<Complex>* dataPtr=getDataPointer(centerLoc2D, true);
2591 0 : gridder->setOffset(IPosition(2, offsetLoc(0), offsetLoc(1)));
2592 0 : Int aNx=dataPtr->shape()(0);
2593 0 : Int aNy=dataPtr->shape()(1);
2594 :
2595 : // Now use FORTRAN to do the gridding. Remember to
2596 : // ensure that the shape and offsets of the tile are
2597 : // accounted for.
2598 :
2599 0 : Vector<Double> actualOffset(3);
2600 0 : for (Int i=0;i<2;i++)
2601 0 : actualOffset(i)=uvOffset(i)-Double(offsetLoc(i));
2602 :
2603 0 : actualOffset(2)=uvOffset(2);
2604 0 : IPosition s(vb.visCubeModel().shape());
2605 :
2606 0 : Int ScanNo=0, tmpPAI;
2607 0 : Double area=1.0;
2608 0 : tmpPAI = 1;
2609 0 : runFortranGetGrad(uvw,dphase,visdata,s,
2610 : gradVisAzData,gradVisElData,
2611 : Conj,flags,rowFlags,rownr,
2612 0 : actualOffset,dataPtr,aNx,aNy,npol,nchan,vb,NAnt,ScanNo,sigma,
2613 0 : l_offsets,m_offsets,area,doGrad,tmpPAI);
2614 0 : }
2615 : }
2616 0 : }
2617 : else
2618 : {
2619 :
2620 0 : IPosition s(vb.visCubeModel().shape());
2621 0 : Int ScanNo=0, tmpPAI;
2622 0 : Double area=1.0;
2623 :
2624 0 : tmpPAI = 1;
2625 :
2626 0 : runFortranGetGrad(uvw,dphase,visdata/*vb.modelVisCube()*/,s,
2627 : gradVisAzData, gradVisElData,
2628 : Conj,flags,rowFlags,row,
2629 0 : uvOffset,&griddedData,nx,ny,npol,nchan,vb,NAnt,ScanNo,sigma,
2630 0 : l_offsets,m_offsets,area,doGrad,tmpPAI);
2631 0 : }
2632 0 : }
2633 : //---------------------------------------------------------------
2634 : //---------------------------------------------------------------
2635 : //---------------------------------------------------------------
2636 : //---------------------------------------------------------------
2637 : // THIS IS FOR NON-PRODUCTION WORK. HERE SINCE IT DOES LINK TO POINTING SELFCAL CODE
2638 :
2639 : } //# NAMESPACE CASA - END
2640 : };
|
