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