Line data Source code
1 : // -*- C++ -*-
2 : //# AWConvFunc.cc: Implementation of the AWConvFunc 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 :
30 :
31 : #include <synthesis/TransformMachines2/AWConvFunc.h>
32 : #include <synthesis/TransformMachines2/AWProjectFT.h>
33 : #include <synthesis/TransformMachines/SynthesisError.h>
34 : #include <casacore/images/Images/ImageInterface.h>
35 : #include <synthesis/TransformMachines2/Utils.h>
36 : #include <synthesis/TransformMachines/BeamCalc.h>
37 : #include <synthesis/TransformMachines2/CFStore.h>
38 : #include <synthesis/TransformMachines2/CFStore2.h>
39 : #include <synthesis/TransformMachines2/VB2CFBMap.h>
40 : #include <synthesis/TransformMachines2/PSTerm.h>
41 : #include <synthesis/TransformMachines2/WTerm.h>
42 : #include <synthesis/TransformMachines2/ATerm.h>
43 : #include <synthesis/TransformMachines2/VLACalcIlluminationConvFunc.h>
44 : #include <synthesis/TransformMachines2/ConvolutionFunction.h>
45 : #include <synthesis/TransformMachines2/PolOuterProduct.h>
46 :
47 : #include <synthesis/TransformMachines2/EVLAAperture.h>
48 : #include <synthesis/TransformMachines2/WPConvFunc.h>
49 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
50 : #include <casacore/coordinates/Coordinates/LinearCoordinate.h>
51 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
52 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
53 : #include <casacore/casa/System/ProgressMeter.h>
54 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
55 : #include <casacore/casa/Utilities/CompositeNumber.h>
56 : #include <casacore/casa/OS/Directory.h>
57 : #include <casacore/casa/OS/Timer.h>
58 : #include <ostream>
59 : #ifdef _OPENMP
60 : #include <omp.h>
61 : #endif
62 :
63 : #define MAX_FREQ 1e30
64 :
65 : using namespace casacore;
66 : namespace casa{
67 : namespace refim{
68 : using namespace casacore;
69 : using namespace casa::vi;
70 :
71 0 : AWConvFunc::AWConvFunc(const casacore::CountedPtr<ATerm> aTerm,
72 : const casacore::CountedPtr<PSTerm> psTerm,
73 : const casacore::CountedPtr<WTerm> wTerm,
74 : const casacore::Bool wbAWP,
75 0 : const casacore::Bool conjPB):
76 0 : ConvolutionFunction(),aTerm_p(aTerm),psTerm_p(psTerm), wTerm_p(wTerm), pixFieldGrad_p(),
77 0 : wbAWP_p(wbAWP), conjPB_p(conjPB), baseCFB_p(), atermMaker_p(nullptr), wtermMaker_p(nullptr)
78 : {
79 0 : if (psTerm->isNoOp() && aTerm->isNoOp())
80 : {
81 0 : LogIO log_l(LogOrigin("AWConvFunc", "AWConvFunc"));
82 0 : log_l << "Both, psterm and aterm cannot be set to NoOp. " << LogIO::EXCEPTION;
83 0 : }
84 0 : if (wbAWP && aTerm->isNoOp())
85 : {
86 : //log_l << "wbawp=True is ineffective when aterm is OFF. Setting wbawp to False." << LogIO::NORMAL1;
87 0 : wbAWP_p=false;
88 : }
89 :
90 0 : pixFieldGrad_p.resize(2);pixFieldGrad_p(0)=0.0; pixFieldGrad_p(1)=0.0;
91 0 : }
92 :
93 0 : AWConvFunc::AWConvFunc(std::shared_ptr<EVLAAperture> aterm, std::shared_ptr<WPConvFunc> wterm){
94 0 : atermMaker_p=aterm;
95 0 : wtermMaker_p=wterm;
96 :
97 :
98 0 : }
99 :
100 :
101 : //
102 : //----------------------------------------------------------------------
103 : //
104 0 : AWConvFunc& AWConvFunc::operator=(const AWConvFunc& other)
105 : {
106 0 : if(this!=&other)
107 : {
108 0 : aTerm_p = other.aTerm_p;
109 0 : psTerm_p = other.psTerm_p;
110 0 : wTerm_p = other.wTerm_p;
111 0 : atermMaker_p=other.atermMaker_p;
112 0 : wtermMaker_p=other.wtermMaker_p;
113 :
114 : }
115 0 : return *this;
116 :
117 : }
118 : //
119 : //----------------------------------------------------------------------
120 : //
121 0 : void AWConvFunc::makePBSq(ImageInterface<Complex>& PB)
122 : {
123 0 : IPosition pbShape=PB.shape();
124 0 : IPosition cursorShape(4, pbShape(0), pbShape(1), 1, 1), axisPath(4,0,1,2,3);
125 0 : Array<Complex> buf; PB.get(buf,false);
126 0 : ArrayLattice<Complex> lat(buf, true);
127 0 : LatticeStepper latStepper(lat.shape(), cursorShape,axisPath);
128 0 : LatticeIterator<Complex> latIter(lat, latStepper);
129 :
130 0 : IPosition start0(4,0,0,0,0), start1(4,0,0,1,0), length(4, pbShape(0), pbShape(1),1,1);
131 0 : Slicer slicePol0(start0, length), slicePol1(start1, length);
132 0 : if (pbShape(2) > 1)
133 : {
134 0 : Array<Complex> pol0, pol1,tmp;
135 :
136 0 : lat.getSlice(pol0, slicePol0);
137 0 : lat.getSlice(pol1, slicePol1);
138 0 : tmp = pol0;
139 0 : pol0 = pol0*conj(pol1);
140 0 : pol1 = tmp*conj(pol1);
141 0 : lat.putSlice(pol0,start0);
142 0 : lat.putSlice(pol1,start1);
143 0 : }
144 : else
145 : {
146 : // Array<Complex> pol0;
147 : // lat.getSlice(pol0,slicePol0);
148 : // pol0 = pol0*conj(pol0);
149 0 : buf = buf * conj(buf);
150 : }
151 0 : }
152 : //
153 : //----------------------------------------------------------------------
154 : //
155 0 : void AWConvFunc::makeConjPolAxis(CoordinateSystem& cs,
156 : Int conjStokes_in)
157 : {
158 : //LogIO log_l(LogOrigin("AWConvFunc2", "makeConjPolAxis[R&D]"));
159 0 : IPosition dummy;
160 0 : Vector<String> csList;
161 0 : Vector<Int> stokes, conjStokes;
162 :
163 : // cout << "CoordSys: ";
164 : // csList = cs.list(log_l,MDoppler::RADIO,dummy,dummy);
165 : // cout << csList << endl;
166 0 : Int stokesIndex=cs.findCoordinate(Coordinate::STOKES);
167 0 : StokesCoordinate sc=cs.stokesCoordinate(stokesIndex);
168 :
169 0 : if (conjStokes_in == -1)
170 : {
171 0 : stokes=sc.stokes();
172 0 : conjStokes.resize(stokes.shape());
173 0 : for (uInt i=0; i<stokes.nelements(); i++)
174 : {
175 0 : if (stokes(i) == Stokes::RR) conjStokes(i) = Stokes::LL;
176 0 : if (stokes(i) == Stokes::LL) conjStokes(i) = Stokes::RR;
177 0 : if (stokes(i) == Stokes::LR) conjStokes(i) = Stokes::RL;
178 0 : if (stokes(i) == Stokes::RL) conjStokes(i) = Stokes::LR;
179 :
180 0 : if (stokes(i) == Stokes::XX) conjStokes(i) = Stokes::YY;
181 0 : if (stokes(i) == Stokes::YY) conjStokes(i) = Stokes::XX;
182 0 : if (stokes(i) == Stokes::YX) conjStokes(i) = Stokes::XY;
183 0 : if (stokes(i) == Stokes::XY) conjStokes(i) = Stokes::YX;
184 : }
185 : }
186 : else
187 : {
188 0 : conjStokes.resize(1);
189 0 : conjStokes[0]=conjStokes_in;
190 : }
191 0 : sc.setStokes(conjStokes);
192 0 : cs.replaceCoordinate(sc,stokesIndex);
193 0 : }
194 : //
195 : //----------------------------------------------------------------------
196 : //
197 0 : void AWConvFunc::fillConvFuncBuffer(CFBuffer& cfb, CFBuffer& cfWtb,
198 : const Int&,// skyNX,
199 : const Int&,// skyNY,
200 : const Vector<Double>&,// skyIncr,
201 : const Int& nx, const Int& ny,
202 : const Vector<Double>& freqValues,
203 : const Vector<Double>& wValues,
204 : const Double& wScale,
205 : const Double& vbPA, const Double& freqHi,
206 : const PolMapType& muellerElements,
207 : const PolMapType& muellerElementsIndex,
208 : const VisBuffer2& vb,
209 : const Float& psScale,
210 : PSTerm& psTerm, WTerm& wTerm, ATerm& aTerm,
211 : Bool isDryRun)
212 : {
213 : // Unused variable from the dark-ages era interface that should ultimately go.
214 : (void)psScale;
215 : (void)muellerElementsIndex;
216 : (void)freqHi;
217 : // Int ttt=0;
218 0 : Complex cfNorm, cfWtNorm;
219 : //Double vbPA = getPA(vb);
220 0 : Complex cpeak,wtcpeak;
221 0 : aTerm.cacheVBInfo(vb);
222 0 : Int totalCFs=muellerElements.shape().product()*freqValues.shape().product()*wValues.shape().product()*2,
223 0 : cfsDone=0;
224 :
225 0 : ProgressMeter pm(1.0, Double(totalCFs), "fillCF", "","","",true);
226 :
227 0 : for (uInt imx=0;imx<muellerElements.nelements();imx++) // Loop over all MuellerElements
228 0 : for (uInt imy=0;imy<muellerElements(imx).nelements();imy++)
229 : {
230 : {
231 0 : for (uInt inu=0;inu<freqValues.nelements();inu++) // All freq. channels
232 : {
233 : Float sampling, samplingWt;
234 : Int xSupport, ySupport, xSupportWt, ySupportWt;
235 0 : CoordinateSystem cs_l;
236 0 : String bandName;
237 : // Extract the parameters index by (MuellerElement, Freq, W)
238 0 : cfWtb.getParams(cs_l, samplingWt, xSupportWt, ySupportWt, bandName,
239 0 : freqValues(inu),
240 : // wValues(iw),
241 0 : wValues(0),
242 0 : muellerElements(imx)(imy));
243 0 : cfb.getParams(cs_l, sampling, xSupport, ySupport, bandName,
244 0 : freqValues(inu),
245 0 : wValues(0),
246 0 : muellerElements(imx)(imy));
247 0 : aTerm.setBandName(bandName);
248 :
249 : // {
250 : // Double lambdaByD = 1.22*C::c/freqValues(inu)/25.0;
251 : // Double FoV_x = fabs(skyNX*skyIncr(0));
252 : // Double FoV_y = fabs(skyNY*skyIncr(1));
253 : // Vector<Double> uvScale_l(3);
254 : // uvScale_l(0) = (FoV_x < lambdaByD) ? FoV_x : lambdaByD;
255 : // uvScale_l(1) = (FoV_y < lambdaByD) ? FoV_y : lambdaByD;
256 : // uvScale_l(2) = 0.0;
257 :
258 : // //Hints that only uvScale needs to be updated in PSTerm.
259 : // IPosition dummy;
260 : // Vector<Double> dummyoffset;
261 : // Double pss = -1;
262 : // // << "############ " << freqValues(inu) << " " << skyIncr << skyNX << " " << uvScale_l << endl;
263 : // psTerm.reinit(dummy, uvScale_l, dummyoffset,pss);
264 : // }
265 :
266 :
267 0 : IPosition pbshp(4,nx, ny,1,1);
268 : // Set the shape to 2x2 pixel images for dry gridding
269 0 : if (isDryRun) pbshp[0]=pbshp[1]=2;
270 :
271 : //
272 : // Cache the A-Term for this polarization and frequency
273 : //
274 0 : Double conjFreq=SynthesisUtils::conjFreq(freqValues(inu),imRefFreq_p);
275 : Int conjFreqIndex;
276 0 : conjFreq=SynthesisUtils::nearestValue(freqValues, conjFreq, conjFreqIndex);
277 :
278 : // cout<<"Muller Array = "<<muellerElements(imx)(imy)<<"\n" ;
279 : // USEFUL DEBUG MESSAGE
280 : // cerr << "Freq. values: "
281 : // << freqValues(inu) << " "
282 : // << imRefFreq_p << " "
283 : // << conjFreq << " "
284 : // << endl;
285 :
286 0 : CoordinateSystem conjPolCS_l=cs_l; AWConvFunc::makeConjPolAxis(conjPolCS_l);
287 0 : TempImage<Complex> ftATerm_l(pbshp, cs_l), ftATermSq_l(pbshp,conjPolCS_l);
288 : Int index;
289 0 : Vector<Int> conjPol;
290 0 : index = conjPolCS_l.findCoordinate(Coordinate::STOKES);
291 0 : conjPol = conjPolCS_l.stokesCoordinate(index).stokes();
292 : //cerr << "ConjPol = " << conjPol << endl;
293 :
294 : // {
295 : // // Vector<Double> chanFreq = vb.frequency();
296 : // CoordinateSystem skyCS(ftATerm_l.coordinates());
297 : // Int index = skyCS.findCoordinate(Coordinate::SPECTRAL);
298 : // SpectralCoordinate SpC = skyCS.spectralCoordinate(index);
299 : // Vector<Double> refVal = SpC.referenceValue();
300 :
301 : // Double ff = refVal[0];
302 : // cerr << "Freq, ConjFreq: " << freqValues(inu) << " " << conjFreq << " " << ff << endl;
303 : // }
304 :
305 :
306 0 : Bool doSquint=true;
307 : // Bool doSquint=false; Complex tt;
308 0 : ftATerm_l.set(Complex(1.0,0.0)); ftATermSq_l.set(Complex(1.0,0.0));
309 :
310 0 : Int me=muellerElements(imx)(imy);
311 0 : if (!isDryRun)
312 : {
313 0 : aTerm.applySky(ftATerm_l, vb, doSquint, 0, me, freqValues(inu));//freqHi);
314 : // {
315 : // ostringstream name;
316 : // name << "ftATerm" << "_" << inu << "_" << muellerElements(imx)(imy) <<".im";
317 : // storeImg(name,ftATerm_l);
318 : // }
319 : //tt=max(ftATerm_l.get()); ftATerm_l.put(ftATerm_l.get()/tt);
320 0 : if (conjPB_p) aTerm.applySky(ftATermSq_l, vb, doSquint, 0,me,conjFreq);
321 0 : else aTerm.applySky(ftATermSq_l, vb, doSquint, 0,me,freqValues(inu));
322 :
323 : }
324 :
325 : //tt=max(ftATermSq_l.get()); ftATermSq_l.put(abs(ftATermSq_l.get()/tt));
326 :
327 : //{
328 : // ostringstream name;
329 : // name << "ftTermSq" << "_" << muellerElements(imx)(imy) <<".im";
330 : // storeImg(name,ftATermSq_l);
331 : //}
332 : // TempImage<Complex> ftATermSq_l(pbshp,cs_l);
333 : // ftATermSq_l.set(Complex(1.0,0.0));
334 : // aTerm.applySky(ftATermSq_l, vb, false, 0);
335 : // tt=max(ftATermSq_l.get());
336 : // ftATermSq_l.put(ftATermSq_l.get()/tt);
337 :
338 0 : Int directionIndex=cs_l.findCoordinate(Coordinate::DIRECTION);
339 0 : DirectionCoordinate dc=cs_l.directionCoordinate(directionIndex);
340 0 : Vector<Double> cellSize;
341 0 : cellSize = dc.increment();
342 :
343 : //
344 : // Now compute the PS x W-Term and apply the cached
345 : // A-Term to build the full CF.
346 : //
347 0 : for (uInt iw=0;iw<wValues.nelements();iw++) // All w-planes
348 : {
349 0 : if (!isDryRun)
350 : {
351 0 : LogIO log_l(LogOrigin("AWConvFunc2", "fillConvFuncBuffer[R&D]"));
352 :
353 : log_l << " CF("
354 0 : << "M:"<<muellerElements(imx)(imy)
355 : << ",C:" << inu
356 0 : << ",W:" << iw << "): ";
357 0 : }
358 : // {
359 : // CountedPtr<CFCell> thisCell=cfb.getCFCellPtr(freqValues(inu), wValues(iw), muellerElements(imx)(imy));
360 : // thisCell->conjFreq_p = conjFreq;
361 : // cerr << "ConjFreq: " << thisCell->conjFreq_p << " " << inu << " " << iw << " " << muellerElements(imx)(imy) << endl;
362 : // }
363 :
364 0 : Array<Complex> &cfWtBuf=(*(cfWtb.getCFCellPtr(freqValues(inu), wValues(iw),
365 0 : muellerElements(imx)(imy))->storage_p));
366 0 : Array<Complex> &cfBuf=(*(cfb.getCFCellPtr(freqValues(inu), wValues(iw),
367 0 : muellerElements(imx)(imy))->storage_p));
368 : // IPosition cfWtBufShape= cfWtb.getCFCellPtr(freqValues(inu), wValues(iw),
369 : // muellerElements(imx)(imy))->shape_p;
370 : // IPosition cfBufShape=cfb.getCFCellPtr(freqValues(inu), wValues(iw),
371 : // muellerElements(imx)(imy))->shape_p;
372 :
373 0 : cfWtBuf.resize(pbshp);
374 0 : cfBuf.resize(pbshp);
375 0 : const Vector<Double> sampling_l(2,sampling);
376 : // Double wval = wValues[iw];
377 0 : Matrix<Complex> cfBufMat(cfBuf.nonDegenerate()),
378 0 : cfWtBufMat(cfWtBuf.nonDegenerate());
379 : //
380 : // Apply the Prolate Spheroidal and W-Term kernels
381 : //
382 :
383 0 : Vector<Double> s(2); s=sampling;
384 : // Int inner = cfBufMat.shape()(0)/aTerm.getOversampling();
385 : // Float inner = 2.0*aTerm.getOversampling()/cfBufMat.shape()(0);
386 :
387 : //Timer tim;
388 : //tim.mark();
389 0 : if (psTerm.isNoOp() || isDryRun)
390 0 : cfBufMat = cfWtBufMat = 1.0;
391 : else
392 : {
393 : //psTerm.applySky(cfBufMat, false); // Assign (psScale set in psTerm.init()
394 : //psTerm.applySky(cfWtBufMat, false); // Assign
395 0 : psTerm.applySky(cfBufMat, s, cfBufMat.shape()(0)/s(0)); // Assign (psScale set in psTerm.init()
396 0 : psTerm.applySky(cfWtBufMat, s, cfWtBufMat.shape()(0)/s(0)); // Assign
397 :
398 0 : cfWtBuf *= cfWtBuf;
399 : }
400 : //tim.show("PSTerm*2: ");
401 :
402 : // WBAWP CODE BEGIN -- make PS*PS for Weights
403 : // psTerm.applySky(cfWtBufMat, true); // Multiply
404 : // WBAWP CODE END
405 :
406 : // psTerm.applySky(cfBufMat, s, inner/2.0);//pbshp(0)/(os));
407 : // psTerm.applySky(cfWtBufMat, s, inner/2.0);//pbshp(0)/(os));
408 :
409 : // W-term is a unit-amplitude term in the image
410 : // doimain. No need to apply it to the
411 : // wt-functions.
412 :
413 : //tim.mark();
414 0 : if (!isDryRun)
415 : {
416 0 : wTerm.applySky(cfBufMat, iw, cellSize, wScale, cfBuf.shape()(0));///4);
417 : //cerr << iw << " " << cellSize << " " << iw*iw/wScale << endl;
418 : }
419 : //tim.show("WTerm: ");
420 : // wTerm.applySky(cfWtBufMat, iw, cellSize, wScale, cfWtBuf.shape()(0)/4);
421 :
422 0 : IPosition PolnPlane(4,0,0,0,0),
423 0 : pbShape(4, cfBuf.shape()(0), cfBuf.shape()(1), 1, 1);
424 : //
425 : // Make TempImages and copy the buffers with PS *
426 : // WKernel applied (too bad that TempImages can't be
427 : // made with existing buffers)
428 : //
429 : //-------------------------------------------------------------
430 0 : TempImage<Complex> twoDPB_l(pbShape, cs_l);
431 0 : TempImage<Complex> twoDPBSq_l(pbShape,cs_l);
432 : //-------------------------------------------------------------
433 : // WBAWP CODE BEGIN -- ftATermSq_l has conj. PolCS
434 0 : cfWtBuf *= ftATerm_l.get()*conj(ftATermSq_l.get());
435 : //tim.mark();
436 : //////TESTOO/////////////
437 : /* {
438 : String tmpname=File::newUniqueName("./", "WTerm").baseName();
439 : cerr << "WTERM image " << tmpname << endl;
440 : PagedImage<Complex> tempB(pbShape, cs_l, tmpname);
441 : tempB.putSlice(cfBufMat, PolnPlane);
442 :
443 :
444 : }*/
445 : //////////////////////
446 :
447 :
448 :
449 : //UUU cfWtBuf *= ftATerm_l.get();
450 0 : cfBuf *= ftATerm_l.get();
451 : //tim.show("W*A*2: ");
452 : // WBAWP CODE END
453 :
454 :
455 :
456 : // cfWtBuf = sqrt(cfWtBuf);
457 : // psTerm.applySky(cfWtBufMat,true);
458 :
459 : //tim.mark();
460 0 : twoDPB_l.putSlice(cfBuf, PolnPlane);
461 0 : twoDPBSq_l.putSlice(cfWtBuf, PolnPlane);
462 : //tim.show("putSlice:");
463 : // WBAWP CODE BEGIN
464 : // twoDPB_l *= ftATerm_l;
465 : // WBAWP CODE END
466 :
467 : // twoDPBSq_l *= ftATermSq_l;//*conj(ftATerm_l);
468 :
469 : // To accumulate avgPB2, call this function.
470 : // PBSQWeight
471 0 : Bool PBSQ = false;
472 0 : if(PBSQ) makePBSq(twoDPBSq_l);
473 :
474 :
475 : //
476 : // Set the ref. freq. of the co-ordinate system to
477 : // that set by ATerm::applySky().
478 : //
479 : //tim.mark();
480 0 : CoordinateSystem cs=twoDPB_l.coordinates();
481 0 : Int index= twoDPB_l.coordinates().findCoordinate(Coordinate::SPECTRAL);
482 0 : SpectralCoordinate SpCS = twoDPB_l.coordinates().spectralCoordinate(index);
483 :
484 0 : Double cfRefFreq=SpCS.referenceValue()(0);
485 0 : Vector<Double> refValue; refValue.resize(1); refValue(0)=cfRefFreq;
486 0 : SpCS.setReferenceValue(refValue);
487 0 : cs.replaceCoordinate(SpCS,index);
488 : //tim.show("CSStuff:");
489 : // {
490 : // ostringstream name;
491 : // name << "twoDPB.before" << iw << "_" << inu << "_" << muellerElements(imx)(imy) <<".im";
492 : // storeImg(name,twoDPB_l);
493 : // name << "twoDPBSq.before" << iw << "_" << inu << "_" << muellerElements(imx)(imy) <<".im";
494 : // storeImg(name,twoDPBSq_l);
495 : // }
496 : //
497 : // Now FT the function and copy the data from
498 : // TempImages back to the CFBuffer buffers
499 : //
500 : //tim.mark();
501 0 : if (!isDryRun)
502 : {
503 0 : LatticeFFT::cfft2d(twoDPB_l);
504 0 : LatticeFFT::cfft2d(twoDPBSq_l);
505 : }
506 : //tim.show("FFT*2:");
507 : // Array<Complex> t0;
508 : // twoDPBSq_l.get(t0); t0 = abs(t0);
509 : // twoDPBSq_l.put(t0);
510 :
511 :
512 : //tim.mark();
513 0 : IPosition shp(twoDPB_l.shape());
514 0 : IPosition start(4, 0, 0, 0, 0), pbSlice(4, shp[0]-1, shp[1]-1,1/*polInUse*/, 1),
515 0 : sliceLength(4,cfBuf.shape()[0]-1,cfBuf.shape()[1]-1,1,1);
516 0 : cfBuf(Slicer(start,sliceLength)).nonDegenerate()
517 0 : =(twoDPB_l.getSlice(start, pbSlice, true));
518 :
519 0 : shp = twoDPBSq_l.shape();
520 0 : IPosition pbSqSlice(4, shp[0]-1, shp[1]-1, 1, 1),
521 0 : sqSliceLength(4,cfWtBuf.shape()(0)-1,cfWtBuf.shape()[1]-1,1,1);
522 :
523 0 : cfWtBuf(Slicer(start,sqSliceLength)).nonDegenerate()
524 0 : =(twoDPBSq_l.getSlice(start, pbSqSlice, true));
525 : //tim.show("Slicer*2:");
526 : //
527 : // Finally, resize the buffers, limited to the
528 : // support size determined by the threshold
529 : // suppled by the ATerm (done internally in
530 : // resizeCF()). Transform the co-ord. system to
531 : // the FT domain set the co-ord. sys. and modified
532 : // support sizes.
533 : //
534 : //tim.mark();
535 0 : Int supportBuffer = (Int)(getOversampling(psTerm, wTerm, aTerm)*2.0);
536 0 : if (!isDryRun)
537 : {
538 0 : if (iw==0) wtcpeak = max(cfWtBuf);
539 0 : cfWtBuf /= wtcpeak;
540 : }
541 : //tim.show("Norm");
542 :
543 : //tim.mark();
544 0 : if (!isDryRun)
545 0 : AWConvFunc::resizeCF(cfWtBuf, xSupportWt, ySupportWt, supportBuffer, samplingWt,0.0);
546 : //log_l << "CF WT Support: " << xSupport << " (" << xSupportWt << ") " << "pixels" << LogIO::POST;
547 : //tim.show("Resize:");
548 :
549 : //tim.mark();
550 0 : Vector<Double> ftRef(2);
551 : // ftRef(0)=cfWtBuf.shape()(0)/2-1;
552 : // ftRef(1)=cfWtBuf.shape()(1)/2-1;
553 0 : ftRef(0)=cfWtBuf.shape()(0)/2.0;
554 0 : ftRef(1)=cfWtBuf.shape()(1)/2.0;
555 0 : CoordinateSystem ftCoords=cs_l;
556 0 : if (isDryRun)
557 : {
558 0 : ftRef(0)=nx/2.0;
559 0 : ftRef(1)=ny/2.0;
560 0 : SynthesisUtils::makeFTCoordSys(cs_l, nx,ftRef , ftCoords);
561 : }
562 : else
563 0 : SynthesisUtils::makeFTCoordSys(cs_l, cfWtBuf.shape()(0), ftRef, ftCoords);
564 :
565 0 : CountedPtr<CFCell> cfCellPtr;
566 0 : cfWtb.setParams(inu,iw,imx,imy,//muellerElements(imx)(imy),
567 0 : freqValues(inu), String(""), wValues(iw), muellerElements(imx)(imy),
568 : ftCoords, samplingWt, xSupportWt, ySupportWt,
569 0 : String(""), // Default ==> don't set it in the CFCell
570 0 : conjFreq, conjPol[0]);
571 :
572 0 : cfCellPtr = cfWtb.getCFCellPtr(freqValues(inu), wValues(iw),
573 0 : muellerElements(imx)(imy));
574 0 : cfCellPtr->pa_p=Quantity(vbPA,"rad");
575 0 : cfCellPtr->telescopeName_p = aTerm.getTelescopeName();
576 0 : cfCellPtr->isRotationallySymmetric_p = aTerm.isNoOp();
577 : //cerr << "AWConvFunc: Telescope name = " << cfCellPtr->telescopeName_p << " " << aTerm.getTelescopeName() << endl;
578 : //tim.show("CSStuff:");
579 : // setUpCFSupport(cfBuf, xSupport, ySupport, sampling);
580 : // if (iw==0)
581 : //tim.mark();
582 : //Int supportBuffer = (Int)(aTerm->getOversampling()*1.5);
583 :
584 0 : if (!isDryRun)
585 : {
586 0 : cpeak = max(cfBuf);
587 0 : cfBuf /= cpeak;
588 : }
589 : //tim.show("Peaknorm:");
590 : // {
591 : // ostringstream name;
592 : // name << "twoDPB.after" << iw << "_" << inu << "_" << muellerElements(imx)(imy) << ".im";
593 : // storeImg(name,twoDPB_l);
594 : // // name << "twoDPBSq.after" << iw << "_" << inu << "_" << muellerElements(imx)(imy) << ".im";
595 : // // storeImg(name,twoDPBSq_l);
596 : // }
597 :
598 0 : if (!isDryRun)
599 0 : AWConvFunc::resizeCF(cfBuf, xSupport, ySupport, supportBuffer, sampling,0.0);
600 :
601 :
602 :
603 0 : if (!isDryRun)
604 : {
605 0 : LogIO log_l(LogOrigin("AWConvFunc2", "fillConvFuncBuffer[R&D]"));
606 0 : log_l << "CF Support: " << xSupport << " (" << xSupportWt << ") " << "pixels" << LogIO::POST;
607 0 : }
608 :
609 : // cfb.getCFCellPtr(freqValues(inu), wValues(iw), muellerElement)->storage_p->assign(cfBuf);
610 : // ftRef(0)=cfBuf.shape()(0)/2-1;
611 : // ftRef(1)=cfBuf.shape()(1)/2-1;
612 0 : ftRef(0)=cfBuf.shape()(0)/2.0;
613 0 : ftRef(1)=cfBuf.shape()(1)/2.0;
614 :
615 : //tim.mark();
616 0 : cfNorm=cfWtNorm=1.0;
617 : //if ((iw == 0) && (!isDryRun))
618 0 : if (!isDryRun)
619 : {
620 0 : cfNorm=0; cfWtNorm=0;
621 0 : cfNorm = AWConvFunc::cfArea(cfBufMat, xSupport, ySupport, sampling);
622 0 : cfWtNorm = AWConvFunc::cfArea(cfWtBufMat, xSupportWt, ySupportWt, sampling);
623 : }
624 : //tim.show("Area*2:");
625 :
626 : //tim.mark();
627 0 : if (cfNorm != Complex(0.0)) cfBuf /= cfNorm;
628 0 : if (cfWtNorm != Complex(0.0)) cfWtBuf /= cfWtNorm;
629 :
630 : //tim.show("cfNorm*2:");
631 :
632 : //tim.mark();
633 0 : ftCoords=cs_l;
634 0 : if (isDryRun)
635 : {
636 0 : ftRef(0) = nx/2.0;
637 0 : ftRef(1) = ny/2.0;
638 0 : SynthesisUtils::makeFTCoordSys(cs_l, nx, ftRef, ftCoords);
639 : }
640 : else
641 0 : SynthesisUtils::makeFTCoordSys(cs_l, cfBuf.shape()(0), ftRef, ftCoords);
642 :
643 0 : cfb.setParams(inu,iw,imx,imy,//muellerElements(imx)(imy),
644 0 : freqValues(inu), String(""), wValues(iw), muellerElements(imx)(imy),
645 : ftCoords, sampling, xSupport, ySupport,
646 0 : String(""), // Default ==> Don't set in the CFCell
647 0 : conjFreq, conjPol[0]);
648 0 : cfCellPtr=cfb.getCFCellPtr(freqValues(inu), wValues(iw),
649 0 : muellerElements(imx)(imy));
650 0 : cfCellPtr->pa_p=Quantity(vbPA,"rad");
651 0 : cfCellPtr->telescopeName_p = aTerm.getTelescopeName();
652 0 : cfCellPtr->isRotationallySymmetric_p = aTerm.isNoOp();
653 : //
654 : // Now tha the CFs have been computed, cache its
655 : // paramters in CFCell for quick access in tight
656 : // loops (in the re-sampler, e.g.).
657 : //
658 :
659 0 : (cfWtb.getCFCellPtr(freqValues(inu), wValues(iw),
660 0 : muellerElements(imx)(imy)))->initCache(isDryRun);
661 0 : (cfb.getCFCellPtr(freqValues(inu), wValues(iw),
662 0 : muellerElements(imx)(imy)))->initCache(isDryRun);
663 :
664 0 : pm.update((Double)cfsDone++);
665 : //tim.show("End*2:");
666 0 : }
667 0 : }
668 : }
669 : }
670 0 : }
671 : //
672 : //----------------------------------------------------------------------
673 : //
674 0 : Complex AWConvFunc::cfArea(Matrix<Complex>& cf,
675 : const Int& xSupport, const Int& ySupport,
676 : const Float& sampling)
677 : {
678 0 : Complex cfNorm=0;
679 0 : Int origin=cf.shape()(0)/2;
680 0 : Float peak=0;
681 0 : IPosition ndx(4,0,0,0,0);
682 0 : IPosition peakPix(ndx);
683 0 : peakPix = 0;
684 0 : for(ndx(1)=0;ndx(1)<cf.shape()(1);ndx(1)++)
685 0 : for(ndx(0)=0;ndx(0)<cf.shape()(0);ndx(0)++)
686 0 : if (abs(cf(ndx)) > peak) {peakPix = ndx;peak=abs(cf(ndx));}
687 : // origin = peakPix(0);
688 0 : if (origin != peakPix(0))
689 : {
690 0 : LogIO log_l(LogOrigin("AWConvFunc2","cfArea"));
691 0 : log_l << "Peak not at the center " << origin << " " << cf(IPosition(4,origin,origin,0,0)) << " " << peakPix << " " << peak << LogIO::NORMAL1;
692 : // peakNIC=1e7;
693 0 : }
694 0 : for (Int ix=-xSupport;ix<xSupport;ix++)
695 0 : for (int iy=-ySupport;iy<ySupport;iy++)
696 : {
697 : //cfNorm += Complex(real(cf(ix*(Int)sampling+origin, iy*(Int)sampling+origin)),0.0);
698 0 : cfNorm += (cf(ix*(Int)sampling+origin, iy*(Int)sampling+origin));
699 : // cerr << cfNorm << " " << ix << " " << iy << " " << ix*(Int)sampling+origin << " " << iy*(Int)sampling+origin
700 : // << real(cf(ix*(Int)sampling+origin, iy*(Int)sampling+origin)) << endl;
701 : }
702 : // cf /= cfNorm;
703 0 : return cfNorm;
704 0 : }
705 : //
706 : //----------------------------------------------------------------------
707 : //
708 0 : Vector<Double> AWConvFunc::makeWValList(const Double &dW, const Int &nW)
709 : {
710 0 : Vector<Double> wValues(nW);
711 : // for (Int iw=0;iw<nW;iw++) wValues[iw]=iw*dW;
712 0 : wValues = 0.0;
713 0 : if (dW > 0.0)
714 0 : for (Int iw=0;iw<nW;iw++) wValues[iw]=iw*iw/dW;
715 0 : return wValues;
716 0 : }
717 :
718 : // This methods is depcricated. Keeping it here since it *might*
719 : // have use sometime later and therefore want to push it on to SVN
720 : // before deleting it form the active version of this file.
721 0 : Matrix<Double> AWConvFunc::getFreqRangePerSpw(const VisBuffer2& vb)
722 : {
723 : //
724 : // Find the total effective bandwidth
725 : //
726 0 : Cube<Double> fminmax;
727 0 : Double fMax=0, fMin=MAX_FREQ;
728 0 : ArrayColumn<Double> spwCol=vb.subtableColumns().spectralWindow().chanFreq();
729 0 : fminmax.resize(spwChanSelFlag_p.shape()(0),spwChanSelFlag_p.shape()(1),2);
730 0 : fminmax=0;
731 0 : for (uInt ims=0; ims<spwChanSelFlag_p.shape()(0); ims++)
732 0 : for(uInt ispw=0; ispw<spwChanSelFlag_p.shape()(1); ispw++)
733 : {
734 0 : fMax=0, fMin=MAX_FREQ;
735 0 : for(uInt ichan=0; ichan<spwChanSelFlag_p.shape()(2); ichan++)
736 : {
737 0 : if (spwChanSelFlag_p(ims,ispw,ichan)==1)
738 : {
739 0 : Slicer slicer(IPosition(1,ichan), IPosition(1,1));
740 0 : Vector<Double> freq = spwCol(ispw)(slicer);
741 0 : if (freq(0) < fMin) fMin = freq(0);
742 0 : if (freq(0) > fMax) fMax = freq(0);
743 0 : }
744 : }
745 0 : fminmax(ims,ispw,0)=fMin;
746 0 : fminmax(ims,ispw,1)=fMax;
747 : }
748 :
749 0 : Matrix<Double> freqRangePerSpw(fminmax.shape()(1),2);
750 0 : for (uInt j=0;j<fminmax.shape()(1);j++) // SPW
751 : {
752 0 : freqRangePerSpw(j,0)=0;
753 0 : freqRangePerSpw(j,1)=MAX_FREQ;
754 0 : for (uInt i=0;i<fminmax.shape()(0);i++) //MSes
755 : {
756 0 : if (freqRangePerSpw(j,0) < fminmax(i,j,0)) freqRangePerSpw(j,0)=fminmax(i,j,0);
757 0 : if (freqRangePerSpw(j,1) > fminmax(i,j,1)) freqRangePerSpw(j,1)=fminmax(i,j,1);
758 : }
759 : }
760 0 : for(uInt i=0;i<freqRangePerSpw.shape()(0);i++)
761 : {
762 0 : if (freqRangePerSpw(i,0) == MAX_FREQ) freqRangePerSpw(i,0)=-1;
763 0 : if (freqRangePerSpw(i,1) == 0) freqRangePerSpw(i,1)=-1;
764 : }
765 :
766 0 : return freqRangePerSpw;
767 0 : }
768 : //
769 : //----------------------------------------------------------------------
770 : // Given the VB and the uv-grid, make a list of frequency values to
771 : // sample the frequency axis of the CFBuffer. Typically, this will
772 : // be determined by the bandwidth-smearning limit.
773 : //
774 : // This limit is (deltaNu/Nu) * sqrt(l^2 + m^2) < ResolutionElement.
775 : // Translating max. distance from the phase center to field-of-view
776 : // of the supplied image, and converting Resolution Element to
777 : // 1/Cellsize, this expression translates to deltaNU<FMin/Nx (!)
778 0 : Vector<Double> AWConvFunc::makeFreqValList(Double &dNU,
779 : const VisBuffer2& vb,
780 : const ImageInterface<Complex>& uvGrid,
781 : Vector<String>& bandNames)
782 : {
783 : (void)uvGrid; (void)dNU; (void)vb;
784 0 : Vector<Double> fValues;
785 0 : Int nSpw = spwFreqSelection_p.shape()(0);
786 0 : if (wbAWP_p==false)
787 : {
788 : // Return the sky-image ref. freq.
789 0 : fValues.resize(1);
790 0 : fValues[0]=imRefFreq_p;
791 : }
792 : else
793 : {
794 0 : fValues.resize(nSpw);
795 0 : for(Int i=0;i<nSpw;i++)
796 0 : fValues(i)=spwFreqSelection_p(i,2);
797 : }
798 :
799 0 : bandNames.resize(nSpw);
800 0 : ScalarColumn<String> spwNames=vb.subtableColumns().spectralWindow().name();
801 0 : for(Int i=0;i<nSpw;i++)
802 : {
803 0 : int s=spwFreqSelection_p(i,0);
804 : // LogIO os;
805 : // os << "Spw"<<s<<": " << spwNames.getColumn()[s]
806 : // << " " << s << " " << nSpw
807 : // << LogIO::WARN;
808 :
809 0 : bandNames(i)=spwNames.getColumn()[s];
810 : }
811 0 : return fValues;
812 0 : }
813 : //
814 : //----------------------------------------------------------------------
815 : //
816 0 : void AWConvFunc::makeConvFunction(const ImageInterface<Complex>& image,
817 : const VisBuffer2& vb,
818 : const Int wConvSize,
819 : const CountedPtr<PolOuterProduct>& pop,
820 : const Float pa,
821 : const Float dpa,
822 : const Vector<Double>& uvScale,
823 : const Vector<Double>& ,//uvOffset,
824 : const Matrix<Double>& ,//vbFreqSelection,
825 : CFStore2& cfs2,
826 : CFStore2& cfwts2,
827 : Bool fillCF)
828 : {
829 0 : LogIO log_l(LogOrigin("AWConvFunc2", "makeConvFunction[R&D]"));
830 : Int convSize, convSampling, polInUse;
831 0 : Double wScale=0.0;
832 0 : Array<Complex> convFunc_l, convWeights_l;
833 0 : Double cfRefFreq=-1, freqScale=1e8;
834 0 : Quantity paQuant(pa,"rad");
835 :
836 :
837 0 : Int nx=image.shape()(0);//, ny=image.shape()(1);
838 0 : Vector<Double> skyIncr;
839 :
840 : log_l << "Making a new convolution function for PA="
841 0 : << pa*(180/C::pi) << "deg"
842 0 : << " for field ID " << vb.fieldId()(0);
843 : // log_l << "TimeStamps(0-10) ";
844 : // for (Int i=0;i<10;i++)
845 : // //log_l << MVTime(vb.time()(i)).string(MVTime::TIME) << " ";
846 : // log_l << vb.time()(i)/1e8 << " ";
847 0 : log_l << LogIO::NORMAL << LogIO::POST;
848 :
849 0 : if(wConvSize>0)
850 : {
851 0 : log_l << "Using " << wConvSize << " planes for W-projection" << LogIO::POST;
852 : Double maxUVW;
853 0 : float WFUDGE=4.0;
854 0 : WFUDGE=refim::SynthesisUtils::getenv("WTerm.WFUDGE",WFUDGE);
855 :
856 : //maxUVW=0.25/abs(image.coordinates().increment()(0));
857 0 : maxUVW=1.0/abs(image.coordinates().increment()(0)*WFUDGE);
858 : log_l << "Estimating maximum possible W = " << maxUVW
859 0 : << " (wavelengths)" << LogIO::POST;
860 :
861 : // Double invLambdaC=vb.getFrequencies(0)(0)/C::c;
862 : // Int nFreq = (vb.getFrequencies(0).nelements())-1;
863 : // Double invMinL = vb.getFrequencies(0)(nFreq)/C::c;
864 : // log_l << "wavelength range = " << 1.0/invLambdaC << " (m) to "
865 : // << 1.0/invMinL << " (m)" << LogIO::POST;
866 0 : if (wConvSize > 1)
867 : {
868 0 : wScale=Float((wConvSize-1)*(wConvSize-1))/maxUVW;
869 : log_l << "Scaling in W (at maximum W) = " << 1.0/wScale
870 0 : << " wavelengths per pixel" << LogIO::POST;
871 : }
872 : }
873 : //
874 : // Get the coordinate system
875 : //
876 0 : CoordinateSystem coords(image.coordinates());
877 : //
878 : // Set up the convolution function.
879 : //
880 0 : convSampling=getOversampling(*psTerm_p, *wTerm_p, *aTerm_p);
881 0 : convSize=aTerm_p->getConvSize();
882 : // cout<<"Conv Sampling listed in aipsrc is : "<<convSampling<<endl;
883 : // cout<<"Conv Size is : "<<convSize<<endl;
884 : //
885 : // Make a two dimensional image to calculate auto-correlation of
886 : // the ideal illumination pattern. We want this on a fine grid in
887 : // the UV plane
888 : //
889 0 : Int index= coords.findCoordinate(Coordinate::SPECTRAL);
890 0 : SpectralCoordinate spCS = coords.spectralCoordinate(index);
891 0 : imRefFreq_p=spCS.referenceValue()(0);
892 :
893 0 : index=coords.findCoordinate(Coordinate::DIRECTION);
894 0 : AlwaysAssert(index>=0, AipsError);
895 0 : DirectionCoordinate dc=coords.directionCoordinate(index);
896 0 : Vector<Double> sampling;
897 0 : skyIncr = sampling = dc.increment();
898 : // cout<<"The image sampling is set to :"<<sampling<<endl;
899 0 : sampling*=Double(convSampling);
900 0 : sampling*=Double(nx)/Double(convSize);
901 : // cout<<"The resampled increment is :"<<sampling<<endl;
902 0 : dc.setIncrement(sampling);
903 :
904 0 : Vector<Double> unitVec(2);
905 0 : unitVec=convSize/2;
906 0 : dc.setReferencePixel(unitVec);
907 :
908 : // Set the reference value to that of the image
909 0 : coords.replaceCoordinate(dc, index);
910 : //
911 : // Make an image with circular polarization axis. Return the
912 : // no. of vis. poln. planes that will be used in making the user
913 : // defined Stokes image.
914 : //
915 0 : polInUse=aTerm_p->makePBPolnCoords(vb, convSize, convSampling,
916 : image.coordinates(),nx,nx,
917 : coords);//,feedStokes_l);
918 : //------------------------------------------------------------------
919 : // Make the sky Stokes PB. This will be used in the gridding
920 : // correction
921 : //------------------------------------------------------------------
922 0 : IPosition pbShape(4, convSize, convSize, polInUse, 1);
923 0 : TempImage<Complex> twoDPB(pbShape, coords);
924 0 : IPosition pbSqShp(pbShape);
925 :
926 0 : unitVec=pbSqShp[0]/2;
927 0 : dc.setReferencePixel(unitVec);
928 0 : coords.replaceCoordinate(dc, index);
929 :
930 0 : TempImage<Complex> twoDPBSq(pbSqShp,coords);
931 0 : twoDPB.set(Complex(1.0,0.0));
932 0 : twoDPBSq.set(Complex(1.0,0.0));
933 : //
934 : // Accumulate the various terms that constitute the gridding
935 : // convolution function.
936 : //
937 : //------------------------------------------------------------------
938 : // Int inner=convSize/convSampling;
939 : // CFStore2 cfs2_p, cfwts2_p;
940 0 : CountedPtr<CFBuffer> cfb_p, cfwtb_p;
941 : // cfs2.rememberATerm(aTerm_p);
942 : // cfwts2.rememberATerm(aTerm_p);
943 :
944 0 : Vector<Quantity> paList(1); paList[0]=paQuant;
945 : //
946 : // Determine the "Mueller Matrix" (called PolOuterProduct here for
947 : // a better name) elements to use based on the sky-Stokes planes
948 : // requested. PolOuterProduct::makePolMap() makes a
949 : // Matrix<Int>. The elements of this matrix has the index of the
950 : // convolution function for the pol. product. Unused elements are
951 : // set to -1. The physical definition of the PolOuterProduct
952 : // elements are as defined in Eq. 4 in A&A 487, 419-429 (2008)
953 : // (http://arxiv.org/abs/0805.0834).
954 : //
955 : // First detemine the list of Stokes requested. Then convert the
956 : // requested Stokes to the appropriate Pol cross-product. When
957 : // the off-diagonal elements of the outer-product are significant,
958 : // this will lead to more than one outer-product element per
959 : // Stokes.
960 : //
961 : // The code below still assume a diagonally dominant
962 : // outer-product. This is probably OK for antenna arrays. After the
963 : // debugging phase is over, the
964 : // Vector<PolOuterProduct::CrossCircular> should become
965 : // Matrix<PolOuterProduct> and PolOuterProduct should be
966 : // "templated" to be of type Circular or Linear.
967 : //
968 0 : StokesCoordinate skyStokesCo=coords.stokesCoordinate(coords.findCoordinate(Coordinate::STOKES));
969 0 : Vector<Int> skyStokes=skyStokesCo.stokes();
970 : //Vector<PolOuterProduct::CrossPolCircular> pp(skyStokes.nelements());
971 0 : PolMapType polMap, polIndexMap, conjPolMap, conjPolIndexMap;
972 0 : polMap = pop->getPolMat();
973 0 : polIndexMap = pop->getPol2CFMat();
974 0 : conjPolMap = pop->getConjPolMat();
975 0 : conjPolIndexMap = pop->getConjPol2CFMat();
976 :
977 : //cerr << "AWCF: " << polMap << endl << polIndexMap << endl << conjPolMap << endl << conjPolIndexMap << endl;
978 :
979 : // for(uInt ip=0;ip<pp.nelements();ip++)
980 : // pp(ip)=translateStokesToCrossPol(skyStokes(ip));
981 :
982 : // PolOuterProduct pOP; pOP.makePolMap(pp);
983 : // const Matrix<Int> muellerMatrix=pOP.getPolMap();
984 :
985 0 : Vector<Double> wValues = makeWValList(wScale, wConvSize);
986 0 : Vector<String> bandNames;
987 0 : Vector<Double> freqValues = makeFreqValList(freqScale,vb,image,bandNames);
988 0 : log_l << "Making " << wValues.nelements() << " w plane(s). " << LogIO::POST;
989 0 : log_l << "Making " << freqValues.nelements() << " frequency plane(s)." << LogIO::POST;
990 : //
991 : // If w-term is unity, we can scale the A-term with frequency. So
992 : // compute it only for the highest frequency involved.
993 : //
994 : //log_l << "Disabled scaling of CFs" << LogIO::WARN << LogIO::POST;
995 : // if (wConvSize <= 1)
996 : // {
997 : // Double rFreq = max(freqValues);
998 : // if (freqValues.nelements() > 1)
999 : // freqScale=2*(rFreq-min(freqValues));
1000 : // freqValues.resize(1);freqValues(0)=rFreq;
1001 : // }
1002 : log_l << "CFB Freq. axis [N, Min, Max, Incr. (GHz)]: "
1003 : << freqValues.nelements() << " "
1004 0 : << min(freqValues)/1e9 << " "
1005 0 : << max(freqValues)/1e9 << " "
1006 : << freqScale/1e9
1007 0 : << LogIO::POST;
1008 : //
1009 : // Re-size the CFStore object. It holds CFBuffers index by PA and
1010 : // list of unique baselines (all possible pairs of unique antenna
1011 : // types).
1012 : //
1013 0 : Matrix<Int> uniqueBaselineTypeList=makeBaselineList(aTerm_p->getAntTypeList());
1014 : //Quantity dPA(360.0,"deg");
1015 0 : Quantity dPA(dpa,"rad");
1016 0 : Int totalCFs=uniqueBaselineTypeList.shape().product()*wConvSize*freqValues.nelements()*polMap.shape().product()*2;
1017 0 : ProgressMeter pm(1.0, Double(totalCFs), "makeCF", "","","",true);
1018 0 : int cfDone=0;
1019 0 : for(Int ib=0;ib<uniqueBaselineTypeList.shape()(0);ib++)
1020 : {
1021 0 : Vector<Int> pos;
1022 0 : pos=cfs2.resize(paQuant, dPA, uniqueBaselineTypeList(ib,0), uniqueBaselineTypeList(ib,1));
1023 0 : pos=cfwts2.resize(paQuant, dPA, uniqueBaselineTypeList(ib,0), uniqueBaselineTypeList(ib,1));
1024 : //
1025 : // Re-size the CFBuffer object. It holds the 2D convolution
1026 : // functions index by (FreqValue, WValue, MuellerElement).
1027 : //
1028 0 : cfb_p=cfs2.getCFBuffer(paQuant, dPA, uniqueBaselineTypeList(ib,0),uniqueBaselineTypeList(ib,1));
1029 0 : cfwtb_p=cfwts2.getCFBuffer(paQuant, dPA, uniqueBaselineTypeList(ib,0),uniqueBaselineTypeList(ib,1));
1030 0 : cfb_p->setPointingOffset(pixFieldGrad_p);
1031 : // cfb_p->resize(wValues,freqValues,muellerMatrix);
1032 : // cfwtb_p->resize(wValues,freqValues,muellerMatrix);
1033 :
1034 0 : cfb_p->resize(wScale, freqScale, wValues,freqValues,polMap, polIndexMap,conjPolMap, conjPolIndexMap);
1035 0 : cfwtb_p->resize(wScale, freqScale, wValues,freqValues,polMap, polIndexMap,conjPolMap, conjPolIndexMap);
1036 :
1037 0 : IPosition start(4, 0, 0, 0, 0);
1038 0 : IPosition pbSlice(4, convSize, convSize, 1, 1);
1039 :
1040 0 : Matrix<Complex> screen(convSize, convSize);
1041 : // WTerm wterm_l;
1042 : // PSTerm psTerm_l;
1043 :
1044 : //Initiate construction of the "ConvolveGridder" object inside
1045 : //PSTerm. This is, for historical reasons, used to access the
1046 : //"standard" Prolate Spheroidal function implementaion. This
1047 : //however should be replaced with a simpler access, direct
1048 : //access the PS function implementation (in Utils.h
1049 : //SynthesisUtils::libreSpheroidal() - but this needs more
1050 : //testing).
1051 0 : Float psScale = (2.0*coords.increment()(0))/(nx*image.coordinates().increment()(0)),
1052 0 : innerQuaterFraction=1.0;
1053 : {
1054 0 : Int inner=convSize/(convSampling);
1055 : // Float psScale= (image.coordinates().increment()(0)*nx) /
1056 : // (coords.increment()(0)*screen.shape()(0));
1057 :
1058 : // psScale when using SynthesisUtils::libreSpheroidal() is
1059 : // 2.0/nSupport. nSupport is in pixels and the 2.0 is due to
1060 : // the center being at Nx/2. Here the nSupport is determined
1061 0 : innerQuaterFraction=refim::SynthesisUtils::getenv("AWCF.FUDGE",innerQuaterFraction);
1062 :
1063 0 : Double lambdaByD = innerQuaterFraction*1.22*C::c/min(freqValues)/25.0;
1064 0 : Double FoV_x = fabs(nx*skyIncr(0));
1065 0 : Double FoV_y = fabs(nx*skyIncr(1));
1066 0 : Vector<Double> uvScale_l(3);
1067 0 : uvScale_l(0) = (FoV_x < lambdaByD) ? FoV_x : lambdaByD;
1068 0 : uvScale_l(1) = (FoV_y < lambdaByD) ? FoV_y : lambdaByD;
1069 0 : uvScale_l(2) = uvScale(2);
1070 : // by the sky-image and is equal to convSize/convSampling.
1071 0 : psScale = 2.0/(innerQuaterFraction*convSize/convSampling);// nx*image.coordinates().increment()(0)*convSampling/2;
1072 0 : Vector<Double> uvOffset_cf(3,0); uvOffset_cf(0)=uvOffset_cf(2)=convSize/2;
1073 : // psTerm_p->init(IPosition(2,inner,inner), uvScale, uvOffset_cf,psScale);
1074 0 : psTerm_p->init(IPosition(2,inner,inner), uvScale_l, uvOffset_cf,psScale);
1075 0 : }
1076 :
1077 0 : MuellerElementType muellerElement(0,0);
1078 0 : CoordinateSystem cfb_cs=coords;
1079 : //
1080 : // Set up the Mueller matrix, the co-ordinate system, freq, and
1081 : // wvalues in the CFBuffer for the currenct CFStore object.
1082 : //
1083 : //cerr<<"Mueller matrix of row length:"<<polMap.nelements()<<" at the start of the CFBuf Loop" <<endl;
1084 0 : for (Int iw=0;iw<wConvSize;iw++)
1085 : {
1086 0 : for(uInt inu=0;inu<freqValues.nelements(); inu++)
1087 : {
1088 0 : Int npol=0;
1089 0 : for (uInt ipolx=0;ipolx<polMap.nelements();ipolx++)
1090 : {
1091 0 : npol=0;
1092 0 : for (uInt ipoly=0;ipoly<polMap(ipolx).nelements();ipoly++)
1093 : {
1094 : // Now make a CS with a single appropriate
1095 : // polarization axis per Mueller element
1096 0 : Vector<Int> whichStokes(1,skyStokes(npol++));
1097 0 : Int sIndex=cfb_cs.findCoordinate(Coordinate::STOKES);
1098 0 : StokesCoordinate stokesCS=cfb_cs.stokesCoordinate(sIndex);
1099 0 : Int fIndex=coords.findCoordinate(Coordinate::SPECTRAL);
1100 0 : SpectralCoordinate spCS = coords.spectralCoordinate(fIndex);
1101 0 : Vector<Double> refValue, incr;
1102 0 : refValue = spCS.referenceValue();
1103 0 : incr = spCS.increment();
1104 0 : cfRefFreq=freqValues(inu);
1105 0 : refValue=cfRefFreq;
1106 0 : spCS.setReferenceValue(refValue);
1107 :
1108 0 : stokesCS.setStokes(whichStokes);
1109 0 : cfb_cs.replaceCoordinate(stokesCS,sIndex);
1110 0 : cfb_cs.replaceCoordinate(spCS,fIndex);
1111 : //
1112 : // Set the various axis-parameters for the CFBuffer.
1113 : //
1114 0 : Float s=convSampling;
1115 : // cfb_p->setParams(convSize,convSize,cfb_cs,s,
1116 : // convSize, convSize,
1117 : // freqValues(inu), wValues(iw), polMap(ipolx)(ipoly));
1118 : // cfwtb_p->setParams(convSize,convSize,cfb_cs,s,
1119 : // convSize, convSize,
1120 : // freqValues(inu), wValues(iw), polMap(ipolx)(ipoly));
1121 0 : cfb_p->setParams(inu, iw, ipolx,ipoly,//polMap(ipolx)(ipoly),
1122 0 : freqValues(inu), bandNames(inu), wValues(iw), polMap(ipolx)(ipoly),
1123 : cfb_cs, s, convSize, convSize);
1124 0 : cfwtb_p->setParams(inu, iw, ipolx,ipoly,//polMap(ipolx)(ipoly),
1125 0 : freqValues(inu), bandNames(inu), wValues(iw), polMap(ipolx)(ipoly),
1126 : cfb_cs, s, convSize, convSize);
1127 0 : pm.update((Double)cfDone++);
1128 0 : }
1129 : }
1130 : } // End of loop over Mueller elements.
1131 : } // End of loop over w
1132 : //
1133 : // By this point, the all the 4 axis (Time/PA, Freq, Pol,
1134 : // Baseline) of the CFBuffer objects have been setup. The CFs
1135 : // will now be filled using the supplied PS-, W- ad A-term objects.
1136 : //
1137 0 : if (fillCF) log_l << "Making CFs for baseline type " << ib << LogIO::POST;
1138 0 : else log_l << "Making empty CFs for baseline type " << ib << LogIO::POST;
1139 : {
1140 0 : Double vbPA = getPA(vb), freqHi;
1141 :
1142 :
1143 0 : Vector<Double> chanFreq = vb.getFrequencies(0);
1144 0 : index = image.coordinates().findCoordinate(Coordinate::SPECTRAL);
1145 0 : SpectralCoordinate SpC = cfb_cs.spectralCoordinate(index);
1146 0 : Vector<Double> refVal = SpC.referenceValue();
1147 :
1148 0 : freqHi = refVal[0];
1149 0 : fillConvFuncBuffer(*cfb_p, *cfwtb_p, nx, nx, skyIncr, convSize, convSize, freqValues, wValues, wScale,
1150 : vbPA, freqHi,
1151 : polMap, polIndexMap, vb, psScale,
1152 0 : *psTerm_p, *wTerm_p, *aTerm_p, !fillCF);
1153 0 : }
1154 : // cfb_p->show(NULL,cerr);
1155 : //cfb_p->makePersistent("test.cf");
1156 : // cfwtb_p->makePersistent("test.wtcf");
1157 :
1158 0 : } // End of loop over baselines
1159 :
1160 0 : index=coords.findCoordinate(Coordinate::SPECTRAL);
1161 0 : spCS = coords.spectralCoordinate(index);
1162 0 : Vector<Double> refValue; refValue.resize(1);refValue(0)=cfRefFreq;
1163 0 : spCS.setReferenceValue(refValue);
1164 0 : coords.replaceCoordinate(spCS,index);
1165 :
1166 : // cfs.coordSys=coords; cfwts.coordSys=coords;
1167 : // cfs.pa=paQuant; cfwts.pa=paQuant;
1168 :
1169 : // aTerm_p->makeConvFunction(image,vb,wConvSize,pa,cfs,cfwts);
1170 0 : }
1171 : //
1172 : //----------------------------------------------------------------------
1173 : //
1174 0 : Bool AWConvFunc::setUpCFSupport(Array<Complex>& cffunc, Int& xSupport, Int& ySupport,
1175 : const Float& sampling, const Complex& peak)
1176 : {
1177 : //
1178 : // Find the convolution function support size. No assumption
1179 : // about the symmetry of the conv. func. can be made (except that
1180 : // they are same for all poln. planes).
1181 : //
1182 0 : xSupport = ySupport = -1;
1183 0 : Int convFuncOrigin=cffunc.shape()[0]/2, R;
1184 0 : Bool found=false;
1185 : Float threshold;
1186 : // Threshold as a fraction of the peak (presumed to be the center pixel).
1187 0 : if (abs(peak) != 0) threshold = real(abs(peak));
1188 : else
1189 0 : threshold = real(abs(cffunc(IPosition(4,convFuncOrigin,convFuncOrigin,0,0))));
1190 :
1191 : //threshold *= aTerm_p->getSupportThreshold();
1192 0 : threshold *= 1e-3;
1193 : //threshold *= 7.5e-2;
1194 :
1195 : // threshold *= 0.1;
1196 : // if (aTerm_p->isNoOp())
1197 : // threshold *= 1e-3; // This is the threshold used in "standard" FTMchines
1198 : // else
1199 :
1200 : //
1201 : // Find the support size of the conv. function in pixels
1202 : //
1203 : // Timer tim;
1204 : // tim.mark();
1205 0 : if ((found = AWConvFunc::awFindSupport(cffunc,threshold,convFuncOrigin,R)))
1206 0 : xSupport=ySupport=Int(0.5+Float(R)/sampling)+1;
1207 : // tim.show("findSupport:");
1208 :
1209 : // If the support size overflows, give a warning and set the
1210 : // support size to be convFuncSize/2 + the max. possible offset in
1211 : // the oversamplied grid. The max. possible offset would 0.5
1212 : // pixels on the sky, which would be sampling/2.0.
1213 : //
1214 : // If the extra buffer (max(offset)) is not included, the problem
1215 : // will show up when gridding the data or weights. It will not
1216 : // show up when making the avgPB since the gridding for that is
1217 : // always centered on the center of the image.
1218 0 : if ((xSupport*sampling + int(sampling/2.0+0.5)) > convFuncOrigin)
1219 : {
1220 0 : LogIO log_l(LogOrigin("AWConvFunc2", "setUpCFSupport[R&D]"));
1221 :
1222 : log_l << "Convolution function support size > N/2. Limiting it to N/2 "
1223 : << "(threshold = " << threshold << ")."
1224 0 : << LogIO::WARN;
1225 0 : xSupport = ySupport = (Int)(convFuncOrigin/sampling-1);
1226 0 : }
1227 :
1228 0 : if(xSupport<1)
1229 : {
1230 0 : LogIO log_l(LogOrigin("AWConvFunc2", "setUpCFSupport[R&D]"));
1231 :
1232 : log_l << "Convolution function is misbehaved - support seems to be zero"
1233 0 : << LogIO::EXCEPTION;
1234 0 : }
1235 0 : return found;
1236 : }
1237 : //
1238 : //----------------------------------------------------------------------
1239 : //
1240 0 : Bool AWConvFunc::resizeCF(Array<Complex>& func, Int& xSupport, Int& ySupport,
1241 : const Int& supportBuffer, const Float& sampling, const Complex& peak)
1242 : {
1243 : //LogIO log_l(LogOrigin("AWConvFunc2", "resizeCF[R&D]"));
1244 0 : Int ConvFuncOrigin=func.shape()[0]/2; // Conv. Func. is half that size of convSize
1245 :
1246 0 : Bool found = setUpCFSupport(func, xSupport, ySupport, sampling,peak);
1247 :
1248 : //Int supportBuffer = (Int)(aTerm_p->getOversampling()*1.5);
1249 : ///Make the cutout have even number of pixels...odd numbers are a pest !
1250 0 : Int bot=(Int)((ConvFuncOrigin-sampling*xSupport-supportBuffer)/2)*2; //-convSampling/2,
1251 0 : Int top=(Int)((ConvFuncOrigin+sampling*xSupport+supportBuffer)/2)*2-1; //+convSampling/2;
1252 : // bot *= 2; top *= 2;
1253 0 : bot = max(0,bot);
1254 0 : top = min(top, func.shape()(0)-1);
1255 :
1256 0 : Array<Complex> tmp;
1257 0 : IPosition blc(4,bot,bot,0,0), trc(4,top,top,0,0);
1258 : //
1259 : // Cut out the conv. func., copy in a temp. array, resize the
1260 : // CFStore.data, and copy the cutout version to CFStore.data.
1261 : //
1262 0 : tmp = func(blc,trc);
1263 0 : func.resize(tmp.shape());
1264 0 : func = tmp;
1265 0 : return found;
1266 0 : }
1267 : //
1268 : //----------------------------------------------------------------------
1269 : // A global method for use in OMP'ed findSupport() below
1270 : //
1271 0 : void archPeak(const Float& threshold, const Int& origin, const Block<Int>& cfShape, const Complex* funcPtr,
1272 : const Int& nCFS, const Int& PixInc,const Int& th, const Int& R, Block<Int>& maxR)
1273 : {
1274 0 : Block<Complex> vals;
1275 0 : Block<Int> ndx(nCFS); ndx=0;
1276 : Int NSteps;
1277 : //Check every PixInc pixel along a circle of radius R
1278 0 : NSteps = 90*R/PixInc;
1279 0 : vals.resize((Int)(NSteps+0.5));
1280 0 : uInt valsNelements=vals.nelements();
1281 0 : vals=0;
1282 :
1283 0 : for(Int pix=0;pix<NSteps;pix++)
1284 : {
1285 0 : ndx[0]=(int)(origin + R*sin(2.0*M_PI*pix*PixInc/R));
1286 0 : ndx[1]=(int)(origin + R*cos(2.0*M_PI*pix*PixInc/R));
1287 :
1288 0 : if ((ndx[0] < cfShape[0]) && (ndx[1] < cfShape[1]))
1289 : //vals[pix]=func(ndx);
1290 0 : vals[pix]=funcPtr[ndx[0]+ndx[1]*cfShape[1]+ndx[2]*cfShape[2]+ndx[3]*cfShape[3]];
1291 : }
1292 :
1293 0 : maxR[th]=-R;
1294 0 : for (uInt i=0;i<valsNelements;i++)
1295 0 : if (fabs(vals[i]) > threshold)
1296 : {
1297 0 : maxR[th]=R;
1298 0 : break;
1299 : }
1300 : // th++;
1301 0 : }
1302 : //
1303 : //----------------------------------------------------------------------
1304 : //
1305 0 : Bool AWConvFunc::findSupport(Array<Complex>& func, Float& threshold,
1306 : Int& origin, Int& radius)
1307 : {
1308 0 : return awFindSupport(func, threshold, origin, radius);
1309 : }
1310 0 : Bool AWConvFunc::awFindSupport(Array<Complex>& func, Float& threshold,
1311 : Int& origin, Int& radius)
1312 : {
1313 : //LogIO log_l(LogOrigin("AWConvFunc2", "findSupport[R&D]"));
1314 :
1315 0 : Int nCFS=func.shape().nelements(),
1316 0 : PixInc=1, R0, R1, R, convSize;
1317 0 : Block<Int> cfShape(nCFS);
1318 0 : Bool found=false;
1319 : Complex *funcPtr;
1320 : Bool dummy;
1321 0 : uInt Nth=1, threadID=0;
1322 :
1323 0 : for (Int i=0;i<nCFS;i++)
1324 0 : cfShape[i]=func.shape()[i];
1325 0 : convSize = cfShape[0];
1326 :
1327 : #ifdef _OPENMP
1328 0 : Nth = max(omp_get_max_threads()-2,1);
1329 : #endif
1330 :
1331 0 : Block<Int> maxR(Nth);
1332 :
1333 0 : funcPtr = func.getStorage(dummy);
1334 :
1335 0 : R1 = convSize/2-2;
1336 :
1337 0 : while (R1 > 1)
1338 : {
1339 0 : R0 = R1; R1 -= Nth;
1340 :
1341 : //#pragma omp parallel default(none) firstprivate(R0,R1) private(R,threadID) shared(origin,threshold,PixInc,maxR,cfShape,nCFS,funcPtr) num_threads(Nth)
1342 0 : #pragma omp parallel firstprivate(R0,R1) private(R,threadID) shared(PixInc,maxR,cfShape,nCFS,funcPtr) num_threads(Nth)
1343 : {
1344 : #pragma omp for
1345 : for(R=R0;R>R1;R--)
1346 : {
1347 : #ifdef _OPENMP
1348 : threadID=omp_get_thread_num();
1349 : #endif
1350 : archPeak(threshold, origin, cfShape, funcPtr, nCFS, PixInc, threadID, R, maxR);
1351 : }
1352 : }///omp
1353 :
1354 0 : for (uInt th=0;th<Nth;th++)
1355 0 : if (maxR[th] > 0)
1356 0 : {found=true; radius=maxR[th]; return found;}
1357 : }
1358 0 : return found;
1359 0 : }
1360 : //
1361 : //----------------------------------------------------------------------
1362 : //
1363 : // Bool AWConvFunc::findSupport(Array<Complex>& func, Float& threshold,
1364 : // Int& origin, Int& R)
1365 : // {
1366 : // LogIO log_l(LogOrigin("AWConvFunc2", "findSupport[R&D]"));
1367 : // Double NSteps;
1368 : // Int PixInc=1;
1369 : // Vector<Complex> vals;
1370 : // IPosition ndx(4,origin,0,0,0);
1371 : // Bool found=false;
1372 : // IPosition cfShape=func.shape();
1373 : // Int convSize = cfShape(0);
1374 :
1375 : // for(R=convSize/2-2;R>1;R--)
1376 : // {
1377 : // //Check every PixInc pixel along a circle of radius R
1378 : // NSteps = 90*R/PixInc;
1379 : // vals.resize((Int)(NSteps+0.5));
1380 : // vals=0;
1381 : // for(Int th=0;th<NSteps;th++)
1382 : // {
1383 : // ndx(0)=(int)(origin + R*sin(2.0*M_PI*th*PixInc/R));
1384 : // ndx(1)=(int)(origin + R*cos(2.0*M_PI*th*PixInc/R));
1385 :
1386 : // if ((ndx(0) < cfShape(0)) && (ndx(1) < cfShape(1)))
1387 : // vals(th)=func(ndx);
1388 : // }
1389 :
1390 : // if (max(abs(vals)) > threshold)
1391 : // {found=true;break;}
1392 : // }
1393 : // return found;
1394 : // }
1395 : //
1396 : //----------------------------------------------------------------------
1397 : //
1398 0 : Bool AWConvFunc::makeAverageResponse(const VisBuffer2& vb,
1399 : const ImageInterface<Complex>& image,
1400 : ImageInterface<Float>& theavgPB,
1401 : Bool reset)
1402 : {
1403 0 : TempImage<Complex> complexPB;
1404 : Bool pbMade;
1405 0 : pbMade = makeAverageResponse(vb, image, complexPB,reset);
1406 0 : normalizeAvgPB(complexPB, theavgPB);
1407 0 : return pbMade;
1408 0 : }
1409 : //
1410 : //----------------------------------------------------------------------
1411 : //
1412 0 : Bool AWConvFunc::makeAverageResponse(const VisBuffer2& vb,
1413 : const ImageInterface<Complex>& image,
1414 : ImageInterface<Complex>& theavgPB,
1415 : Bool reset)
1416 : {
1417 0 : LogIO log_l(LogOrigin("AWConvFunc2","makeAverageResponse(Complex)[R&D]"));
1418 :
1419 0 : log_l << "Making the average response for " << aTerm_p->name()
1420 0 : << LogIO::NORMAL << LogIO::POST;
1421 :
1422 0 : if (reset)
1423 : {
1424 : log_l << "Initializing the average PBs"
1425 0 : << LogIO::NORMAL << LogIO::POST;
1426 0 : theavgPB.resize(image.shape());
1427 0 : theavgPB.setCoordinateInfo(image.coordinates());
1428 0 : theavgPB.set(1.0);
1429 : }
1430 :
1431 0 : aTerm_p->applySky(theavgPB, vb, true, 0);
1432 :
1433 0 : return true; // i.e., an average PB was made
1434 0 : }
1435 : //
1436 : //----------------------------------------------------------------------
1437 : //
1438 0 : void AWConvFunc::normalizeAvgPB(ImageInterface<Complex>& inImage,
1439 : ImageInterface<Float>& outImage)
1440 : {
1441 0 : LogIO log_l(LogOrigin("AWConvFunc2", "normalizeAvgPB[R&D]"));
1442 :
1443 0 : String name("avgpb.im");
1444 0 : storeImg(name,inImage);
1445 0 : IPosition inShape(inImage.shape()),ndx(4,0,0,0,0);
1446 0 : Vector<Complex> peak(inShape(2));
1447 :
1448 0 : outImage.resize(inShape);
1449 0 : outImage.setCoordinateInfo(inImage.coordinates());
1450 :
1451 : Bool isRefIn;
1452 0 : Array<Complex> inBuf;
1453 0 : Array<Float> outBuf;
1454 :
1455 0 : isRefIn = inImage.get(inBuf);
1456 : //isRefOut = outImage.get(outBuf);
1457 : log_l << "Normalizing the average PBs to unity"
1458 0 : << LogIO::NORMAL << LogIO::POST;
1459 : //
1460 : // Normalize each plane of the inImage separately to unity.
1461 : //
1462 0 : Complex inMax = max(inBuf);
1463 0 : if (abs(inMax)-1.0 > 1E-3)
1464 : {
1465 0 : for(ndx(3)=0;ndx(3)<inShape(3);ndx(3)++)
1466 0 : for(ndx(2)=0;ndx(2)<inShape(2);ndx(2)++)
1467 : {
1468 0 : peak(ndx(2)) = 0;
1469 0 : for(ndx(1)=0;ndx(1)<inShape(1);ndx(1)++)
1470 0 : for(ndx(0)=0;ndx(0)<inShape(0);ndx(0)++)
1471 0 : if (abs(inBuf(ndx)) > peak(ndx(2)))
1472 0 : peak(ndx(2)) = inBuf(ndx);
1473 :
1474 0 : for(ndx(1)=0;ndx(1)<inShape(1);ndx(1)++)
1475 0 : for(ndx(0)=0;ndx(0)<inShape(0);ndx(0)++)
1476 : // avgPBBuf(ndx) *= (pbPeaks(ndx(2))/peak(ndx(2)));
1477 0 : inBuf(ndx) /= peak(ndx(2));
1478 : }
1479 0 : if (isRefIn) inImage.put(inBuf);
1480 : }
1481 :
1482 0 : ndx=0;
1483 0 : for(ndx(1)=0;ndx(1)<inShape(1);ndx(1)++)
1484 0 : for(ndx(0)=0;ndx(0)<inShape(0);ndx(0)++)
1485 : {
1486 0 : IPosition plane1(ndx);
1487 0 : plane1=ndx;
1488 0 : plane1(2)=1; // The other poln. plane
1489 : // avgPBBuf(ndx) = (avgPBBuf(ndx) + avgPBBuf(plane1))/2.0;
1490 0 : outBuf(ndx) = sqrt(real(inBuf(ndx) * inBuf(plane1)));
1491 0 : }
1492 : //
1493 : // Rather convoluted way of copying Pol. plane-0 to Pol. plane-1!!!
1494 : //
1495 0 : for(ndx(1)=0;ndx(1)<inShape(1);ndx(1)++)
1496 0 : for(ndx(0)=0;ndx(0)<inShape(0);ndx(0)++)
1497 : {
1498 0 : IPosition plane1(ndx);
1499 0 : plane1=ndx;
1500 0 : plane1(2)=1; // The other poln. plane
1501 0 : outBuf(plane1) = real(outBuf(ndx));
1502 0 : }
1503 0 : }
1504 : //
1505 : //-------------------------------------------------------------------------
1506 : // Legacy code. Should ultimately be deleteted after re-facatoring
1507 : // is finished.
1508 : //
1509 0 : Bool AWConvFunc::makeAverageResponse_org(const VisBuffer2& vb,
1510 : const ImageInterface<Complex>& image,
1511 : ImageInterface<Float>& theavgPB,
1512 : Bool reset)
1513 : {
1514 0 : LogIO log_l(LogOrigin("AWConvFunc2", "makeAverageResponse_org[R&D]"));
1515 0 : TempImage<Float> localPB;
1516 :
1517 : log_l << "Making the average response for "
1518 0 : << aTerm_p->name()
1519 0 : << LogIO::NORMAL << LogIO::POST;
1520 :
1521 0 : localPB.resize(image.shape()); localPB.setCoordinateInfo(image.coordinates());
1522 0 : if (reset)
1523 : {
1524 0 : log_l << "Initializing the average PBs" << LogIO::NORMAL << LogIO::POST;
1525 0 : theavgPB.resize(localPB.shape());
1526 0 : theavgPB.setCoordinateInfo(localPB.coordinates());
1527 0 : theavgPB.set(0.0);
1528 : }
1529 : //
1530 : // Make the Stokes PB
1531 : //
1532 0 : localPB.set(1.0);
1533 :
1534 : // Block<CountedPtr<ImageInterface<Float > > > tmpBlock(1);
1535 : // tmpBlock[0]=CountedPtr<ImageInterface<Float> >(&localPB, false);
1536 : // aTerm_p->applySky(tmpBlock, vb, 0, false);
1537 0 : aTerm_p->applySky(localPB, vb, false, 0);
1538 :
1539 0 : IPosition twoDPBShape(localPB.shape());
1540 0 : TempImage<Complex> localTwoDPB(twoDPBShape,localPB.coordinates());
1541 : // localTwoDPB.setMaximumCacheSize(cachesize);
1542 : Int NAnt;
1543 0 : NAnt=1;
1544 :
1545 0 : for(Int ant=0;ant<NAnt;ant++)
1546 : { //Ant loop
1547 : {
1548 0 : IPosition ndx(4,0,0,0,0);
1549 0 : for(ndx(0)=0; ndx(0)<twoDPBShape(0); ndx(0)++)
1550 0 : for(ndx(1)=0; ndx(1)<twoDPBShape(1); ndx(1)++)
1551 0 : for(ndx(2)=0; ndx(2)<twoDPBShape(2); ndx(2)++)
1552 0 : for(ndx(3)=0; ndx(3)<twoDPBShape(3); ndx(3)++)
1553 0 : localTwoDPB.putAt(Complex((localPB(ndx)),0.0),ndx);
1554 0 : }
1555 : //
1556 : // Accumulate the shifted PBs
1557 : //
1558 : {
1559 : Bool isRefF;
1560 0 : Array<Float> fbuf;
1561 0 : Array<Complex> cbuf;
1562 0 : isRefF=theavgPB.get(fbuf);
1563 : //isRefC=localTwoDPB.get(cbuf);
1564 :
1565 0 : IPosition fs(fbuf.shape());
1566 0 : IPosition ndx(4,0,0,0,0),avgNDX(4,0,0,0,0);
1567 0 : for(ndx(3)=0,avgNDX(3)=0;ndx(3)<fs(3);ndx(3)++,avgNDX(3)++)
1568 0 : for(ndx(2)=0,avgNDX(2)=0;ndx(2)<twoDPBShape(2);ndx(2)++,avgNDX(2)++)
1569 0 : for(ndx(0)=0,avgNDX(0)=0;ndx(0)<fs(0);ndx(0)++,avgNDX(0)++)
1570 0 : for(ndx(1)=0,avgNDX(1)=0;ndx(1)<fs(1);ndx(1)++,avgNDX(1)++)
1571 : {
1572 : Float val;
1573 0 : val = real(cbuf(ndx));
1574 0 : fbuf(avgNDX) += val;
1575 : }
1576 0 : if (!isRefF) theavgPB.put(fbuf);
1577 0 : }
1578 : }
1579 0 : theavgPB.setCoordinateInfo(localPB.coordinates());
1580 0 : return true; // i.e., an average PB was made
1581 0 : }
1582 : //
1583 : //----------------------------------------------------------------------
1584 : //
1585 : // void AWConvFunc::prepareConvFunction(const VisBuffer2& vb, VBRow2CFBMapType& theMap)
1586 0 : void AWConvFunc::prepareConvFunction(const VisBuffer2& vb, VB2CFBMap& theMap)
1587 : {
1588 0 : if (aTerm_p->rotationallySymmetric() == false) return;
1589 0 : Int nRow=theMap.nelements();
1590 : // CountedPtr<CFBuffer> cfb, cbPtr;
1591 : // CountedPtr<CFCell> cfc;
1592 : // CountedPtr<ATerm> aTerm_l=aTerm_p;
1593 0 : CFBuffer *cfb, *cbPtr=0;
1594 0 : CFCell *cfc, *baseCFC=NULL;
1595 0 : ATerm *aTerm_l=&*aTerm_p;
1596 :
1597 0 : cfb=&*(theMap[0]);
1598 0 : cfc = &*(cfb->getCFCellPtr(0,0,0));
1599 0 : Double actualPA = getPA(vb), currentCFPA = cfc->pa_p.getValue("rad");
1600 0 : Double dPA = currentCFPA-actualPA;
1601 :
1602 0 : if (fabs(dPA) <= fabs(rotateCFOTFAngleRad_p)) return;
1603 :
1604 :
1605 : // Int Nth=1;
1606 : // #ifdef _OPENMP
1607 : // Nth=max(omp_get_max_threads()-2,1);
1608 : // #endif
1609 0 : for (Int irow=0;irow<nRow;irow++)
1610 : {
1611 0 : cfb=&*(theMap[irow]);
1612 : // if ((!cfb.null()) && (cfb != cbPtr))
1613 0 : if ((cfb!=NULL) && (cfb != cbPtr))
1614 : {
1615 : // baseCFB_p = cfb->clone();
1616 : // cerr << "NRef = " << baseCFB_p.nrefs() << endl;
1617 : //
1618 : // If the following messsage is emitted more than once, we
1619 : // are in a heterogeneous-array case
1620 : //
1621 0 : LogIO log_l(LogOrigin("AWConvFunc2", "prepareConvFunction"));
1622 0 : log_l << "Rotating the base CFB from PA=" << cfb->getCFCellPtr(0,0,0)->pa_p.getValue("deg")
1623 : << " to " << actualPA*57.2957795131
1624 0 : << " " << cfb->getCFCellPtr(0,0,0)->shape_p
1625 0 : << LogIO::DEBUG1 << LogIO::POST;
1626 :
1627 0 : IPosition shp(cfb->shape());
1628 0 : cbPtr = cfb;
1629 0 : for(Int k=0;k<shp(2);k++) // Mueller-loop
1630 0 : for(Int j=0;j<shp(1);j++) // W-loop
1631 : // #pragma omp parallel default(none) firstprivate(j,k) shared(shp,cfb,aTerm_l) num_threads(Nth)
1632 : {
1633 : // #pragma omp for
1634 0 : for (Int i=0;i<shp(0);i++) // Chan-loop
1635 : {
1636 0 : cfc = &*(cfb->getCFCellPtr(i,j,k));
1637 : //baseCFC = &*(baseCFB_p->getCFCellPtr(i,j,k));
1638 : // Call this for every VB. Any optimization
1639 : // (e.g. rotating at some increment only) is
1640 : // implemented in the ATerm::rotate().
1641 : // if (rotateCF_p)
1642 : // Rotate the cell only if it has been loaded.
1643 0 : if (cfc->getShape().product() > 0)
1644 0 : aTerm_l->rotate2(vb,*baseCFC, *cfc,rotateCFOTFAngleRad_p);
1645 : }
1646 : }
1647 0 : }
1648 : }
1649 : };
1650 : //
1651 : //----------------------------------------------------------------------
1652 : //
1653 0 : void AWConvFunc::setMiscInfo(const RecordInterface& params)
1654 : {
1655 : (void)params;
1656 0 : }
1657 : //
1658 : // REFACTORED CODE
1659 : //
1660 :
1661 : //
1662 : //----------------------------------------------------------------------
1663 : //
1664 0 : void AWConvFunc::fillConvFuncBuffer2(CFBuffer& cfb, CFBuffer& cfWtb,
1665 : const Int& nx, const Int& ny,
1666 : const ImageInterface<Complex>& skyImage,
1667 : const CFCStruct& miscInfo,
1668 : PSTerm& psTerm, WTerm& wTerm, ATerm& aTerm,
1669 : Bool conjBeams)
1670 :
1671 : {
1672 0 : LogIO log_l(LogOrigin("Moluscs", "fillConvFuncBuffer2[R&D]"));
1673 0 : Complex cfNorm, cfWtNorm;
1674 0 : Complex cpeak;
1675 : {
1676 : Float sampling, samplingWt;
1677 : Int xSupport, ySupport, xSupportWt, ySupportWt;
1678 0 : CoordinateSystem cs_l;
1679 0 : String bandName;
1680 : // Extract the parameters index by (MuellerElement, Freq, W)
1681 0 : cfWtb.getParams(cs_l, samplingWt, xSupportWt, ySupportWt, bandName,
1682 0 : miscInfo.freqValue, miscInfo.wValue, //The address of CFCell as physical co-ords
1683 0 : miscInfo.muellerElement);
1684 0 : cfb.getParams(cs_l, sampling, xSupport, ySupport, bandName,
1685 0 : miscInfo.freqValue,miscInfo.wValue, //The address of CFCell as physical co-ords
1686 0 : miscInfo.muellerElement);
1687 :
1688 : //cerr << "FCFB2: frq " << miscInfo.freqValue << " cs_l " << cs_l.toWorld(IPosition(4, 0,0,0,0)) << endl;
1689 0 : aTerm.setBandName(bandName);
1690 : //
1691 : // Cache the A-Term for this polarization and frequency
1692 : //
1693 : Double conjFreq, vbPA;
1694 0 : CountedPtr<CFCell> thisCell=cfb.getCFCellPtr(miscInfo.freqValue, miscInfo.wValue, miscInfo.muellerElement);
1695 0 : vbPA = thisCell->pa_p.getValue("rad");
1696 0 : conjFreq = thisCell->conjFreq_p;
1697 0 : CoordinateSystem conjPolCS_l=cs_l; AWConvFunc::makeConjPolAxis(conjPolCS_l, thisCell->conjPoln_p);
1698 0 : IPosition pbshp(4,nx,ny,1,1);
1699 0 : TempImage<Complex> ftATerm_l(pbshp, cs_l), ftATermSq_l(pbshp,conjPolCS_l);
1700 0 : Bool doSquint=true;
1701 0 : ftATerm_l.set(Complex(1.0,0.0)); ftATermSq_l.set(Complex(1.0,0.0));
1702 0 : Double freq_l=miscInfo.freqValue;
1703 : //TESTOO
1704 : {
1705 0 : Vector<String> csList;
1706 0 : IPosition dummy;
1707 : // cout << "CoordSys:===================== ";
1708 : // // csList = ftATermSq_l.coordinates().list(log_l,MDoppler::RADIO,dummy,dummy);
1709 :
1710 0 : csList = cs_l.list(log_l,MDoppler::RADIO,dummy,dummy);
1711 : //cerr << csList << endl;
1712 : // csList = conjPolCS_l.list(log_l,MDoppler::RADIO,dummy,dummy);
1713 : // cout << csList << endl;
1714 0 : }
1715 :
1716 : //if (!isDryRun)
1717 : //cerr <<"applying ATERM for " << freq_l << endl;
1718 : //TESTOO
1719 0 : CoordinateSystem lalacs=cs_l;
1720 : //
1721 :
1722 : // cerr << "#########$$$$$$ " << pbshp << " " << nx << " " << freq_l << " " << conjFreq << endl;
1723 : {
1724 :
1725 0 : aTerm.applySky(ftATerm_l, vbPA, doSquint, 0, miscInfo.muellerElement,freq_l);//freqHi);
1726 0 : if (conjBeams) aTerm.applySky(ftATermSq_l, vbPA, doSquint, 0, miscInfo.muellerElement, conjFreq);//freqHi);
1727 0 : else aTerm.applySky(ftATermSq_l, vbPA, doSquint, 0,miscInfo.muellerElement,freq_l);
1728 : }
1729 :
1730 : //TESTOO
1731 : {
1732 0 : PagedImage<Complex> lala(ftATerm_l.shape(), lalacs, "ATERMFCFB2.im");
1733 0 : lala.copyData(ftATerm_l);
1734 :
1735 0 : }
1736 : ////
1737 :
1738 0 : Vector<Double> cellSize;
1739 : // {
1740 : // Int linIndex=cs_l.findCoordinate(Coordinate::LINEAR);
1741 : // LinearCoordinate lc=cs_l.linearCoordinate(linIndex);
1742 : // Vector<Bool> axes(2); axes=true;
1743 : // Vector<Int> dirShape(2); dirShape(0)=nx;dirShape(1)=ny;
1744 : // Coordinate* FTlc=lc.makeFourierCoordinate(axes,dirShape);
1745 : // cellSize = lc.increment();
1746 : // }
1747 : {
1748 0 : CoordinateSystem skyCoords(skyImage.coordinates());
1749 :
1750 0 : Int directionIndex=skyCoords.findCoordinate(Coordinate::DIRECTION);
1751 0 : DirectionCoordinate dc=skyCoords.directionCoordinate(directionIndex);
1752 : //Vector<Double> cellSize;
1753 0 : cellSize = dc.increment()*(Double)(miscInfo.sampling*skyImage.shape()(0)/nx); // nx is the size of the CF buffer
1754 0 : }
1755 : //cerr << "#########$$$$$$ " << cellSize << endl;
1756 :
1757 : // Int directionIndex=cs_l.findCoordinate(Coordinate::DIRECTION);
1758 : // DirectionCoordinate dc=cs_l.directionCoordinate(directionIndex);
1759 : // cellSize = dc.increment();
1760 :
1761 : //
1762 : // Now compute the PS x W-Term and apply the cached
1763 : // A-Term to build the full CF.
1764 : //
1765 : {
1766 : log_l << " CF("
1767 0 : << "M:"<< miscInfo.muellerElement
1768 0 : << ",C:" << miscInfo.freqValue/1e9
1769 0 : << ",W:" << miscInfo.wValue << "): ";
1770 0 : Array<Complex> &cfWtBuf=(*(cfWtb.getCFCellPtr(miscInfo.freqValue, miscInfo.wValue, miscInfo.muellerElement))->storage_p);
1771 0 : Array<Complex> &cfBuf=(*(cfb.getCFCellPtr(miscInfo.freqValue, miscInfo.wValue, miscInfo.muellerElement))->storage_p);
1772 :
1773 0 : cfWtBuf.resize(pbshp);
1774 0 : cfBuf.resize(pbshp);
1775 :
1776 0 : const Vector<Double> sampling_l(2,sampling);
1777 0 : Matrix<Complex> cfBufMat(cfBuf.nonDegenerate()),
1778 0 : cfWtBufMat(cfWtBuf.nonDegenerate());
1779 : //
1780 : // Apply the Prolate Spheroidal and W-Term kernels
1781 : //
1782 0 : Vector<Double> s(2); s=sampling;
1783 : //Timer tim;
1784 : //tim.mark();
1785 : // if (psTerm.isNoOp() || isDryRun)
1786 0 : if (psTerm.isNoOp())
1787 0 : cfBufMat = cfWtBufMat = 1.0;
1788 : else
1789 : {
1790 : // psTerm.applySky(cfBufMat, False); // Assign (psScale set in psTerm.init()
1791 : // psTerm.applySky(cfWtBufMat, False); // Assign
1792 0 : psTerm.applySky(cfBufMat, s, cfBufMat.shape()(0)/s(0)); // Assign (psScale set in psTerm.init()
1793 0 : psTerm.applySky(cfWtBufMat, s, cfWtBufMat.shape()(0)/s(0)); // Assign
1794 0 : cfWtBuf *= cfWtBuf;
1795 : }
1796 :
1797 : //tim.mark();
1798 0 : if (miscInfo.wValue > 0)
1799 0 : wTerm.applySky(cfBufMat, cellSize, miscInfo.wValue, cfBuf.shape()(0));///4);
1800 :
1801 0 : IPosition PolnPlane(4,0,0,0,0),
1802 0 : pbShape(4, cfBuf.shape()(0), cfBuf.shape()(1), 1, 1);
1803 : //
1804 : // Make TempImages and copy the buffers with PS *
1805 : // WKernel applied (too bad that TempImages can't be
1806 : // made with existing buffers)
1807 : //
1808 : //-------------------------------------------------------------
1809 0 : TempImage<Complex> twoDPB_l(pbShape, cs_l);
1810 0 : TempImage<Complex> twoDPBSq_l(pbShape,cs_l);
1811 : //-------------------------------------------------------------
1812 : // WBAWP CODE BEGIN -- ftATermSq_l has conj. PolCS
1813 :
1814 0 : cfWtBuf *= ftATerm_l.get()*conj(ftATermSq_l.get());
1815 :
1816 : //tim.mark();
1817 0 : cfBuf *= ftATerm_l.get();
1818 : //tim.show("W*A*2: ");
1819 : // WBAWP CODE END
1820 : //tim.mark();
1821 0 : twoDPB_l.putSlice(cfBuf, PolnPlane);
1822 0 : twoDPBSq_l.putSlice(cfWtBuf, PolnPlane);
1823 : //tim.show("putSlice:");
1824 :
1825 : // To accumulate avgPB2, call this function.
1826 : // PBSQWeight
1827 : // Bool PBSQ = false;
1828 : // if(PBSQ) makePBSq(twoDPBSq_l);
1829 :
1830 : //
1831 : // Set the ref. freq. of the co-ordinate system to
1832 : // that set by ATerm::applySky().
1833 : //
1834 : //tim.mark();
1835 0 : CoordinateSystem cs=twoDPB_l.coordinates();
1836 0 : Int index= twoDPB_l.coordinates().findCoordinate(Coordinate::SPECTRAL);
1837 0 : SpectralCoordinate SpCS = twoDPB_l.coordinates().spectralCoordinate(index);
1838 :
1839 0 : Double cfRefFreq=SpCS.referenceValue()(0);
1840 0 : Vector<Double> refValue; refValue.resize(1); refValue(0)=cfRefFreq;
1841 0 : SpCS.setReferenceValue(refValue);
1842 0 : cs.replaceCoordinate(SpCS,index);
1843 :
1844 : //tim.mark();
1845 : // if (!isDryRun)
1846 : {
1847 0 : LatticeFFT::cfft2d(twoDPB_l);
1848 0 : LatticeFFT::cfft2d(twoDPBSq_l);
1849 : }
1850 : //tim.show("FFT*2:");
1851 :
1852 : //tim.mark();
1853 0 : IPosition shp(twoDPB_l.shape());
1854 0 : IPosition start(4, 0, 0, 0, 0), pbSlice(4, shp[0]-1, shp[1]-1,1/*polInUse*/, 1),
1855 0 : sliceLength(4,cfBuf.shape()[0]-1,cfBuf.shape()[1]-1,1,1);
1856 :
1857 0 : cfBuf(Slicer(start,sliceLength)).nonDegenerate()
1858 0 : =(twoDPB_l.getSlice(start, pbSlice, true));
1859 :
1860 0 : shp = twoDPBSq_l.shape();
1861 0 : IPosition pbSqSlice(4, shp[0]-1, shp[1]-1, 1, 1),
1862 0 : sqSliceLength(4,cfWtBuf.shape()(0)-1,cfWtBuf.shape()[1]-1,1,1);
1863 :
1864 0 : cfWtBuf(Slicer(start,sqSliceLength)).nonDegenerate()
1865 0 : =(twoDPBSq_l.getSlice(start, pbSqSlice, true));
1866 : //tim.show("Slicer*2:");
1867 : //
1868 : //tim.mark();
1869 : // if (!isDryRun)
1870 : // {
1871 : // if (wValue==0) wtcpeak = max(cfWtBuf);
1872 : // cfWtBuf /= wtcpeak;
1873 : // }
1874 : //tim.show("Norm");
1875 :
1876 : //tim.mark();
1877 : // if (!isDryRun)
1878 0 : Int supportBuffer = (Int)(AWConvFunc::getOversampling(psTerm, wTerm, aTerm)*2.0);
1879 0 : AWConvFunc::resizeCF(cfWtBuf, xSupportWt, ySupportWt, supportBuffer, samplingWt,0.0);
1880 : //tim.show("Resize:");
1881 :
1882 : //tim.mark();
1883 0 : Vector<Double> ftRef(2);
1884 0 : ftRef(0)=cfWtBuf.shape()(0)/2.0;
1885 0 : ftRef(1)=cfWtBuf.shape()(1)/2.0;
1886 0 : CoordinateSystem ftCoords=cs_l;
1887 0 : SynthesisUtils::makeFTCoordSys(cs_l, cfWtBuf.shape()(0), ftRef, ftCoords);
1888 :
1889 0 : thisCell=cfWtb.getCFCellPtr(miscInfo.freqValue, miscInfo.wValue, miscInfo.muellerElement);
1890 0 : thisCell->pa_p=Quantity(vbPA,"rad");
1891 0 : thisCell->coordSys_p = ftCoords;
1892 0 : thisCell->xSupport_p = xSupportWt;
1893 0 : thisCell->ySupport_p = ySupportWt;
1894 0 : thisCell->isRotationallySymmetric_p = aTerm.isNoOp();
1895 : //tim.show("CSStuff:");
1896 :
1897 : //tim.mark();
1898 : // if (!isDryRun)
1899 : {
1900 0 : cpeak = max(cfBuf);
1901 0 : cfBuf /= cpeak;
1902 : }
1903 : //tim.show("Peaknorm:");
1904 :
1905 : // if (!isDryRun)
1906 0 : AWConvFunc::resizeCF(cfBuf, xSupport, ySupport, supportBuffer, sampling,0.0);
1907 :
1908 0 : log_l << "CF Support: " << xSupport << " (" << xSupportWt << ") " << "pixels" << LogIO::POST;
1909 :
1910 0 : ftRef(0)=cfBuf.shape()(0)/2.0;
1911 0 : ftRef(1)=cfBuf.shape()(1)/2.0;
1912 :
1913 : //tim.mark();
1914 : //cfNorm=cfWtNorm=1.0;
1915 : // if ((wValue == 0) && (!isDryRun))
1916 : //if (miscInfo.wValue == 0)
1917 : {
1918 0 : cfNorm=0; cfWtNorm=0;
1919 0 : cfNorm = AWConvFunc::cfArea(cfBufMat, xSupport, ySupport, sampling);
1920 0 : cfWtNorm = AWConvFunc::cfArea(cfWtBufMat, xSupportWt, ySupportWt, sampling);
1921 : }
1922 : //tim.show("Area*2:");
1923 :
1924 : //tim.mark();
1925 0 : if (cfNorm != Complex(0.0)) cfBuf /= cfNorm;
1926 0 : if (cfWtNorm != Complex(0.0)) cfWtBuf /= cfWtNorm;
1927 : //tim.show("cfNorm*2:");
1928 :
1929 : //tim.mark();
1930 0 : ftCoords=cs_l;
1931 0 : SynthesisUtils::makeFTCoordSys(cs_l, cfBuf.shape()(0), ftRef, ftCoords);
1932 :
1933 0 : thisCell=cfb.getCFCellPtr(miscInfo.freqValue, miscInfo.wValue, miscInfo.muellerElement);
1934 0 : thisCell->pa_p=Quantity(vbPA,"rad");
1935 0 : thisCell->coordSys_p = ftCoords;
1936 0 : thisCell->xSupport_p = xSupport;
1937 0 : thisCell->ySupport_p = ySupport;
1938 0 : thisCell->isRotationallySymmetric_p = aTerm.isNoOp();
1939 0 : (cfWtb.getCFCellPtr(miscInfo.freqValue, miscInfo.wValue, miscInfo.muellerElement))->initCache();
1940 0 : (cfb.getCFCellPtr(miscInfo.freqValue, miscInfo.wValue, miscInfo.muellerElement))->initCache();
1941 : //tim.show("End*2:");
1942 0 : }
1943 0 : }
1944 0 : }
1945 :
1946 : // extern casacore::Double casa::EVLABandMinFreqDefaults[EVLABeamCalc_NumBandCodes];
1947 :
1948 : //
1949 : //----------------------------------------------------------------------
1950 : //
1951 0 : void AWConvFunc::makeConvFunction2(const String& cfCachePath,
1952 : const Vector<Double>&,// uvScale,
1953 : const Vector<Double>& uvOffset,
1954 : const Matrix<Double>& ,//vbFreqSelection,
1955 : CFStore2& cfs2,
1956 : CFStore2& cfwts2,
1957 : const Bool psTermOn,
1958 : const Bool aTermOn,
1959 : const Bool conjBeams)
1960 : {
1961 0 : LogIO log_l(LogOrigin("crustaceans", "makeConvFunction2[R&D]"));
1962 : Int convSize, convSampling;//, polInUse;
1963 0 : Array<Complex> convFunc_l, convWeights_l;
1964 : //
1965 : // Get the coordinate system
1966 : //
1967 0 : const String uvGridDiskImage=cfCachePath+"/"+"uvgrid.im";
1968 0 : PagedImage<Complex> skyImage_l(uvGridDiskImage);//cfs2.getCacheDir()+"/uvgrid.im");
1969 : Double skyMinFreq;
1970 0 : Vector<Double> skyIncr;
1971 : Int skyNX,skyNY;
1972 : {
1973 0 : skyNX=skyImage_l.shape()(0);
1974 0 : skyNY=skyImage_l.shape()(1);
1975 0 : CoordinateSystem skyCoords(skyImage_l.coordinates());
1976 0 : skyCoords.list(log_l, MDoppler::RADIO, skyImage_l.shape(), skyImage_l.shape(), True);
1977 0 : Int directionIndex=skyCoords.findCoordinate(Coordinate::DIRECTION);
1978 0 : DirectionCoordinate dc=skyCoords.directionCoordinate(directionIndex);
1979 0 : skyIncr = dc.increment();
1980 0 : }
1981 0 : CountedPtr<CFBuffer> cfb_p, cfwtb_p;
1982 :
1983 0 : IPosition cfsShape = cfs2.getShape();
1984 0 : IPosition wCFStShape = cfwts2.getShape();
1985 :
1986 : //Matrix<Int> uniqueBaselineTypeList=makeBaselineList(aTerm_p->getAntTypeList());
1987 : Bool wbAWP, wTermOn;
1988 :
1989 0 : for (int iPA=0; iPA<cfsShape[0]; iPA++)
1990 0 : for (int iB=0; iB<cfsShape[1]; iB++)
1991 : {
1992 0 : log_l << "Filling CFs for baseline type " << iB << ", PA slot " << iPA << LogIO::WARN << LogIO::POST;
1993 0 : cfb_p=cfs2.getCFBuffer(iPA,iB);
1994 0 : cfwtb_p=cfwts2.getCFBuffer(iPA,iB);
1995 :
1996 0 : IPosition cfbShape = cfb_p->shape();
1997 : //cerr << "@@@CFBSHAPE " <<cfbShape << endl;
1998 0 : for (int iNu=0; iNu<cfbShape(0); iNu++) // Frequency axis
1999 0 : for (int iPol=0; iPol<cfbShape(2); iPol++) // Polarization axis
2000 0 : for (int iW=0; iW<cfbShape(1); iW++) // W axis
2001 : {
2002 0 : CFCStruct miscInfo;
2003 0 : CoordinateSystem cs_l;
2004 : Int xSupport, ySupport;
2005 : Float sampling;
2006 :
2007 0 : CountedPtr<CFCell>& tt=(*cfb_p).getCFCellPtr(iNu, iW, iPol);
2008 : //cerr << "####@#$#@$@ " << iNu << " " << iW << " " << iPol << endl;
2009 : //tt->show("test",cout);
2010 0 : if (tt->cfShape_p.nelements() != 0)
2011 : {
2012 0 : (*cfb_p)(iNu,iW,iPol).getAsStruct(miscInfo); // Get misc. info. for this CFCell
2013 : {
2014 : //This code uses the BeamCalc class to get
2015 : //the nominal min. freq. of the band in
2016 : //use. While not accurate, may be
2017 : //sufficient for the purpose of the
2018 : //anti-aliasing operator.
2019 : ///At this stage if telescopeName is blank or empty spaces
2020 : //then it is EVLA
2021 0 : if(miscInfo.telescopeName.size() < 2)
2022 0 : miscInfo.telescopeName="EVLA";
2023 0 : Int bandID = BeamCalc::Instance()->getBandID(miscInfo.freqValue,miscInfo.telescopeName,miscInfo.bandName);
2024 0 : skyMinFreq = casa::EVLABandMinFreqDefaults[bandID];
2025 : }
2026 0 : wbAWP=True; // Always true since the Freq. value is got from the coord. sys.
2027 0 : wTermOn=(miscInfo.wValue > 0.0);
2028 :
2029 : CountedPtr<ConvolutionFunction> awCF = AWProjectFT::makeCFObject(miscInfo.telescopeName,
2030 0 : aTermOn, psTermOn, wTermOn, True, wbAWP, conjBeams);
2031 0 : (static_cast<AWConvFunc &>(*awCF)).aTerm_p->cacheVBInfo(miscInfo.telescopeName, miscInfo.diameter);
2032 : //aTerm_p->cacheVBInfo(miscInfo.telescopeName, miscInfo.diameter);
2033 :
2034 0 : String bandName;
2035 0 : cfb_p->getParams(cs_l, sampling, xSupport, ySupport,bandName,iNu,iW,iPol);
2036 0 : convSampling=miscInfo.sampling;
2037 :
2038 : //convSize=miscInfo.shape[0];
2039 : // This method loads "empty CFs". Those have
2040 : // support size equal to the CONVBUF size
2041 : // required. So use that, instead of the
2042 : // "shape" information from CFs, since the
2043 : // latter for empty CFs can be small (to save
2044 : // disk space and i/o -- the CFs are supposed
2045 : // to be empty anyway at this stage!)
2046 0 : convSize=xSupport;
2047 :
2048 0 : IPosition start(4, 0, 0, 0, 0);
2049 0 : IPosition pbSlice(4, convSize, convSize, 1, 1);
2050 :
2051 0 : Matrix<Complex> screen(convSize, convSize);
2052 :
2053 : {
2054 : // Set up the anti-aliasing operator (psTerm_p) for this CF.
2055 0 : Int inner=convSize/(convSampling);
2056 :
2057 : //Float psScale = (2*coords.increment()(0))/(nx*image.coordinates().increment()(0));
2058 0 : Float innerQuaterFraction=1.0;
2059 0 : innerQuaterFraction=refim::SynthesisUtils::getenv("AWCF.FUDGE",innerQuaterFraction);
2060 :
2061 0 : Double lambdaByD = innerQuaterFraction*1.22*C::c/skyMinFreq/miscInfo.diameter;
2062 0 : Double FoV_x = fabs(skyNX*skyIncr(0));
2063 0 : Double FoV_y = fabs(skyNY*skyIncr(1));
2064 0 : Vector<Double> uvScale_l(3);
2065 0 : uvScale_l(0) = (FoV_x < lambdaByD) ? FoV_x : lambdaByD;
2066 0 : uvScale_l(1) = (FoV_y < lambdaByD) ? FoV_y : lambdaByD;
2067 0 : uvScale_l(2) = 0.0;
2068 :
2069 0 : Float psScale = 2.0/(innerQuaterFraction*convSize/convSampling);// nx*image.coordinates().increment()(0)*convSampling/2;
2070 0 : ((static_cast<AWConvFunc &>(*awCF)).psTerm_p)->init(IPosition(2,inner,inner), uvScale_l, uvOffset,psScale);
2071 0 : }
2072 :
2073 : //
2074 : // By this point, the all the 4 axis (Time/PA, Freq, Pol,
2075 : // Baseline) of the CFBuffer objects have been setup. The CFs
2076 : // will now be filled using the supplied PS-, W- ad A-term objects.
2077 : //
2078 :
2079 0 : AWConvFunc::fillConvFuncBuffer2(*cfb_p, *cfwtb_p, convSize, convSize,
2080 : skyImage_l,
2081 : miscInfo,
2082 0 : *((static_cast<AWConvFunc &>(*awCF)).psTerm_p),
2083 0 : *((static_cast<AWConvFunc &>(*awCF)).wTerm_p),
2084 0 : *((static_cast<AWConvFunc &>(*awCF)).aTerm_p),
2085 : conjBeams);
2086 :
2087 : // *psTerm_p, *wTerm_p, *aTerm_p);
2088 : //cfb_p->show(NULL,cerr);
2089 : //
2090 : // Make the CFStores persistent.
2091 : //
2092 : // cfs2.makePersistent(cfCachePath.c_str());
2093 : // cfwts2.makePersistent(cfCachePath.c_str(),"WT");
2094 0 : }
2095 0 : }
2096 0 : } // End of loop over baselines
2097 :
2098 0 : cfs2.makePersistent(cfCachePath.c_str());
2099 0 : cfwts2.makePersistent(cfCachePath.c_str(),"","WT");
2100 : // Directory dir(uvGridDiskImage);
2101 : // dir.removeRecursive(false);
2102 : // dir.remove();
2103 0 : }
2104 0 : Int AWConvFunc::getOversampling(PSTerm& psTerm, WTerm& wTerm, ATerm& aTerm)
2105 : {
2106 : Int os;
2107 0 : if (!aTerm.isNoOp()) os=aTerm.getOversampling();
2108 0 : else if (!wTerm.isNoOp()) os=wTerm.getOversampling();
2109 0 : else os=psTerm.getOversampling();
2110 0 : return os;
2111 : }
2112 0 : void AWConvFunc::makeAConvFunc(Array<Complex>& convFunc,
2113 : Array<Complex>& wtConvFunc, CoordinateSystem& csys,
2114 : Vector<Int>& asupport, Int& npix,
2115 : const Vector<Double>& freqlist, const Bool doSquint,
2116 : const Double& pa){
2117 :
2118 : //Assuming first freq is lowest
2119 0 : Quantity cell;
2120 0 : Int convnx=0;
2121 0 : asupport.resize(freqlist.nelements());
2122 0 : if(!atermMaker_p)
2123 0 : throw(AipsError("Programmer Error: AWConvFunc has to be constructed with a EVLAAperture"));
2124 0 : String bandname=EVLAAperture::getVLABandName(freqlist[int(freqlist.nelements()/2)], atermMaker_p->getTelescopeName());
2125 : // cerr << "BANDNAME " << bandname << " telescip " << atermMaker_p->getTelescopeName() << " csys tel " << csys.obsInfo().telescope() << endl;
2126 0 : std::tie(cell,convnx)=getBeamCellSize(bandname);
2127 : // cerr << "@@@cell " << cell << " npix " << convnx << endl;
2128 0 : csys_p=csys;
2129 0 : Vector<String> units=csys_p.worldAxisUnits();
2130 0 : Vector<Double> incr=csys_p.increment();
2131 0 : Double inpFov=fabs(incr[0]*npix);
2132 0 : incr[0]=cell.get(units[0]).getValue();
2133 0 : incr[1]=cell.get(units[1]).getValue();
2134 : //cerr << "###inp fov" << inpFov << " conv fov " << fabs(incr[0]*Double(convnx)) << endl;
2135 0 : Double pbFov= fabs(incr[0]*Double(convnx));
2136 : /*if(inpFov < fabs(incr[0]*Double(convnx))){
2137 : //incr = csys_p.increment();
2138 : npix=int(std::ceil(incr[0]*Double(convnx)/inpFov/2.0))*2;
2139 : //npix remains the same and csys used for beam calc stays
2140 : }
2141 : else{
2142 : npix = convnx; // return the npix used to calc beam
2143 : }*/
2144 0 : npix=convnx;
2145 0 : if((inpFov/pbFov) < 1.0){
2146 0 : npix=int(std::ceil(inpFov/pbFov*Double(convnx)/2.0))*2;
2147 : //cerr << "$$$$ npix " << npix << " cnx " << convnx << endl;
2148 : }
2149 0 : Vector<Int> stoks={Stokes::RR, Stokes::RL, Stokes::LR, Stokes::LL};
2150 0 : StokesCoordinate stokesCoords(stoks);
2151 0 : Quantum<Vector<Double> > freqs(freqlist, "Hz");
2152 0 : SpectralCoordinate specCoord(MFrequency::TOPO, freqs);
2153 0 : CoordinateSystem csysA=csys_p;
2154 0 : csysA.setIncrement(incr);
2155 0 : Vector<Double> refpix=csysA.referencePixel();
2156 0 : refpix[0]=refpix[1]=Double(convnx)/2.0;
2157 0 : csysA.setReferencePixel(refpix);
2158 0 : csysA.replaceCoordinate(stokesCoords, 1);
2159 0 : csysA.replaceCoordinate(specCoord,2);
2160 0 : csys=csysA; //return the coordinate used for doing the beam
2161 0 : IPosition shp(4, convnx, convnx, 4,1);
2162 0 : Int support=0;
2163 : //aa.cacheVBInfo("EVLA", 25.0);
2164 0 : IPosition blc(4,0,0,0,0);
2165 0 : IPosition trc(4,0,0,3,0);
2166 0 : FFT2D ftm;
2167 : Bool isCopy, isWtCopy;
2168 0 : uInt nchans=freqlist.nelements();
2169 0 : MathUtils m;
2170 0 : Record miscInfo;
2171 0 : miscInfo.define("bandname", bandname);
2172 : //cerr << "MISCINFO " << miscInfo << endl;
2173 0 : for (uInt k=0; k < nchans; ++k){
2174 : //higest freq will have largest supp ..so going through freqlist backwards
2175 0 : Quantum<Vector<Double> > lefreq(Vector<Double>(1, freqlist[nchans-k-1]), "Hz");
2176 0 : SpectralCoordinate specplane(MFrequency::TOPO, lefreq);
2177 0 : CoordinateSystem csysPlane=csysA;
2178 0 : csysPlane.replaceCoordinate(specplane,2);
2179 0 : TempImage<Complex> pbim(shp, csysPlane);
2180 0 : pbim.setMiscInfo(miscInfo);
2181 0 : pbim.set(1.0);
2182 0 : if(doSquint)
2183 0 : atermMaker_p->applyDiagSkyJones(pbim, pa);
2184 : else
2185 0 : atermMaker_p->applyAvgSkyJones(pbim);
2186 0 : Array<Complex> arr;
2187 0 : if(inpFov >= pbFov){
2188 0 : arr=pbim.getSlice(IPosition(4,0), IPosition(4, convnx, convnx, 4,1), False);
2189 : }
2190 : else{
2191 0 : arr=pbim.getSlice(IPosition(4,(convnx-npix)/2, (convnx-npix)/2, 0, 0), IPosition(4, npix, npix, 4,1), False);
2192 : }
2193 : //cerr << "MAX arr "<< max(arr) << " min "<< min(arr) << endl;
2194 0 : Array<Complex> wtArr=arr*conj(arr);
2195 0 : Complex * arrptr=arr.getStorage(isCopy);
2196 0 : Complex * wtarrptr=wtArr.getStorage(isWtCopy);
2197 0 : for(uint j=0; j <4; ++j){
2198 0 : Complex* arrplane= arrptr+j*npix*npix;
2199 0 : Complex* wtarrplane=wtarrptr+j*npix*npix;
2200 0 : ftm.c2cFFT(arrplane, npix, npix, True);
2201 0 : ftm.c2cFFT(wtarrplane, npix, npix, True);
2202 : }
2203 :
2204 0 : arr.putStorage(arrptr, isCopy);
2205 :
2206 0 : wtArr.putStorage(wtarrptr, isWtCopy);
2207 : /*if(inpFov < pbFov){
2208 : Double fac= (inpFov/pbFov);
2209 : incr[0]/=fac;
2210 : incr[1]/=fac;
2211 : csys.setIncrement(incr);
2212 : Array<Complex>newArr= m.resampleViaFFT(arr, fac, fac);
2213 : arr.set(0.0);
2214 : cerr << "@@@neArr shape "<< newArr.shape() << endl;
2215 : m.putMiddle(arr, newArr);
2216 : newArr.resize();
2217 : newArr=m.resampleViaFFT(wtArr, fac, fac);
2218 : wtArr.set(0.0);
2219 : m.putMiddle(wtArr, newArr);
2220 :
2221 :
2222 : }*/
2223 : //cerr << "Post FT MAX arr "<< max(wtArr) << " min "<< min(wtArr) << endl;
2224 0 : supportAndNormalizeAFunc(support, arr, wtArr);
2225 : //cerr << "Post Norm MAX arr "<< max(arr) << " min "<< min(arr) << endl;
2226 0 : if(k==0){
2227 :
2228 0 : IPosition tmpShape=shp;
2229 0 : tmpShape[3]=freqlist.nelements();
2230 : //This is going to be returned to allow resacling etc later
2231 : //csys=refim::SynthesisUtils::makeUVCoords(csysA, tmpShape);
2232 0 : IPosition convShp(4, 2*support, 2*support, 4, freqlist.nelements());
2233 0 : refpix[0]=refpix[1]=support;
2234 0 : csys.setReferencePixel(refpix);
2235 0 : convFunc.resize(convShp);
2236 : //convShp[0]=convShp[1]=2*support;
2237 0 : wtConvFunc.resize(convShp);
2238 0 : blc[0]=blc[1]=npix/2-support;
2239 0 : trc[0]=trc[1]=npix/2+support-1;
2240 :
2241 0 : }
2242 0 : asupport[nchans-k-1]=support;
2243 0 : IPosition blcChanplane=IPosition(4, 0,0,0,nchans-k-1);
2244 0 : IPosition trcChanplane=IPosition(4, convFunc.shape()[0]-1,convFunc.shape()[1]-1,3,nchans-k-1);
2245 :
2246 0 : convFunc(blcChanplane, trcChanplane)=arr(blc,trc);
2247 0 : wtConvFunc(blcChanplane, trcChanplane)=wtArr(blc, trc);
2248 :
2249 0 : }
2250 : //cerr << "Post Slicing MAX arr "<< max(wtConvFunc) << " min "<< min(wtConvFunc) << endl;
2251 : //cerr << "@@@support " << max(asupport) << endl;
2252 0 : }
2253 :
2254 0 : void AWConvFunc::makeAWConvFunc(Array<Complex>& convFunc,
2255 : Array<Complex>& wtconv,
2256 : CoordinateSystem& csys,
2257 : Matrix<Int>& awsupport, Int& npix, const Vector<Double>& freqlist,
2258 : const Vector<Double>& wVals, const Bool dosquint, const Double& pa){
2259 :
2260 :
2261 0 : Array<Complex> pbFT;
2262 0 : Array<Complex> pbWFT;
2263 0 : Vector<Int> aTsup;
2264 0 : makeAConvFunc(pbFT, pbWFT, csys, aTsup, npix, freqlist, dosquint, pa);
2265 : //cerr << "In AW wtConv "<< max(pbWFT) << " min "<< min(pbWFT) << endl;
2266 : /////TESTOO
2267 : //{
2268 : // cerr << "ATSUP " << aTsup << endl;
2269 : //PagedImage<Complex> booltaq(pbFT.shape(), csys, "BOOLTAQ");
2270 : //booltaq.put(pbFT);
2271 : //}
2272 : /////
2273 0 : WPConvFunc wpc;
2274 0 : Cube<Complex> wCon;
2275 0 : Vector<int> wTsup;
2276 0 : Int nfreq=freqlist.nelements();
2277 0 : Int nw=wVals.nelements();
2278 0 : wpc.makeWConvFuncs(wCon, wTsup, csys, npix, wVals);
2279 : //cerr << "AWC wcon shape "<< wCon.shape() << " pbFT " << pbFT.shape() << endl;
2280 0 : Int newNx=max(wCon.shape()[0], pbFT.shape()[0]);
2281 : // Let's start with this size..we can reduce this later
2282 0 : convFunc.resize(IPosition(5, newNx, newNx, 4, nfreq, nw));
2283 0 : IPosition blcW(2, (newNx-wCon.shape()[0])/2, (newNx-wCon.shape()[1])/2);
2284 0 : IPosition trcW(2, (newNx+wCon.shape()[0])/2-1, (newNx+wCon.shape()[1])/2-1);
2285 0 : IPosition blcA(4, (newNx-pbFT.shape()[0])/2, (newNx-pbFT.shape()[1])/2, 0, 0);
2286 0 : IPosition trcA(4, (newNx+pbFT.shape()[0])/2-1, (newNx+pbFT.shape()[1])/2-1, 3, nfreq-1);
2287 0 : FFT2D ftm;
2288 : //IPosition blcOut(5, 0, 0, 0, 0, 0);
2289 : //IPosition trcOut(5, newNx-1, newNx-1, 4, nfreq-1, 0);
2290 : ArrayIterator<Complex> itOut(convFunc,
2291 0 : IPosition(4,0,1,2,3));
2292 0 : itOut.origin();
2293 0 : for (Int k=0; k < nw; ++k){
2294 0 : Matrix<Complex> w(newNx, newNx,0);
2295 0 : w(blcW, trcW)=wCon.xyPlane(k);
2296 : Bool isCopy;
2297 0 : Complex* wptr=w.getStorage(isCopy);
2298 0 : ftm.c2cFFT(wptr, newNx, newNx, False);
2299 0 : w.putStorage(wptr, isCopy);
2300 0 : Array<Complex> pb(IPosition(4, newNx, newNx, 4, nfreq));
2301 0 : pb.set(0.0);
2302 0 : pb(blcA, trcA)=pbFT;
2303 0 : ArrayIterator<Complex> it(pb, IPosition(3,0,1,2));
2304 0 : it.origin();
2305 0 : for (Int j = 0; j < nfreq; ++j){
2306 : //IPosition blcf(4, 0, 0, 0, j);
2307 : //IPosition trcf(4, newNx-1, newNx-1, 3,j);
2308 0 : Cube<Complex> aPB(it.array());
2309 0 : for (Int pol=0; pol<4; ++pol){
2310 0 : Matrix<Complex> aPBplane=aPB.xyPlane(pol);
2311 0 : Complex* pbptr=aPBplane.getStorage(isCopy);
2312 0 : ftm.c2cFFT(pbptr, newNx, newNx, False);
2313 0 : aPBplane.putStorage(pbptr, isCopy);
2314 0 : aPBplane *= w;
2315 0 : pbptr=aPBplane.getStorage(isCopy);
2316 0 : ftm.c2cFFT(pbptr, newNx, newNx, True);
2317 0 : aPBplane.putStorage(pbptr, isCopy);
2318 :
2319 0 : }
2320 0 : it.next();
2321 0 : }
2322 : //blcOut[4]=k;
2323 : //trcOut[4]=k;
2324 : //convFunc(blcOut, trcOut)=pb;
2325 0 : itOut.array()=pb;
2326 : /////TESTOO
2327 : //{
2328 : //PagedImage<Complex> booltaq(pb.shape(), csys, "BOOLTAQ2");
2329 : //booltaq.put(pb);
2330 : //}
2331 : /////
2332 0 : itOut.next();
2333 0 : }
2334 : //cerr << "BEFORE convshape "<< convFunc.shape() << endl;
2335 0 : supportResizeAWConv(awsupport, convFunc, aTsup);
2336 : //For now we will copy the weightConvFunc to all W's just to keep indexing the same
2337 0 : wtconv.resize(convFunc.shape());
2338 0 : wtconv.set(0.0);
2339 0 : ArrayIterator<Complex>itw(wtconv, IPosition(4,0,1,2,3));
2340 0 : newNx=wtconv.shape()[0];
2341 0 : Int nx=pbWFT.shape()[0];
2342 0 : itw.origin();
2343 0 : IPosition blcw(4, (newNx-nx)/2, (newNx-nx)/2, 0, 0);
2344 0 : IPosition trcw(4, (newNx+nx)/2-1, (newNx+nx)/2-1, 3, nfreq-1);
2345 0 : for (int k=0; k < nw; ++k){
2346 : //cerr << "@@@ w " << k << " shape " << itw.array()(blcw,trcw).shape() << " " << pbWFT.shape() << endl;
2347 0 : itw.array()(blcw,trcw)=pbWFT;
2348 0 : itw.next();
2349 : }
2350 : //cerr << "AFTER convshape "<< convFunc.shape() << endl;
2351 : //cerr << "wtconv " << max(wtconv) << " " << min(wtconv) << endl;
2352 : ////TESTOO
2353 : /*{
2354 : Vector<Double>pixW(wVals.nelements());
2355 : indgen(pixW);
2356 : CoordinateSystem fiveAxis=csys;
2357 : TabularCoordinate tab(pixW, wVals, "m", "W");
2358 : fiveAxis.addCoordinate(tab);
2359 : PagedImage<Complex> noo(convFunc.shape(), fiveAxis, "AWConvVals");
2360 : noo.put(convFunc);
2361 :
2362 : }*/
2363 : ////
2364 : //cerr << "aWSup " << awsupport << endl;
2365 : //cerr << "aTsup " << aTsup << endl;
2366 :
2367 0 : }
2368 : //
2369 :
2370 0 : Bool AWConvFunc::supportResizeAWConv(Matrix<Int>& sup, Array<Complex>& conv, const Vector<Int>& ATsup){
2371 :
2372 0 : Int convSize=conv.shape()[0];
2373 0 : sup.resize(conv.shape()[3], conv.shape()[4]);
2374 0 : sup.set(0);
2375 : Float maxAbsConvFunc, minAbsConvFunc;
2376 0 : IPosition minpos, maxpos;
2377 0 : ArrayIterator<Complex> it(conv, IPosition(2,0,1));
2378 : //The w=0 will determine the sum under conv func for normalization
2379 0 : Vector<Double> sumUnder(conv.shape()[3], 0.0);
2380 0 : it.origin();
2381 0 : for (Int w=0; w < conv.shape()[4]; ++w){
2382 0 : for (Int chan=0; chan < conv.shape()[3]; ++chan){
2383 0 : Matrix<Complex> convPlane(it.array());
2384 0 : minMax(minAbsConvFunc, maxAbsConvFunc, minpos, maxpos, amplitude(convPlane));
2385 : //cerr << "minMAX " << minAbsConvFunc << " " << maxAbsConvFunc << endl;
2386 0 : Bool found=false;
2387 0 : Int trial=0;
2388 0 : for (trial=convSize/2-2; trial>0; --trial) {
2389 : //Searching down a diagonal
2390 0 : if(abs(convPlane(convSize/2-trial,convSize/2-trial)) > (1.0e-2*maxAbsConvFunc) ) {
2391 0 : found=true;
2392 0 : trial=Int(sqrt(2.0*Float(trial*trial)));
2393 :
2394 0 : break;
2395 : }
2396 : }
2397 0 : if(trial < 5)
2398 0 : trial=5;
2399 0 : if(found){
2400 0 : sup(chan, w)=trial;
2401 : }
2402 : else{
2403 0 : sup(chan,w)=5;
2404 : }
2405 0 : if(sup(chan, w) < ATsup(chan)){
2406 : //cerr << chan << " w " << w << "sup " << sup(chan,w) << " ATsup " << ATsup(chan) << endl;
2407 0 : sup(chan,w)=ATsup(chan);
2408 : }
2409 0 : if(w==0){
2410 0 : IPosition blc(2,-sup(chan,0)+convSize/2, -sup(chan,0)+convSize/2);
2411 0 : IPosition trc(2, sup(chan,0)+convSize/2-1, sup(chan, 0)+convSize/2-1);
2412 : //cerr << "chan " << chan << " blc " << blc << " trc " << trc << " sup "<< sup(chan,0) << endl;
2413 0 : sumUnder(chan)=real(sum(convPlane(blc,trc)));
2414 :
2415 0 : }
2416 0 : for (Int pol=0; pol < 4; ++pol) //skip to the next pol RR plane of next chan, w
2417 0 : it.next();
2418 : //cerr << "POSITION " << it.pos() << endl;
2419 0 : }
2420 : }
2421 : //Some thing to expt here as we are limiting w's in a beam or slightly moer
2422 0 : Int newConvSize = 2*max(sup);
2423 : //Int newConvSize = 2*mean(sup);
2424 0 : if(newConvSize < convSize){
2425 0 : Array<Complex> newConv(IPosition(4, newConvSize, newConvSize, 4, conv.shape()[3], conv.shape()[4]));
2426 0 : IPosition blc(5,(-newConvSize+convSize)/2, (-newConvSize+convSize)/2, 0, 0, 0);
2427 0 : IPosition trc(2, (newConvSize+convSize)/2-1, (newConvSize+convSize)/2-1, 3, conv.shape()[3]-1, conv.shape()[4]-1);
2428 0 : newConv=conv(blc,trc);
2429 0 : conv.resize();
2430 0 : conv.assign(newConv);
2431 0 : }
2432 : else
2433 0 : newConvSize=convSize;
2434 0 : ArrayIterator<Complex> it2(conv, IPosition(3, 0, 1, 2));
2435 0 : it2.origin();
2436 0 : for (Int w=0; w< conv.shape()[4]; ++w){
2437 0 : for (Int chan=0; chan < conv.shape()[3]; ++chan){
2438 0 : if(sumUnder(chan) >0.0){
2439 0 : Complex mult(1.0/sumUnder(chan));
2440 0 : it2.array() *= mult;
2441 : //cerr << std::setprecision(12) << "it2.shape " << it2.array().shape() << " //pos " << it2.pos() << "chan " << chan << " sumUnder " << sumUnder(chan) << " w " << w << endl;
2442 : }
2443 0 : it2.next();
2444 : }
2445 :
2446 : }
2447 0 : return True;
2448 :
2449 0 : }
2450 0 : Bool AWConvFunc::supportAndNormalizeAFunc(Int& sup, Array<Complex>& conv, Array<Complex>& wtconv){
2451 0 : sup=-1;
2452 0 : IPosition begin(4, 0, 0, 0, 0);
2453 0 : IPosition end=conv.shape()-1;
2454 0 : Int convSize=conv.shape()[0];
2455 0 : end[2]=0;
2456 0 : Matrix<Complex> convPlane(wtconv(begin,end).reform(IPosition(2, convSize, convSize)));
2457 : Float maxAbsConvFunc, minAbsConvFunc;
2458 0 : IPosition minpos, maxpos;
2459 0 : minMax(minAbsConvFunc, maxAbsConvFunc, minpos, maxpos, amplitude(convPlane));
2460 0 : Bool found=false;
2461 0 : Int trial=0;
2462 0 : for (trial=convSize/2-2; trial>0; trial--) {
2463 : //Searching down a diagonal
2464 0 : if(abs(convPlane(convSize/2-trial,convSize/2-trial)) > (1.0e-2*maxAbsConvFunc) ) {
2465 0 : found=true;
2466 0 : trial=Int(sqrt(2.0*Float(trial*trial)));
2467 :
2468 0 : break;
2469 : }
2470 : }
2471 0 : if(!found) {
2472 0 : if((maxAbsConvFunc-minAbsConvFunc) > (1.0e-2*maxAbsConvFunc))
2473 0 : found=true;
2474 : // if it drops by more than 2 magnitudes per pixel
2475 0 : trial= (convSize >10) ? 5 : (convSize/2 - 4);
2476 : }
2477 :
2478 :
2479 0 : if(found) {
2480 0 : if(trial < 5)
2481 0 : trial= ( (convSize >10 ) ? 5 : (convSize/2 - 4));
2482 0 : sup=trial+1;
2483 : //support is really over the edge
2484 0 : if( (sup >= convSize/2)) {
2485 0 : sup=convSize/2-1;
2486 : }
2487 :
2488 0 : ArrayIterator<Complex> pbIt(conv, IPosition(2, 0,1));
2489 0 : ArrayIterator<Complex> wtIt(wtconv, IPosition(2, 0,1));
2490 0 : IPosition blc(2,-sup+convSize/2, -sup+convSize/2);
2491 0 : IPosition trc(2, sup+convSize/2, sup+convSize/2);
2492 0 : Double pbSum=0.0;
2493 : //Iterate of pol
2494 0 : while(!pbIt.pastEnd()){
2495 0 : Matrix<Complex> pbplane(pbIt.array());
2496 0 : Matrix<Complex> wtplane(wtIt.array());
2497 0 : pbSum=real(sum(wtplane(blc, trc)));
2498 : //cerr << "pbSumWt "<< pbSum << endl;
2499 0 : if(pbSum > 0.0){
2500 0 : wtplane=wtplane*Complex(1.0/pbSum, 0.0);
2501 0 : pbSum=real(sum(pbplane(blc,trc)));
2502 : //cerr << "pbSum "<< pbSum << endl;
2503 0 : pbplane=pbplane*Complex(1.0/pbSum,0.0);
2504 : }
2505 0 : pbIt.next();
2506 0 : wtIt.next();
2507 :
2508 0 : }
2509 0 : }
2510 :
2511 0 : return found;
2512 0 : }
2513 0 : std::pair<Quantity, int> AWConvFunc::getBeamCellSize(const String& band){
2514 :
2515 0 : Quantity fov(0.024,"rad"); //fov at 1 GHz for VLA
2516 0 : Quantity cell=fov/256;
2517 0 : if(band=="EVLA_S")
2518 0 : cell=cell/2.0;
2519 0 : else if(band=="EVLA_C" || band=="VLA_C")
2520 0 : cell=cell/4.0;
2521 0 : else if(band=="EVLA_X" || band=="VLA_X")
2522 0 : cell=cell/8.0;
2523 0 : else if(band=="EVLA_U" || band=="VLA_U")
2524 0 : cell=cell/12.0;
2525 0 : else if(band=="EVLA_K" || band == "VLA_K")
2526 0 : cell = cell/18.0;
2527 0 : else if(band == "EVLA_A")
2528 0 : cell = cell/26.0;
2529 0 : else if(band == "EVLA_Q" || band == "VLA_Q")
2530 0 : cell=cell/40.0;
2531 :
2532 : //we can optimize this later but we'll use 512 pixels for nRow
2533 : //ahem....BeamCalc starts deviating if we want to make VP and PB
2534 : //that are highly oversampled
2535 :
2536 0 : return std::pair<Quantity, int>(cell, 512);
2537 :
2538 0 : }
2539 :
2540 : //
2541 : //----------------------------------------------------------------------
2542 : //
2543 : // Vector<Vector<Double> > AWConvFunc::findPointingOffset(const ImageInterface<Complex>& image,
2544 : // const VisBuffer2& vb, const Bool& doPointing)
2545 : // {
2546 : // Assert(po_p.null()==False && "Pointingoffset call has not been initialized in AWProjectFT call being made");
2547 : // return po_p->findPointingOffset(image,vb,doPointing);
2548 : // // if (!doPointing)
2549 : // // {cerr<<"AWCF: Using mosaic pointing \n";return po_p->findMosaicPointingOffset(image,vb);}
2550 : // // else
2551 : // // {cerr<<"AWCF: Using antenna pointing table \n";return po_p->findAntennaPointingOffset(image,vb);}
2552 : // }
2553 :
2554 :
2555 :
2556 : };
2557 : };
|
