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 968 : Array<Complex> tt=SynthesisUtils::getCFPixels(cfb->getCFCacheDir(), cfcell->fileName_p);
149 968 : cfcell->setStorage(tt);
150 :
151 : //cerr << (cfcell->isRotationallySymmetric_p?"o":"+");
152 :
153 : // No rotation necessary if the CF is rotationally symmetric
154 968 : if (!(cfcell->isRotationallySymmetric_p))
155 : {
156 920 : CFCell *baseCFC=NULL;
157 : // Rotate only if the difference between CF PA and VB PA
158 : // is greater than paTolerance.
159 920 : SynthesisUtils::rotate2(vbPA, *baseCFC, *cfcell, paTolerance_p);
160 : }
161 968 : convFuncV = &(*cfcell->getStorage());
162 968 : }
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 : //
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 0 : }
228 0 : }
229 :
230 : //
231 : //-----------------------------------------------------------------------------------
232 : // Template implementation for DataToGrid
233 : //
234 : template <class T>
235 28083 : 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 28083 : Vector<Float> sampling(2),scaledSampling(2);
245 28083 : Vector<Int> support(2),loc(3), iloc(4),scaledSupport(2);
246 28083 : Vector<Double> pos(3), off(3);
247 28083 : Vector<Int> igrdpos(4);
248 :
249 28083 : Complex phasor, nvalue /*, wt */;
250 28083 : DComplex norm;
251 28083 : Vector<Int> cfShape;
252 : // cfShape=(*vb2CFBMap_p)[0]->getStorage()(0,0,0)->getStorage()->shape().asVector();
253 :
254 28083 : Vector<Int> convOrigin;// = (cfShape)/2;
255 : Double cfRefFreq;
256 28083 : const Matrix<Double> UVW=vbs.uvw_p;
257 :
258 : // Double cfScale=1.0;
259 :
260 : //timer_p.mark();
261 28083 : rbeg = 0; rend = vbs.nRow_p;
262 28083 : rbeg = vbs.beginRow_p;
263 28083 : rend = vbs.endRow_p;
264 :
265 28083 : nx = grid.shape()[0]; ny = grid.shape()[1];
266 28083 : nGridPol = grid.shape()[2]; nGridChan = grid.shape()[3];
267 :
268 28083 : nDataPol = vbs.flagCube_p.shape()[0];
269 28083 : nDataChan = vbs.flagCube_p.shape()[1];
270 :
271 : Bool Dummy, gDummy,
272 28083 : accumCFs=((UVW.nelements() == 0) && dopsf);
273 :
274 28083 : T* __restrict__ gridStore = grid.getStorage(gDummy);
275 :
276 28083 : Double *freq=vbs.freq_p.getStorage(Dummy);
277 :
278 28083 : cacheAxisIncrements(grid.shape().asVector(), gridInc_p);
279 :
280 28083 : Bool * __restrict__ flagCube_ptr=vbs.flagCube_p.getStorage(Dummy);
281 28083 : Bool * __restrict__ rowFlag_ptr = vbs.rowFlag_p.getStorage(Dummy);;
282 28083 : Float * __restrict__ imgWts_ptr = vbs.imagingWeight_p.getStorage(Dummy);
283 28083 : Complex * __restrict__ visCube_ptr = vbs.visCube_p.getStorage(Dummy);
284 :
285 28083 : Vector<Double> wVals, fVals; PolMapType mVals, mNdx, conjMVals, conjMNdx;
286 : Double fIncr, wIncr;
287 28083 : CountedPtr<CFBuffer> cfb = (*vb2CFBMap_p)[0];
288 28083 : bool finitePointingOffsets=cfb->finitePointingOffsets();
289 :
290 28083 : cfb->getCoordList(fVals,wVals,mNdx, mVals, conjMNdx, conjMVals, fIncr, wIncr);
291 :
292 28083 : nw = wVals.nelements();
293 28083 : iloc = 0;
294 :
295 28083 : IPosition shp=vbs.flagCube_p.shape();
296 28083 : if (accumCFs)
297 : {
298 7150 : startChan = vbs.startChan_p;
299 7150 : endChan = vbs.endChan_p;
300 : }
301 : else
302 : {
303 20933 : startChan = 0;
304 20933 : endChan = nDataChan;
305 : }
306 :
307 28083 : Int vbSpw = (vbs.vb_p)->spectralWindows()(0);
308 28083 : Double vbPA = vbs.paQuant_p.getValue("rad");
309 :
310 3301599 : for(Int irow=rbeg; irow< rend; irow++){
311 :
312 3273516 : if(!(*(rowFlag_ptr+irow)))
313 : {
314 3269652 : setFieldPhaseGrad((vb2CFBMap_p->vectorPhaseGradCalculator_p[vb2CFBMap_p->vbRow2BLMap_p[irow]])->field_phaseGrad_p);
315 3269652 : cfb = (*vb2CFBMap_p)[irow];
316 :
317 6539304 : for(Int ichan=startChan; ichan< endChan; ichan++)
318 : {
319 3269652 : if (*(imgWts_ptr + ichan+irow*nDataChan)!=0.0)
320 : {
321 3269652 : targetIMChan=chanMap_p[ichan];
322 :
323 3269652 : if((targetIMChan>=0) && (targetIMChan<nGridChan))
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 6501368 : if(cfShape[0]%2==0 && cfShape[1]%2==0)
392 1155072 : convOrigin=cfShape/2;
393 : else
394 5346296 : convOrigin=cfShape/2+1;
395 6501368 : cacheAxisIncrements(cfShape, cfInc_p);
396 6501368 : nVisGridded_p++;
397 : #include <synthesis/TransformMachines2/accumulateToGrid.inc>
398 : }
399 6539304 : sumwt(targetIMPol,targetIMChan) += vbs.imagingWeight_p(ichan, irow)*fabs(norm);
400 : // *(sumWt_ptr+apol+achan*nGridChan)+= *(imgWts_ptr+ichan+irow*nDataChan);
401 6539304 : }
402 : }
403 : } // End poln-loop
404 : }
405 : }
406 : }
407 : } // End chan-loop
408 : }
409 : } // End row-loop
410 :
411 28083 : T *tt=(T *)gridStore;
412 28083 : grid.putStorage(tt,gDummy);
413 28083 : }
414 : //
415 : //-----------------------------------------------------------------------------------
416 : // Re-sample VisBuffer to a regular grid (griddedData) (a.k.a. de-gridding)
417 : //
418 6102 : 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 6102 : Vector<Float> sampling(2);//scaledSampling(2);
423 6102 : Vector<Int> support(2),loc(3), iloc(4);
424 6102 : Vector<Double> pos(3), off(3);
425 :
426 6102 : IPosition grdpos(4);
427 :
428 6102 : Vector<Complex> norm(4,0.0);
429 6102 : Complex phasor, nvalue;
430 6102 : CountedPtr<CFBuffer> cfb=(*vb2CFBMap_p)[0];
431 6102 : Vector<Int> cfShape=cfb->getStorage()(0,0,0)->getStorage()->shape().asVector();
432 6102 : Bool finitePointingOffset=cfb->finitePointingOffsets();
433 :
434 6102 : Vector<Int> convOrigin; // = ((cfShape[0]%2==0) && (cfShape[1]%2==0)) ? (cfShape)/2 : cfShape/2+1;
435 : Double cfRefFreq;//cfScale=1.0
436 :
437 6102 : rbeg=0;
438 6102 : rend=vbs.nRow_p;
439 6102 : rbeg = vbs.beginRow_p;
440 6102 : rend = vbs.endRow_p;
441 6102 : nx = grid.shape()[0]; ny = grid.shape()[1];
442 6102 : nGridPol = grid.shape()[2]; nGridChan = grid.shape()[3];
443 :
444 6102 : nDataPol = vbs.flagCube_p.shape()[0];
445 6102 : 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 6102 : const Complex* __restrict__ gridStore = grid.getStorage(Dummy);
455 : (void)gridStore;
456 6102 : Vector<Int> igrdpos(4);
457 6102 : Double *freq=vbs.freq_p.getStorage(Dummy);
458 6102 : Bool *rowFlag=vbs.rowFlag_p.getStorage(Dummy);
459 :
460 6102 : Matrix<Double>& uvw=vbs.uvw_p;
461 6102 : Cube<Complex>& visCube=vbs.visCube_p;
462 6102 : Cube<Bool>& flagCube=vbs.flagCube_p;
463 :
464 6102 : Vector<Int> gridInc, cfInc;
465 :
466 6102 : 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 6102 : iloc = 0;
473 6102 : Int vbSpw = (vbs.vb_p)->spectralWindows()(0);
474 6102 : Double vbPA = vbs.paQuant_p.getValue("rad");
475 :
476 6102 : Vector<Double> wVals, fVals; PolMapType mVals, mNdx, conjMVals, conjMNdx;
477 : Double fIncr, wIncr;
478 6102 : cfb->getCoordList(fVals,wVals,mNdx, mVals, conjMNdx, conjMVals, fIncr, wIncr);
479 6102 : nw = wVals.nelements();
480 :
481 504720 : for(Int irow=rbeg; irow<rend; irow++) {
482 498618 : if(!rowFlag[irow]) {
483 :
484 497652 : setFieldPhaseGrad((vb2CFBMap_p->vectorPhaseGradCalculator_p[vb2CFBMap_p->vbRow2BLMap_p[irow]])->field_phaseGrad_p);
485 497652 : cfb = (*vb2CFBMap_p)[irow];
486 :
487 995304 : for (Int ichan=0; ichan < nDataChan; ichan++) {
488 497652 : achan=chanMap_p[ichan];
489 :
490 497652 : if((achan>=0) && (achan<nGridChan)) {
491 497652 : Double dataWVal = (vbs.vb_p->uvw()(2,irow));
492 497652 : Int wndx = cfb->nearestWNdx(abs(dataWVal)*freq[ichan]/C::c);
493 497652 : Int fndx = cfb->nearestFreqNdx(vbSpw,ichan);
494 : Float s;
495 :
496 497652 : cfb->getParams(cfRefFreq,s,support(0),support(1),fndx,wndx,0);
497 497652 : sampling(0) = sampling(1) = SynthesisUtils::nint(s);
498 :
499 497652 : sgrid(pos,loc,off,phasor,irow,uvw,dphase_p[irow],freq[ichan],
500 497652 : uvwScale_p,offset_p,sampling);
501 :
502 : Bool isOnGrid;
503 :
504 : {
505 1492956 : for(Int ipol=0; ipol < nDataPol; ipol++)
506 : {
507 995304 : if(!flagCube(ipol,ichan,irow))
508 : {
509 995304 : apol=polMap_p[ipol];
510 :
511 995304 : if((apol>=0) && (apol<nGridPol))
512 : {
513 995304 : igrdpos[2]=apol; igrdpos[3]=achan;
514 995304 : nvalue=0.0; norm(ipol)=0.0;
515 :
516 : // With VBRow2CFMap in use, CF for each pol. plane is a separate 2D Array.
517 1984082 : 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 995304 : Complex* convFuncV=NULL;
526 : try
527 : {
528 995304 : 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 0 : }
537 : // Set the polarization plane of the gridded data to use for predicting with the CF from mCols column
538 995304 : visGridElement=(int)(muellerElement%nDataPol);
539 995304 : igrdpos[2]=polMap_p[visGridElement];
540 : //
541 : // Compute the incrmenets and center pixel for the current CF
542 : //
543 995304 : if ((isOnGrid=onGrid(nx, ny, nw, loc, support))==false) break;
544 988778 : cacheAxisIncrements(cfShape, cfInc_p);
545 1300882 : convOrigin = ((cfShape[0]%2==0) && (cfShape[1]%2==0)) ? (cfShape)/2 : cfShape/2+1;
546 :
547 : #include <synthesis/TransformMachines2/accumulateFromGrid.inc>
548 : }
549 995304 : 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 6102 : }
559 : //
560 : //-----------------------------------------------------------------------------------
561 : //
562 3767304 : 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 3767304 : Vector<Double> uvw_l(3,0); // This allows gridding of weights
573 : // centered on the uv-origin
574 12913401 : if (uvw.nelements() > 0) for(Int i=0;i<3;i++) uvw_l[i]=uvw(i,irow);
575 :
576 3767304 : pos(2)=sqrt(abs(scale[2]*uvw_l(2)*freq/C::c))+offset[2];
577 : //loc(2)=SynthesisUtils::nint(pos[2]);
578 3767304 : loc(2)=std::lrint(pos[2]);
579 3767304 : off(2)=0;
580 :
581 11301912 : for(Int idim=0;idim<2;idim++)
582 : {
583 7534608 : pos[idim]=scale[idim]*uvw_l(idim)*freq/C::c+(offset[idim]);
584 :
585 7534608 : loc[idim]=std::lrint(pos[idim]);
586 :
587 7534608 : off[idim]=std::lrint((loc[idim]-pos[idim])*sampling[idim]);
588 : }
589 :
590 3767304 : if (dphase != 0.0)
591 : {
592 3733725 : phase=-2.0*C::pi*dphase*freq/C::c;
593 : Double sp,cp;
594 3733725 : SINCOS(phase,sp,cp);
595 3733725 : phasor=Complex(cp,sp);
596 : }
597 : else
598 33579 : phasor=Complex(1.0);
599 3767304 : }
600 : using namespace casacore;
601 : };// end namespace casa
|
