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