Line data Source code
1 : // -*- C++ -*-
2 : //# AWVisResampler.cc: Implementation of the AWVisResampler 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/SynthesisError.h>
30 : #include <synthesis/TransformMachines2/AWVisResampler.h>
31 : #include <synthesis/TransformMachines2/Utils.h>
32 : #include <synthesis/TransformMachines/SynthesisMath.h>
33 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
34 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
35 : #include <casacore/casa/OS/Timer.h>
36 : #include <fstream>
37 : #include <iostream>
38 : #include <typeinfo>
39 : #include <iomanip>
40 : #include <cfenv>
41 : //#include <synthesis/TransformMachines2/FortranizedLoops.h>
42 : //#include <synthesis/TransformMachines2/hpg.hpp>
43 :
44 : #ifdef _OPENMP
45 : #include <omp.h>
46 : #endif
47 : extern "C" {
48 : void clLoopsToGrid();
49 : };
50 : //#include <casa/BasicMath/Functors.h>
51 : using namespace casacore;
52 : //using namespace hpg;
53 : namespace casa{
54 : using namespace refim;
55 : //
56 : //-----------------------------------------------------------------------------------
57 : // Re-sample the griddedData on the VisBuffer (a.k.a gridding)
58 : //
59 : // Template instantiations for re-sampling onto a double precision
60 : // or single precision grid.
61 : //
62 : template
63 : void AWVisResampler::addTo4DArray(DComplex* __restrict__ & store,
64 : const Int* __restrict__ & iPos,
65 : const Vector<Int>& inc,
66 : Complex& nvalue, Complex& wt) __restrict__ ;
67 : template
68 : void AWVisResampler::addTo4DArray(Complex* __restrict__ & store,
69 : const Int* __restrict__ & iPos,
70 : const Vector<Int>& inc,
71 : Complex& nvalue, Complex& wt) __restrict__;
72 :
73 : template
74 : void AWVisResampler::addTo4DArray_ptr(DComplex* __restrict__ & store,
75 : const Int* __restrict__ & iPos,
76 : const Int* __restrict__ & inc,
77 : Complex& nvalue, Complex& wt) __restrict__ ;
78 : template
79 : void AWVisResampler::addTo4DArray_ptr(Complex* __restrict__ & store,
80 : const Int* __restrict__ & iPos,
81 : const Int* __restrict__ & inc,
82 : Complex& nvalue, Complex& wt) __restrict__ ;
83 :
84 : template
85 : void AWVisResampler::DataToGridImpl_p(Array<DComplex>& grid, VBStore& vbs,
86 : Matrix<Double>& sumwt,const Bool& dopsf,
87 : Bool useConjFreqCF); // __restrict__;
88 : template
89 : void AWVisResampler::DataToGridImpl_p(Array<Complex>& grid, VBStore& vbs,
90 : Matrix<Double>& sumwt,const Bool& dopsf,
91 : Bool useConjFreqCF); // __restrict__;
92 :
93 7534608 : Complex* AWVisResampler::getConvFunc_p(const double& vbPA, Vector<Int>& cfShape,
94 : Vector<int>& support,
95 : int& muellerElement,
96 : CountedPtr<CFBuffer>& cfb,
97 : Double& wVal, Int& fndx, Int& wndx,
98 : PolMapType& mNdx, PolMapType& conjMNdx,
99 : Int& ipol, uInt& mRow)
100 : {
101 : Bool Dummy;
102 : Array<Complex> *convFuncV;
103 : CFCell *cfcell;
104 : //
105 : // Since we conjugate the CF depending on the sign of the w-value,
106 : // pick the appropriate element of the Mueller Matrix (see note on
107 : // this for details). Without reading the note/understanding,
108 : // fiddle with this logic at your own risk (can easily lead to a
109 : // lot of grief. --Sanjay).
110 : //
111 : //timer_p.mark();
112 :
113 : int pndx;
114 7534608 : if (wVal > 0.0)
115 : {
116 2259366 : pndx=mNdx[ipol][mRow];
117 : // cfcell=&(*(cfb->getCFCellPtr(fndx,wndx,mNdx[ipol][mRow])));
118 : // CFCell& cfO=cfb(fndx,wndx,mNdx[ipol][mRow]);
119 : // convFuncV = &(*cfO.getStorage());
120 : // support(0)=support(1)=cfO.xSupport_p;
121 : }
122 : else
123 : {
124 5275242 : pndx=conjMNdx[ipol][mRow];
125 : // cfcell=&(*(cfb->getCFCellPtr(fndx,wndx,conjMNdx[ipol][mRow])));
126 : // CFCell& cfO=cfb(fndx,wndx,conjMNdx[ipol][mRow]);
127 : // convFuncV = &(*cfO.getStorage());
128 : // support(0)=support(1)=cfO.xSupport_p;
129 : }
130 : //pndx=1;
131 : // cerr << "CFC indexes(w,f,p),ipol: " << wndx << " " << fndx << " " << pndx << " " << ipol << endl;
132 7534608 : cfcell=&(*(cfb->getCFCellPtr(fndx,wndx,pndx)));
133 :
134 : //cerr << "CF Name: " << cfcell->fileName_p << endl;
135 :
136 7534608 : convFuncV = &(*cfcell->getStorage());
137 7534608 : support(0)=support(1)=cfcell->xSupport_p;
138 :
139 : // Get the pointer to the CFCell storage (a single CF)
140 : // if ((convFuncV = &(*cfcell->getStorage())) == NULL)
141 7534608 : if (convFuncV == NULL)
142 0 : throw(SynthesisFTMachineError("cfcell->getStorage() == null"));
143 :
144 : // Load the CF if it not already loaded. If a new CF is loaded,
145 : // check if it needs to be rotated.
146 7534608 : if (convFuncV->shape().product() == 0)
147 : {
148 0 : Array<Complex> tt=SynthesisUtils::getCFPixels(cfb->getCFCacheDir(), cfcell->fileName_p);
149 0 : cfcell->setStorage(tt);
150 :
151 : //cerr << (cfcell->isRotationallySymmetric_p?"o":"+");
152 :
153 : // No rotation necessary if the CF is rotationally symmetric
154 0 : if (!(cfcell->isRotationallySymmetric_p))
155 : {
156 0 : CFCell *baseCFC=NULL;
157 : // Rotate only if the difference between CF PA and VB PA
158 : // is greater than paTolerance.
159 0 : SynthesisUtils::rotate2(vbPA, *baseCFC, *cfcell, paTolerance_p);
160 : }
161 0 : convFuncV = &(*cfcell->getStorage());
162 0 : }
163 :
164 : //cfShape.reference(cfcell->cfShape_p);
165 7534608 : cfShape.assign(convFuncV->shape().asVector());
166 :
167 : // Always extract the Mueller element value from mNdx. mNdx
168 : // carries the direct mapping between Mueller Matrix and
169 : // Visibility vector.
170 : // muellerElement=cfb->getCFCellPtr(fndx,wndx,mNdx[ipol][mRow])->muellerElement_p;
171 7534608 : muellerElement=cfcell->muellerElement_p;
172 :
173 : // cfShape.assign(cfcell->cfShape_p);
174 : //runTimeG1_p += timer_p.real();
175 :
176 :
177 15069216 : return convFuncV->getStorage(Dummy);
178 : };
179 :
180 : //
181 : //-----------------------------------------------------------------------------------
182 : //
183 0 : void AWVisResampler::cachePhaseGrad_p(const Vector<Double>& pointingOffset,
184 : const Vector<Int>&cfShape,
185 : const Vector<Int>& convOrigin,
186 : const Double& /*cfRefFreq*/,
187 : const Double& /*imRefFreq*/,
188 : const Int& spwID, const Int& fieldId)
189 : {
190 0 : if (
191 0 : ((fabs(pointingOffset[0]-cached_PointingOffset_p[0])) > 1e-6) ||
192 0 : ((fabs(pointingOffset[1]-cached_PointingOffset_p[1])) > 1e-6) ||
193 0 : (cached_phaseGrad_p.shape()[0] < cfShape[0]) ||
194 0 : (cached_phaseGrad_p.shape()[1] < cfShape[1])
195 : )
196 : {
197 0 : LogIO log_l(LogOrigin("AWVisResampler","cachePhaseGrad[R&D]"));
198 : log_l << "Computing phase gradiant for pointing offset "
199 : << pointingOffset << cfShape << " " << cached_phaseGrad_p.shape()
200 : << "(SPW: " << spwID << " Field: " << fieldId << ")"
201 : << LogIO::DEBUGGING
202 0 : << LogIO::POST;
203 0 : Int nx=cfShape(0), ny=cfShape(1);
204 : Double grad;
205 0 : Complex phx,phy;
206 :
207 0 : cached_phaseGrad_p.resize(nx,ny);
208 0 : cached_PointingOffset_p = pointingOffset;
209 :
210 0 : for(Int ix=0;ix<nx;ix++)
211 : {
212 0 : grad = (ix-convOrigin[0])*pointingOffset[0];
213 : Double sx,cx;
214 0 : SINCOS(grad,sx,cx);
215 :
216 0 : phx = Complex(cx,sx);
217 0 : for(Int iy=0;iy<ny;iy++)
218 : {
219 0 : grad = (iy-convOrigin[1])*pointingOffset[1];
220 : Double sy,cy;
221 0 : SINCOS(grad,sy,cy);
222 0 : phy = Complex(cy,sy);
223 0 : cached_phaseGrad_p(ix,iy)=phx*phy;
224 : }
225 : }
226 0 : }
227 0 : }
228 :
229 : //
230 : //-----------------------------------------------------------------------------------
231 : // Template implementation for DataToGrid
232 : //
233 : template <class T>
234 28083 : void AWVisResampler::DataToGridImpl_p(Array<T>& grid, VBStore& vbs,
235 : Matrix<Double>& sumwt,const Bool& dopsf,
236 : Bool /*useConjFreqCF*/)
237 : {
238 : Int nDataChan, nDataPol, nGridPol, nGridChan, nx, ny, nw;//, nCFFreq;
239 : Int targetIMChan, targetIMPol, rbeg, rend;//, PolnPlane, ConjPlane;
240 : Int startChan, endChan;
241 :
242 :
243 28083 : Vector<Float> sampling(2),scaledSampling(2);
244 28083 : Vector<Int> support(2),loc(3), iloc(4),scaledSupport(2);
245 28083 : Vector<Double> pos(3), off(3);
246 28083 : Vector<Int> igrdpos(4);
247 :
248 28083 : Complex phasor, nvalue /*, wt */;
249 28083 : DComplex norm;
250 28083 : Vector<Int> cfShape;
251 : // cfShape=(*vb2CFBMap_p)[0]->getStorage()(0,0,0)->getStorage()->shape().asVector();
252 :
253 28083 : Vector<Int> convOrigin;// = (cfShape)/2;
254 : Double cfRefFreq;
255 28083 : const Matrix<Double> UVW=vbs.uvw_p;
256 :
257 : // Double cfScale=1.0;
258 :
259 : //timer_p.mark();
260 28083 : rbeg = 0; rend = vbs.nRow_p;
261 28083 : rbeg = vbs.beginRow_p;
262 28083 : rend = vbs.endRow_p;
263 :
264 28083 : nx = grid.shape()[0]; ny = grid.shape()[1];
265 28083 : nGridPol = grid.shape()[2]; nGridChan = grid.shape()[3];
266 :
267 28083 : nDataPol = vbs.flagCube_p.shape()[0];
268 28083 : nDataChan = vbs.flagCube_p.shape()[1];
269 :
270 : Bool Dummy, gDummy,
271 28083 : accumCFs=((UVW.nelements() == 0) && dopsf);
272 :
273 28083 : T* __restrict__ gridStore = grid.getStorage(gDummy);
274 :
275 28083 : Double *freq=vbs.freq_p.getStorage(Dummy);
276 :
277 28083 : cacheAxisIncrements(grid.shape().asVector(), gridInc_p);
278 :
279 28083 : Bool * __restrict__ flagCube_ptr=vbs.flagCube_p.getStorage(Dummy);
280 28083 : Bool * __restrict__ rowFlag_ptr = vbs.rowFlag_p.getStorage(Dummy);;
281 28083 : Float * __restrict__ imgWts_ptr = vbs.imagingWeight_p.getStorage(Dummy);
282 28083 : Complex * __restrict__ visCube_ptr = vbs.visCube_p.getStorage(Dummy);
283 :
284 28083 : Vector<Double> wVals, fVals; PolMapType mVals, mNdx, conjMVals, conjMNdx;
285 : Double fIncr, wIncr;
286 28083 : CountedPtr<CFBuffer> cfb = (*vb2CFBMap_p)[0];
287 28083 : bool finitePointingOffsets=cfb->finitePointingOffsets();
288 :
289 28083 : cfb->getCoordList(fVals,wVals,mNdx, mVals, conjMNdx, conjMVals, fIncr, wIncr);
290 :
291 28083 : nw = wVals.nelements();
292 28083 : iloc = 0;
293 :
294 28083 : IPosition shp=vbs.flagCube_p.shape();
295 28083 : if (accumCFs)
296 : {
297 7150 : startChan = vbs.startChan_p;
298 7150 : endChan = vbs.endChan_p;
299 : }
300 : else
301 : {
302 20933 : startChan = 0;
303 20933 : endChan = nDataChan;
304 : }
305 :
306 28083 : Int vbSpw = (vbs.vb_p)->spectralWindows()(0);
307 28083 : Double vbPA = vbs.paQuant_p.getValue("rad");
308 :
309 3301599 : for(Int irow=rbeg; irow< rend; irow++){
310 :
311 3273516 : if(!(*(rowFlag_ptr+irow)))
312 : {
313 3269652 : setFieldPhaseGrad((vb2CFBMap_p->vectorPhaseGradCalculator_p[vb2CFBMap_p->vbRow2BLMap_p[irow]])->field_phaseGrad_p);
314 3269652 : cfb = (*vb2CFBMap_p)[irow];
315 :
316 6539304 : for(Int ichan=startChan; ichan< endChan; ichan++)
317 : {
318 3269652 : if (*(imgWts_ptr + ichan+irow*nDataChan)!=0.0)
319 : {
320 3269652 : targetIMChan=chanMap_p[ichan];
321 :
322 3269652 : if((targetIMChan>=0) && (targetIMChan<nGridChan))
323 : {
324 :
325 :
326 3269652 : Double dataWVal = (UVW.nelements() > 0) ? UVW(2,irow) : 0.0;
327 3269652 : Int wndx = cfb->nearestWNdx(dataWVal*freq[ichan]/C::c);
328 3269652 : Int cfFreqNdx = cfb->nearestFreqNdx(vbSpw,ichan,vbs.conjBeams_p);
329 : Float s;
330 : //
331 : //------------------------------------------------------------------------------
332 : //
333 : // Using the int-index version for Freq, W and Muellerelements
334 : //
335 : //------------------------------------------------------------------------------
336 : //
337 :
338 3269652 : cfb->getParams(cfRefFreq, s, support(0), support(1),cfFreqNdx,wndx,0);
339 :
340 3269652 : sampling(0) = sampling(1) = SynthesisUtils::nint(s);
341 :
342 3269652 : sgrid(pos,loc,off, phasor, irow, UVW, dphase_p[irow], freq[ichan],
343 3269652 : uvwScale_p, offset_p, sampling);
344 : {
345 : // Loop over all image-plane polarization planes.
346 :
347 9808956 : for(Int ipol=0; ipol< nDataPol; ipol++)
348 : {
349 6539304 : if((!(*(flagCube_ptr + ipol + ichan*nDataPol + irow*nDataPol*nDataChan))))
350 : {
351 6539304 : targetIMPol=polMap_p(ipol);
352 6539304 : if ((targetIMPol>=0) && (targetIMPol<nGridPol))
353 : {
354 6539304 : igrdpos[2]=targetIMPol; igrdpos[3]=targetIMChan;
355 :
356 6539304 : norm = 0.0;
357 : // Loop over all relevant elements of the Mueller matrix for the polarization
358 : // ipol.
359 6539304 : Vector<int> conjMRow = conjMNdx[ipol];
360 :
361 13040672 : for (uInt mCols=0;mCols<conjMRow.nelements(); mCols++)
362 : {
363 6539304 : int visVecElement=mCols, muellerElement;
364 :
365 6539304 : Complex* convFuncV=NULL;
366 : //timer_p.mark();
367 : try
368 : {
369 6539304 : convFuncV=getConvFunc_p(vbPA,
370 : cfShape, support,muellerElement,
371 : cfb, dataWVal, cfFreqNdx,
372 : wndx, mNdx, conjMNdx, ipol, mCols);
373 : }
374 0 : catch (SynthesisFTMachineError& x)
375 : {
376 0 : LogIO log_l(LogOrigin("AWVisResampler[R&D]","DataToGridImpl_p"));
377 0 : log_l << x.getMesg() << LogIO::EXCEPTION;
378 0 : }
379 : // Extract the vis. vector element corresponding to the mCols column of the conjMRow row of the Mueller matrix.
380 :
381 6539304 : visVecElement=(int)(muellerElement%nDataPol);
382 : // If the vis. vector element is flagged, don't grid it.
383 6577240 : if(((*(flagCube_ptr + visVecElement + ichan*nDataPol + irow*nDataPol*nDataChan)))) break;
384 :
385 6539304 : if(dopsf) nvalue=Complex(*(imgWts_ptr + ichan + irow*nDataChan));
386 2755092 : else nvalue=Complex(*(imgWts_ptr+ichan+irow*nDataChan))*
387 5510184 : (*(visCube_ptr+visVecElement+ichan*nDataPol+irow*nDataChan*nDataPol)*phasor);
388 :
389 6539304 : if (!onGrid(nx, ny, nw, loc, support)) break;
390 :
391 :
392 6501368 : if(cfShape[0]%2==0 && cfShape[1]%2==0)
393 1155072 : convOrigin=cfShape/2;
394 : else
395 5346296 : convOrigin=cfShape/2+1;
396 6501368 : cacheAxisIncrements(cfShape, cfInc_p);
397 6501368 : nVisGridded_p++;
398 : #include <synthesis/TransformMachines2/accumulateToGrid.inc>
399 : }
400 6539304 : sumwt(targetIMPol,targetIMChan) += vbs.imagingWeight_p(ichan, irow)*fabs(norm);
401 : // *(sumWt_ptr+apol+achan*nGridChan)+= *(imgWts_ptr+ichan+irow*nDataChan);
402 6539304 : }
403 : }
404 : } // End poln-loop
405 : }
406 : }
407 : }
408 : } // End chan-loop
409 : }
410 : } // End row-loop
411 :
412 28083 : T *tt=(T *)gridStore;
413 28083 : grid.putStorage(tt,gDummy);
414 28083 : }
415 : //
416 : //-----------------------------------------------------------------------------------
417 : // Re-sample VisBuffer to a regular grid (griddedData) (a.k.a. de-gridding)
418 : //
419 6102 : void AWVisResampler::GridToData(VBStore& vbs, const Array<Complex>& grid)
420 : {
421 : Int nDataChan, nDataPol, nGridPol, nGridChan, nx, ny,nw;//, nCFFreq;
422 : Int achan, apol, rbeg, rend;//, PolnPlane, ConjPlane;
423 6102 : Vector<Float> sampling(2);//scaledSampling(2);
424 6102 : Vector<Int> support(2),loc(3), iloc(4);
425 6102 : Vector<Double> pos(3), off(3);
426 :
427 6102 : IPosition grdpos(4);
428 :
429 6102 : Vector<Complex> norm(4,0.0);
430 6102 : Complex phasor, nvalue;
431 6102 : CountedPtr<CFBuffer> cfb=(*vb2CFBMap_p)[0];
432 6102 : Vector<Int> cfShape=cfb->getStorage()(0,0,0)->getStorage()->shape().asVector();
433 6102 : Bool finitePointingOffset=cfb->finitePointingOffsets();
434 :
435 :
436 6102 : Vector<Int> convOrigin; // = ((cfShape[0]%2==0) && (cfShape[1]%2==0)) ? (cfShape)/2 : cfShape/2+1;
437 : Double cfRefFreq;//cfScale=1.0
438 :
439 6102 : rbeg=0;
440 6102 : rend=vbs.nRow_p;
441 6102 : rbeg = vbs.beginRow_p;
442 6102 : rend = vbs.endRow_p;
443 6102 : nx = grid.shape()[0]; ny = grid.shape()[1];
444 6102 : nGridPol = grid.shape()[2]; nGridChan = grid.shape()[3];
445 :
446 6102 : nDataPol = vbs.flagCube_p.shape()[0];
447 6102 : nDataChan = vbs.flagCube_p.shape()[1];
448 :
449 : //
450 : // The following code reduces most array accesses to the simplest
451 : // possible to improve performance. However this made no
452 : // difference in the run-time performance compared to Vector,
453 : // Matrix and Cube indexing.
454 : //
455 : Bool Dummy;
456 6102 : const Complex* __restrict__ gridStore = grid.getStorage(Dummy);
457 : (void)gridStore;
458 6102 : Vector<Int> igrdpos(4);
459 6102 : Double *freq=vbs.freq_p.getStorage(Dummy);
460 6102 : Bool *rowFlag=vbs.rowFlag_p.getStorage(Dummy);
461 :
462 6102 : Matrix<Double>& uvw=vbs.uvw_p;
463 6102 : Cube<Complex>& visCube=vbs.visCube_p;
464 6102 : Cube<Bool>& flagCube=vbs.flagCube_p;
465 :
466 6102 : Vector<Int> gridInc, cfInc;
467 :
468 6102 : cacheAxisIncrements(grid.shape().asVector(), gridInc_p);
469 : //cacheAxisIncrements(cfShape, cfInc_p);
470 : // Initialize the co-ordinates used for reading the CF values The
471 : // CFs are 4D arrays, with the last two axis degenerate (of length
472 : // 1). The last two axis were formerly the W-, and
473 : // Polarization-axis.
474 6102 : iloc = 0;
475 6102 : Int vbSpw = (vbs.vb_p)->spectralWindows()(0);
476 6102 : Double vbPA = vbs.paQuant_p.getValue("rad");
477 :
478 6102 : Vector<Double> wVals, fVals; PolMapType mVals, mNdx, conjMVals, conjMNdx;
479 : Double fIncr, wIncr;
480 6102 : cfb->getCoordList(fVals,wVals,mNdx, mVals, conjMNdx, conjMVals, fIncr, wIncr);
481 6102 : nw = wVals.nelements();
482 :
483 504720 : for(Int irow=rbeg; irow<rend; irow++) {
484 498618 : if(!rowFlag[irow]) {
485 :
486 497652 : setFieldPhaseGrad((vb2CFBMap_p->vectorPhaseGradCalculator_p[vb2CFBMap_p->vbRow2BLMap_p[irow]])->field_phaseGrad_p);
487 497652 : cfb = (*vb2CFBMap_p)[irow];
488 :
489 995304 : for (Int ichan=0; ichan < nDataChan; ichan++) {
490 497652 : achan=chanMap_p[ichan];
491 :
492 497652 : if((achan>=0) && (achan<nGridChan)) {
493 497652 : Double dataWVal = (vbs.vb_p->uvw()(2,irow));
494 497652 : Int wndx = cfb->nearestWNdx(abs(dataWVal)*freq[ichan]/C::c);
495 497652 : Int fndx = cfb->nearestFreqNdx(vbSpw,ichan);
496 : Float s;
497 :
498 497652 : cfb->getParams(cfRefFreq,s,support(0),support(1),fndx,wndx,0);
499 497652 : sampling(0) = sampling(1) = SynthesisUtils::nint(s);
500 :
501 497652 : sgrid(pos,loc,off,phasor,irow,uvw,dphase_p[irow],freq[ichan],
502 497652 : uvwScale_p,offset_p,sampling);
503 :
504 : Bool isOnGrid;
505 :
506 : {
507 1492956 : for(Int ipol=0; ipol < nDataPol; ipol++)
508 : {
509 995304 : if(!flagCube(ipol,ichan,irow))
510 : {
511 995304 : apol=polMap_p[ipol];
512 :
513 995304 : if((apol>=0) && (apol<nGridPol))
514 : {
515 995304 : igrdpos[2]=apol; igrdpos[3]=achan;
516 995304 : nvalue=0.0; norm(ipol)=0.0;
517 :
518 : // With VBRow2CFMap in use, CF for each pol. plane is a separate 2D Array.
519 1984082 : for (uInt mCol=0; mCol<conjMNdx[ipol].nelements(); mCol++)
520 : {
521 : //
522 : int visGridElement, muellerElement;
523 : // Get the pointer to the storage for the CF
524 : // indexed by the Freq, W-term and Mueller
525 : // Element.
526 : //
527 995304 : Complex* convFuncV=NULL;
528 : try
529 : {
530 995304 : convFuncV = getConvFunc_p(vbPA,cfShape, support, muellerElement,
531 : cfb, dataWVal, fndx, wndx, conjMNdx,mNdx,
532 : ipol, mCol);
533 : }
534 0 : catch (SynthesisFTMachineError& x)
535 : {
536 0 : LogIO log_l(LogOrigin("AWVisResampler[R&D]","GridToData"));
537 0 : log_l << x.getMesg() << LogIO::EXCEPTION;
538 0 : }
539 : // Set the polarization plane of the gridded data to use for predicting with the CF from mCols column
540 995304 : visGridElement=(int)(muellerElement%nDataPol);
541 995304 : igrdpos[2]=polMap_p[visGridElement];
542 : //
543 : // Compute the incrmenets and center pixel for the current CF
544 : //
545 995304 : if ((isOnGrid=onGrid(nx, ny, nw, loc, support))==false) break;
546 988778 : cacheAxisIncrements(cfShape, cfInc_p);
547 :
548 1300882 : convOrigin = ((cfShape[0]%2==0) && (cfShape[1]%2==0)) ? (cfShape)/2 : cfShape/2+1;
549 :
550 : #include <synthesis/TransformMachines2/accumulateFromGrid.inc>
551 : }
552 995304 : if (norm[ipol] != Complex(0.0)) visCube(ipol,ichan,irow)=nvalue/norm[ipol]; // Goes with FortranizedLoopsFromGrid.cc
553 : }
554 : }
555 : }
556 : }
557 : }
558 : }
559 : }
560 : } // End row-loop
561 6102 : }
562 : //
563 : //-----------------------------------------------------------------------------------
564 : //
565 3767304 : void AWVisResampler::sgrid(Vector<Double>& pos, Vector<Int>& loc,
566 : Vector<Double>& off, Complex& phasor,
567 : const Int& irow, const Matrix<Double>& uvw,
568 : const Double& dphase, const Double& freq,
569 : const Vector<Double>& scale, //float in HPG
570 : const Vector<Double>& offset,
571 : const Vector<Float>& sampling)
572 : {
573 : // inv_lambda float in HPG Double here
574 : Double phase;
575 3767304 : Vector<Double> uvw_l(3,0); // This allows gridding of weights
576 : // centered on the uv-origin
577 12913401 : if (uvw.nelements() > 0) for(Int i=0;i<3;i++) uvw_l[i]=uvw(i,irow);
578 :
579 3767304 : pos(2)=sqrt(abs(scale[2]*uvw_l(2)*freq/C::c))+offset[2];
580 : //loc(2)=SynthesisUtils::nint(pos[2]);
581 3767304 : loc(2)=std::lrint(pos[2]);
582 3767304 : off(2)=0;
583 :
584 11301912 : for(Int idim=0;idim<2;idim++)
585 : {
586 7534608 : pos[idim]=scale[idim]*uvw_l(idim)*freq/C::c+(offset[idim]);
587 :
588 7534608 : loc[idim]=std::lrint(pos[idim]);
589 :
590 7534608 : off[idim]=std::lrint((loc[idim]-pos[idim])*sampling[idim]);
591 : }
592 :
593 3767304 : if (dphase != 0.0)
594 : {
595 3733725 : phase=-2.0*C::pi*dphase*freq/C::c;
596 : Double sp,cp;
597 3733725 : SINCOS(phase,sp,cp);
598 3733725 : phasor=Complex(cp,sp);
599 : }
600 : else
601 33579 : phasor=Complex(1.0);
602 3767304 : }
603 : using namespace casacore;
604 : };// end namespace casa
|
