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