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