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: aips2-request@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : //#
27 : //# $Id$
28 :
29 : #include <synthesis/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 7256616 : 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 7256616 : if (wVal > 0.0)
115 : {
116 2170326 : 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 5086290 : 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 7256616 : cfcell=&(*(cfb->getCFCellPtr(fndx,wndx,pndx)));
133 :
134 : //cerr << "CF Name: " << cfcell->fileName_p << endl;
135 :
136 7256616 : convFuncV = &(*cfcell->getStorage());
137 7256616 : 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 7256616 : 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 7256616 : if (convFuncV->shape().product() == 0)
147 : {
148 932 : Array<Complex> tt=SynthesisUtils::getCFPixels(cfb->getCFCacheDir(), cfcell->fileName_p);
149 932 : cfcell->setStorage(tt);
150 :
151 : //cerr << (cfcell->isRotationallySymmetric_p?"o":"+");
152 :
153 : // No rotation necessary if the CF is rotationally symmetric
154 932 : if (!(cfcell->isRotationallySymmetric_p))
155 : {
156 884 : CFCell *baseCFC=NULL;
157 : // Rotate only if the difference between CF PA and VB PA
158 : // is greater than paTolerance.
159 884 : SynthesisUtils::rotate2(vbPA, *baseCFC, *cfcell, paTolerance_p);
160 : }
161 932 : convFuncV = &(*cfcell->getStorage());
162 : }
163 :
164 : //cfShape.reference(cfcell->cfShape_p);
165 7256616 : 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 7256616 : muellerElement=cfcell->muellerElement_p;
172 :
173 : // cfShape.assign(cfcell->cfShape_p);
174 : //runTimeG1_p += timer_p.real();
175 :
176 :
177 14513232 : return convFuncV->getStorage(Dummy);
178 : };
179 :
180 :
181 : //
182 : //-----------------------------------------------------------------------------------
183 : //
184 0 : void AWVisResampler::cachePhaseGrad_p(const Vector<Double>& pointingOffset,
185 : const Vector<Int>&cfShape,
186 : const Vector<Int>& convOrigin,
187 : const Double& /*cfRefFreq*/,
188 : const Double& /*imRefFreq*/,
189 : const Int& spwID, const Int& fieldId)
190 : {
191 0 : if (
192 0 : ((fabs(pointingOffset[0]-cached_PointingOffset_p[0])) > 1e-6) ||
193 0 : ((fabs(pointingOffset[1]-cached_PointingOffset_p[1])) > 1e-6) ||
194 0 : (cached_phaseGrad_p.shape()[0] < cfShape[0]) ||
195 0 : (cached_phaseGrad_p.shape()[1] < cfShape[1])
196 : )
197 : {
198 0 : LogIO log_l(LogOrigin("AWVisResampler","cachePhaseGrad[R&D]"));
199 : log_l << "Computing phase gradiant for pointing offset "
200 : << pointingOffset << cfShape << " " << cached_phaseGrad_p.shape()
201 : << "(SPW: " << spwID << " Field: " << fieldId << ")"
202 : << LogIO::DEBUGGING
203 0 : << LogIO::POST;
204 0 : Int nx=cfShape(0), ny=cfShape(1);
205 : Double grad;
206 0 : Complex phx,phy;
207 :
208 0 : cached_phaseGrad_p.resize(nx,ny);
209 0 : cached_PointingOffset_p = pointingOffset;
210 :
211 0 : for(Int ix=0;ix<nx;ix++)
212 : {
213 0 : grad = (ix-convOrigin[0])*pointingOffset[0];
214 : Double sx,cx;
215 0 : SINCOS(grad,sx,cx);
216 :
217 0 : phx = Complex(cx,sx);
218 0 : for(Int iy=0;iy<ny;iy++)
219 : {
220 0 : grad = (iy-convOrigin[1])*pointingOffset[1];
221 : Double sy,cy;
222 0 : SINCOS(grad,sy,cy);
223 0 : phy = Complex(cy,sy);
224 0 : cached_phaseGrad_p(ix,iy)=phx*phy;
225 : }
226 : }
227 : }
228 0 : }
229 :
230 : //
231 : //-----------------------------------------------------------------------------------
232 : // Template implementation for DataToGrid
233 : //
234 : template <class T>
235 27831 : void AWVisResampler::DataToGridImpl_p(Array<T>& grid, VBStore& vbs,
236 : Matrix<Double>& sumwt,const Bool& dopsf,
237 : Bool /*useConjFreqCF*/)
238 : {
239 : Int nDataChan, nDataPol, nGridPol, nGridChan, nx, ny, nw;//, nCFFreq;
240 : Int targetIMChan, targetIMPol, rbeg, rend;//, PolnPlane, ConjPlane;
241 : Int startChan, endChan;
242 :
243 :
244 55662 : Vector<Float> sampling(2),scaledSampling(2);
245 55662 : Vector<Int> support(2),loc(3), iloc(4),scaledSupport(2);
246 55662 : Vector<Double> pos(3), off(3);
247 55662 : Vector<Int> igrdpos(4);
248 :
249 27831 : Complex phasor, nvalue /*, wt */;
250 27831 : DComplex norm;
251 55662 : Vector<Int> cfShape;
252 : // cfShape=(*vb2CFBMap_p)[0]->getStorage()(0,0,0)->getStorage()->shape().asVector();
253 :
254 55662 : Vector<Int> convOrigin;// = (cfShape)/2;
255 : Double cfRefFreq;
256 55662 : const Matrix<Double> UVW=vbs.uvw_p;
257 :
258 : // Double cfScale=1.0;
259 :
260 : //timer_p.mark();
261 27831 : rbeg = 0; rend = vbs.nRow_p;
262 27831 : rbeg = vbs.beginRow_p;
263 27831 : rend = vbs.endRow_p;
264 :
265 27831 : nx = grid.shape()[0]; ny = grid.shape()[1];
266 27831 : nGridPol = grid.shape()[2]; nGridChan = grid.shape()[3];
267 :
268 27831 : nDataPol = vbs.flagCube_p.shape()[0];
269 27831 : nDataChan = vbs.flagCube_p.shape()[1];
270 :
271 : Bool Dummy, gDummy,
272 27831 : accumCFs=((UVW.nelements() == 0) && dopsf);
273 :
274 27831 : T* __restrict__ gridStore = grid.getStorage(gDummy);
275 :
276 27831 : Double *freq=vbs.freq_p.getStorage(Dummy);
277 :
278 27831 : cacheAxisIncrements(grid.shape().asVector(), gridInc_p);
279 :
280 27831 : Bool * __restrict__ flagCube_ptr=vbs.flagCube_p.getStorage(Dummy);
281 27831 : Bool * __restrict__ rowFlag_ptr = vbs.rowFlag_p.getStorage(Dummy);;
282 27831 : Float * __restrict__ imgWts_ptr = vbs.imagingWeight_p.getStorage(Dummy);
283 27831 : Complex * __restrict__ visCube_ptr = vbs.visCube_p.getStorage(Dummy);
284 :
285 55662 : Vector<Double> wVals, fVals; PolMapType mVals, mNdx, conjMVals, conjMNdx;
286 : Double fIncr, wIncr;
287 55662 : CountedPtr<CFBuffer> cfb = (*vb2CFBMap_p)[0];
288 27831 : bool finitePointingOffsets=cfb->finitePointingOffsets();
289 :
290 27831 : cfb->getCoordList(fVals,wVals,mNdx, mVals, conjMNdx, conjMVals, fIncr, wIncr);
291 :
292 27831 : nw = wVals.nelements();
293 27831 : iloc = 0;
294 :
295 55662 : IPosition shp=vbs.flagCube_p.shape();
296 27831 : if (accumCFs)
297 : {
298 7114 : startChan = vbs.startChan_p;
299 7114 : endChan = vbs.endChan_p;
300 : }
301 : else
302 : {
303 20717 : startChan = 0;
304 20717 : endChan = nDataChan;
305 : }
306 :
307 27831 : Int vbSpw = (vbs.vb_p)->spectralWindows()(0);
308 27831 : Double vbPA = vbs.paQuant_p.getValue("rad");
309 :
310 3212895 : for(Int irow=rbeg; irow< rend; irow++){
311 :
312 3185064 : if(!(*(rowFlag_ptr+irow)))
313 : {
314 3181200 : setFieldPhaseGrad((vb2CFBMap_p->vectorPhaseGradCalculator_p[vb2CFBMap_p->vbRow2BLMap_p[irow]])->field_phaseGrad_p);
315 3181200 : cfb = (*vb2CFBMap_p)[irow];
316 :
317 6362400 : for(Int ichan=startChan; ichan< endChan; ichan++)
318 : {
319 3181200 : if (*(imgWts_ptr + ichan+irow*nDataChan)!=0.0)
320 : {
321 3181200 : targetIMChan=chanMap_p[ichan];
322 :
323 3181200 : if((targetIMChan>=0) && (targetIMChan<nGridChan))
324 : {
325 :
326 3181200 : Double dataWVal = (UVW.nelements() > 0) ? UVW(2,irow) : 0.0;
327 3181200 : Int wndx = cfb->nearestWNdx(dataWVal*freq[ichan]/C::c);
328 3181200 : 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 3181200 : cfb->getParams(cfRefFreq, s, support(0), support(1),cfFreqNdx,wndx,0);
339 :
340 3181200 : sampling(0) = sampling(1) = SynthesisUtils::nint(s);
341 :
342 3181200 : sgrid(pos,loc,off, phasor, irow, UVW, dphase_p[irow], freq[ichan],
343 3181200 : uvwScale_p, offset_p, sampling);
344 : {
345 : // Loop over all image-plane polarization planes.
346 :
347 9543600 : for(Int ipol=0; ipol< nDataPol; ipol++)
348 : {
349 6362400 : if((!(*(flagCube_ptr + ipol + ichan*nDataPol + irow*nDataPol*nDataChan))))
350 : {
351 6362400 : targetIMPol=polMap_p(ipol);
352 6362400 : if ((targetIMPol>=0) && (targetIMPol<nGridPol))
353 : {
354 6362400 : igrdpos[2]=targetIMPol; igrdpos[3]=targetIMChan;
355 :
356 6362400 : norm = 0.0;
357 : // Loop over all relevant elements of the Mueller matrix for the polarization
358 : // ipol.
359 6362400 : Vector<int> conjMRow = conjMNdx[ipol];
360 :
361 12688160 : for (uInt mCols=0;mCols<conjMRow.nelements(); mCols++)
362 : {
363 6362400 : int visVecElement=mCols, muellerElement;
364 :
365 6362400 : Complex* convFuncV=NULL;
366 : //timer_p.mark();
367 : try
368 : {
369 6362400 : 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 : }
379 : // Extract the vis. vector element corresponding to the mCols column of the conjMRow row of the Mueller matrix.
380 :
381 6362400 : visVecElement=(int)(muellerElement%nDataPol);
382 : // If the vis. vector element is flagged, don't grid it.
383 6399040 : if(((*(flagCube_ptr + visVecElement + ichan*nDataPol + irow*nDataPol*nDataChan)))) break;
384 :
385 6362400 : if(dopsf) nvalue=Complex(*(imgWts_ptr + ichan + irow*nDataChan));
386 2628732 : else nvalue=Complex(*(imgWts_ptr+ichan+irow*nDataChan))*
387 5257464 : (*(visCube_ptr+visVecElement+ichan*nDataPol+irow*nDataChan*nDataPol)*phasor);
388 :
389 6362400 : if (!onGrid(nx, ny, nw, loc, support)) break;
390 :
391 6325760 : if(cfShape[0]%2==0 && cfShape[1]%2==0)
392 979464 : convOrigin=cfShape/2;
393 : else
394 5346296 : convOrigin=cfShape/2+1;
395 6325760 : cacheAxisIncrements(cfShape, cfInc_p);
396 6325760 : nVisGridded_p++;
397 : #include <synthesis/TransformMachines2/accumulateToGrid.inc>
398 : }
399 6362400 : sumwt(targetIMPol,targetIMChan) += vbs.imagingWeight_p(ichan, irow)*fabs(norm);
400 : // *(sumWt_ptr+apol+achan*nGridChan)+= *(imgWts_ptr+ichan+irow*nDataChan);
401 : }
402 : }
403 : } // End poln-loop
404 : }
405 : }
406 : }
407 : } // End chan-loop
408 : }
409 : } // End row-loop
410 :
411 27831 : T *tt=(T *)gridStore;
412 27831 : grid.putStorage(tt,gDummy);
413 27831 : }
414 : //
415 : //-----------------------------------------------------------------------------------
416 : // Re-sample VisBuffer to a regular grid (griddedData) (a.k.a. de-gridding)
417 : //
418 5958 : void AWVisResampler::GridToData(VBStore& vbs, const Array<Complex>& grid)
419 : {
420 : Int nDataChan, nDataPol, nGridPol, nGridChan, nx, ny,nw;//, nCFFreq;
421 : Int achan, apol, rbeg, rend;//, PolnPlane, ConjPlane;
422 11916 : Vector<Float> sampling(2);//scaledSampling(2);
423 11916 : Vector<Int> support(2),loc(3), iloc(4);
424 11916 : Vector<Double> pos(3), off(3);
425 :
426 11916 : IPosition grdpos(4);
427 :
428 11916 : Vector<Complex> norm(4,0.0);
429 5958 : Complex phasor, nvalue;
430 11916 : CountedPtr<CFBuffer> cfb=(*vb2CFBMap_p)[0];
431 11916 : Vector<Int> cfShape=cfb->getStorage()(0,0,0)->getStorage()->shape().asVector();
432 5958 : Bool finitePointingOffset=cfb->finitePointingOffsets();
433 :
434 11916 : Vector<Int> convOrigin; // = ((cfShape[0]%2==0) && (cfShape[1]%2==0)) ? (cfShape)/2 : cfShape/2+1;
435 : Double cfRefFreq;//cfScale=1.0
436 :
437 5958 : rbeg=0;
438 5958 : rend=vbs.nRow_p;
439 5958 : rbeg = vbs.beginRow_p;
440 5958 : rend = vbs.endRow_p;
441 5958 : nx = grid.shape()[0]; ny = grid.shape()[1];
442 5958 : nGridPol = grid.shape()[2]; nGridChan = grid.shape()[3];
443 :
444 5958 : nDataPol = vbs.flagCube_p.shape()[0];
445 5958 : nDataChan = vbs.flagCube_p.shape()[1];
446 :
447 : //
448 : // The following code reduces most array accesses to the simplest
449 : // possible to improve performance. However this made no
450 : // difference in the run-time performance compared to Vector,
451 : // Matrix and Cube indexing.
452 : //
453 : Bool Dummy;
454 5958 : const Complex* __restrict__ gridStore = grid.getStorage(Dummy);
455 : (void)gridStore;
456 11916 : Vector<Int> igrdpos(4);
457 5958 : Double *freq=vbs.freq_p.getStorage(Dummy);
458 5958 : Bool *rowFlag=vbs.rowFlag_p.getStorage(Dummy);
459 :
460 5958 : Matrix<Double>& uvw=vbs.uvw_p;
461 5958 : Cube<Complex>& visCube=vbs.visCube_p;
462 5958 : Cube<Bool>& flagCube=vbs.flagCube_p;
463 :
464 11916 : Vector<Int> gridInc, cfInc;
465 :
466 5958 : cacheAxisIncrements(grid.shape().asVector(), gridInc_p);
467 : //cacheAxisIncrements(cfShape, cfInc_p);
468 : // Initialize the co-ordinates used for reading the CF values The
469 : // CFs are 4D arrays, with the last two axis degenerate (of length
470 : // 1). The last two axis were formerly the W-, and
471 : // Polarization-axis.
472 5958 : iloc = 0;
473 5958 : Int vbSpw = (vbs.vb_p)->spectralWindows()(0);
474 5958 : Double vbPA = vbs.paQuant_p.getValue("rad");
475 :
476 11916 : Vector<Double> wVals, fVals; PolMapType mVals, mNdx, conjMVals, conjMNdx;
477 : Double fIncr, wIncr;
478 5958 : cfb->getCoordList(fVals,wVals,mNdx, mVals, conjMNdx, conjMVals, fIncr, wIncr);
479 5958 : nw = wVals.nelements();
480 :
481 454032 : for(Int irow=rbeg; irow<rend; irow++) {
482 448074 : if(!rowFlag[irow]) {
483 :
484 447108 : setFieldPhaseGrad((vb2CFBMap_p->vectorPhaseGradCalculator_p[vb2CFBMap_p->vbRow2BLMap_p[irow]])->field_phaseGrad_p);
485 447108 : cfb = (*vb2CFBMap_p)[irow];
486 :
487 894216 : for (Int ichan=0; ichan < nDataChan; ichan++) {
488 447108 : achan=chanMap_p[ichan];
489 :
490 447108 : if((achan>=0) && (achan<nGridChan)) {
491 447108 : Double dataWVal = (vbs.vb_p->uvw()(2,irow));
492 447108 : Int wndx = cfb->nearestWNdx(abs(dataWVal)*freq[ichan]/C::c);
493 447108 : Int fndx = cfb->nearestFreqNdx(vbSpw,ichan);
494 : Float s;
495 :
496 447108 : cfb->getParams(cfRefFreq,s,support(0),support(1),fndx,wndx,0);
497 447108 : sampling(0) = sampling(1) = SynthesisUtils::nint(s);
498 :
499 447108 : sgrid(pos,loc,off,phasor,irow,uvw,dphase_p[irow],freq[ichan],
500 447108 : uvwScale_p,offset_p,sampling);
501 :
502 : Bool isOnGrid;
503 :
504 : {
505 1341324 : for(Int ipol=0; ipol < nDataPol; ipol++)
506 : {
507 894216 : if(!flagCube(ipol,ichan,irow))
508 : {
509 894216 : apol=polMap_p[ipol];
510 :
511 894216 : if((apol>=0) && (apol<nGridPol))
512 : {
513 894216 : igrdpos[2]=apol; igrdpos[3]=achan;
514 894216 : nvalue=0.0; norm(ipol)=0.0;
515 :
516 : // With VBRow2CFMap in use, CF for each pol. plane is a separate 2D Array.
517 1782762 : for (uInt mCol=0; mCol<conjMNdx[ipol].nelements(); mCol++)
518 : {
519 : //
520 : int visGridElement, muellerElement;
521 : // Get the pointer to the storage for the CF
522 : // indexed by the Freq, W-term and Mueller
523 : // Element.
524 : //
525 894216 : Complex* convFuncV=NULL;
526 : try
527 : {
528 894216 : convFuncV = getConvFunc_p(vbPA,cfShape, support, muellerElement,
529 : cfb, dataWVal, fndx, wndx, conjMNdx,mNdx,
530 : ipol, mCol);
531 : }
532 0 : catch (SynthesisFTMachineError& x)
533 : {
534 0 : LogIO log_l(LogOrigin("AWVisResampler[R&D]","GridToData"));
535 0 : log_l << x.getMesg() << LogIO::EXCEPTION;
536 : }
537 : // Set the polarization plane of the gridded data to use for predicting with the CF from mCols column
538 894216 : visGridElement=(int)(muellerElement%nDataPol);
539 894216 : igrdpos[2]=polMap_p[visGridElement];
540 : //
541 : // Compute the incrmenets and center pixel for the current CF
542 : //
543 894216 : if ((isOnGrid=onGrid(nx, ny, nw, loc, support))==false) break;
544 888546 : cacheAxisIncrements(cfShape, cfInc_p);
545 888546 : convOrigin = ((cfShape[0]%2==0) && (cfShape[1]%2==0)) ? (cfShape)/2 : cfShape/2+1;
546 :
547 : #include <synthesis/TransformMachines2/accumulateFromGrid.inc>
548 : }
549 894216 : if (norm[ipol] != Complex(0.0)) visCube(ipol,ichan,irow)=nvalue/norm[ipol]; // Goes with FortranizedLoopsFromGrid.cc
550 : }
551 : }
552 : }
553 : }
554 : }
555 : }
556 : }
557 : } // End row-loop
558 5958 : }
559 : //
560 : //-----------------------------------------------------------------------------------
561 : //
562 3628308 : void AWVisResampler::sgrid(Vector<Double>& pos, Vector<Int>& loc,
563 : Vector<Double>& off, Complex& phasor,
564 : const Int& irow, const Matrix<Double>& uvw,
565 : const Double& dphase, const Double& freq,
566 : const Vector<Double>& scale, //float in HPG
567 : const Vector<Double>& offset,
568 : const Vector<Float>& sampling)
569 : {
570 : // inv_lambda float in HPG Double here
571 : Double phase;
572 7256616 : Vector<Double> uvw_l(3,0); // This allows gridding of weights
573 : // centered on the uv-origin
574 12395325 : if (uvw.nelements() > 0) for(Int i=0;i<3;i++) uvw_l[i]=uvw(i,irow);
575 :
576 3628308 : pos(2)=sqrt(abs(scale[2]*uvw_l(2)*freq/C::c))+offset[2];
577 : //loc(2)=SynthesisUtils::nint(pos[2]);
578 3628308 : loc(2)=std::lrint(pos[2]);
579 3628308 : off(2)=0;
580 :
581 10884924 : for(Int idim=0;idim<2;idim++)
582 : {
583 7256616 : pos[idim]=scale[idim]*uvw_l(idim)*freq/C::c+(offset[idim]);
584 :
585 7256616 : loc[idim]=std::lrint(pos[idim]);
586 :
587 7256616 : off[idim]=std::lrint((loc[idim]-pos[idim])*sampling[idim]);
588 : }
589 :
590 3628308 : if (dphase != 0.0)
591 : {
592 3594729 : phase=-2.0*C::pi*dphase*freq/C::c;
593 : Double sp,cp;
594 3594729 : SINCOS(phase,sp,cp);
595 3594729 : phasor=Complex(cp,sp);
596 : }
597 : else
598 33579 : phasor=Complex(1.0);
599 3628308 : }
600 : using namespace casacore;
601 : };// end namespace casa
|
