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 0 : 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 0 : if (wVal > 0.0)
115 : {
116 0 : 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 0 : 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 0 : cfcell=&(*(cfb->getCFCellPtr(fndx,wndx,pndx)));
133 :
134 : //cerr << "CF Name: " << cfcell->fileName_p << endl;
135 :
136 0 : convFuncV = &(*cfcell->getStorage());
137 0 : 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 0 : 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 0 : 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 0 : 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 0 : muellerElement=cfcell->muellerElement_p;
172 :
173 : // cfShape.assign(cfcell->cfShape_p);
174 : //runTimeG1_p += timer_p.real();
175 :
176 :
177 0 : 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 0 : 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 0 : Vector<Float> sampling(2),scaledSampling(2);
244 0 : Vector<Int> support(2),loc(3), iloc(4),scaledSupport(2);
245 0 : Vector<Double> pos(3), off(3);
246 0 : Vector<Int> igrdpos(4);
247 :
248 0 : Complex phasor, nvalue /*, wt */;
249 0 : DComplex norm;
250 0 : Vector<Int> cfShape;
251 : // cfShape=(*vb2CFBMap_p)[0]->getStorage()(0,0,0)->getStorage()->shape().asVector();
252 :
253 0 : Vector<Int> convOrigin;// = (cfShape)/2;
254 : Double cfRefFreq;
255 0 : const Matrix<Double> UVW=vbs.uvw_p;
256 :
257 : // Double cfScale=1.0;
258 :
259 : //timer_p.mark();
260 0 : rbeg = 0; rend = vbs.nRow_p;
261 0 : rbeg = vbs.beginRow_p;
262 0 : rend = vbs.endRow_p;
263 :
264 0 : nx = grid.shape()[0]; ny = grid.shape()[1];
265 0 : nGridPol = grid.shape()[2]; nGridChan = grid.shape()[3];
266 :
267 0 : nDataPol = vbs.flagCube_p.shape()[0];
268 0 : nDataChan = vbs.flagCube_p.shape()[1];
269 :
270 : Bool Dummy, gDummy,
271 0 : accumCFs=((UVW.nelements() == 0) && dopsf);
272 :
273 0 : T* __restrict__ gridStore = grid.getStorage(gDummy);
274 :
275 0 : Double *freq=vbs.freq_p.getStorage(Dummy);
276 :
277 0 : cacheAxisIncrements(grid.shape().asVector(), gridInc_p);
278 :
279 0 : Bool * __restrict__ flagCube_ptr=vbs.flagCube_p.getStorage(Dummy);
280 0 : Bool * __restrict__ rowFlag_ptr = vbs.rowFlag_p.getStorage(Dummy);;
281 0 : Float * __restrict__ imgWts_ptr = vbs.imagingWeight_p.getStorage(Dummy);
282 0 : Complex * __restrict__ visCube_ptr = vbs.visCube_p.getStorage(Dummy);
283 :
284 0 : Vector<Double> wVals, fVals; PolMapType mVals, mNdx, conjMVals, conjMNdx;
285 : Double fIncr, wIncr;
286 0 : CountedPtr<CFBuffer> cfb = (*vb2CFBMap_p)[0];
287 0 : bool finitePointingOffsets=cfb->finitePointingOffsets();
288 :
289 0 : cfb->getCoordList(fVals,wVals,mNdx, mVals, conjMNdx, conjMVals, fIncr, wIncr);
290 :
291 0 : nw = wVals.nelements();
292 0 : iloc = 0;
293 :
294 0 : IPosition shp=vbs.flagCube_p.shape();
295 0 : if (accumCFs)
296 : {
297 0 : startChan = vbs.startChan_p;
298 0 : endChan = vbs.endChan_p;
299 : }
300 : else
301 : {
302 0 : startChan = 0;
303 0 : endChan = nDataChan;
304 : }
305 :
306 0 : Int vbSpw = (vbs.vb_p)->spectralWindows()(0);
307 0 : Double vbPA = vbs.paQuant_p.getValue("rad");
308 :
309 0 : for(Int irow=rbeg; irow< rend; irow++){
310 :
311 0 : if(!(*(rowFlag_ptr+irow)))
312 : {
313 0 : setFieldPhaseGrad((vb2CFBMap_p->vectorPhaseGradCalculator_p[vb2CFBMap_p->vbRow2BLMap_p[irow]])->field_phaseGrad_p);
314 0 : cfb = (*vb2CFBMap_p)[irow];
315 :
316 0 : for(Int ichan=startChan; ichan< endChan; ichan++)
317 : {
318 0 : if (*(imgWts_ptr + ichan+irow*nDataChan)!=0.0)
319 : {
320 0 : targetIMChan=chanMap_p[ichan];
321 :
322 0 : if((targetIMChan>=0) && (targetIMChan<nGridChan))
323 : {
324 :
325 :
326 0 : Double dataWVal = (UVW.nelements() > 0) ? UVW(2,irow) : 0.0;
327 0 : Int wndx = cfb->nearestWNdx(dataWVal*freq[ichan]/C::c);
328 0 : 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 0 : cfb->getParams(cfRefFreq, s, support(0), support(1),cfFreqNdx,wndx,0);
339 :
340 0 : sampling(0) = sampling(1) = SynthesisUtils::nint(s);
341 :
342 0 : sgrid(pos,loc,off, phasor, irow, UVW, dphase_p[irow], freq[ichan],
343 0 : uvwScale_p, offset_p, sampling);
344 : {
345 : // Loop over all image-plane polarization planes.
346 :
347 0 : for(Int ipol=0; ipol< nDataPol; ipol++)
348 : {
349 0 : if((!(*(flagCube_ptr + ipol + ichan*nDataPol + irow*nDataPol*nDataChan))))
350 : {
351 0 : targetIMPol=polMap_p(ipol);
352 0 : if ((targetIMPol>=0) && (targetIMPol<nGridPol))
353 : {
354 0 : igrdpos[2]=targetIMPol; igrdpos[3]=targetIMChan;
355 :
356 0 : norm = 0.0;
357 : // Loop over all relevant elements of the Mueller matrix for the polarization
358 : // ipol.
359 0 : Vector<int> conjMRow = conjMNdx[ipol];
360 :
361 0 : for (uInt mCols=0;mCols<conjMRow.nelements(); mCols++)
362 : {
363 0 : int visVecElement=mCols, muellerElement;
364 :
365 0 : Complex* convFuncV=NULL;
366 : //timer_p.mark();
367 : try
368 : {
369 0 : 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 0 : visVecElement=(int)(muellerElement%nDataPol);
382 : // If the vis. vector element is flagged, don't grid it.
383 0 : if(((*(flagCube_ptr + visVecElement + ichan*nDataPol + irow*nDataPol*nDataChan)))) break;
384 :
385 0 : if(dopsf) nvalue=Complex(*(imgWts_ptr + ichan + irow*nDataChan));
386 0 : else nvalue=Complex(*(imgWts_ptr+ichan+irow*nDataChan))*
387 0 : (*(visCube_ptr+visVecElement+ichan*nDataPol+irow*nDataChan*nDataPol)*phasor);
388 :
389 0 : if (!onGrid(nx, ny, nw, loc, support)) break;
390 :
391 :
392 0 : if(cfShape[0]%2==0 && cfShape[1]%2==0)
393 0 : convOrigin=cfShape/2;
394 : else
395 0 : convOrigin=cfShape/2+1;
396 0 : cacheAxisIncrements(cfShape, cfInc_p);
397 0 : nVisGridded_p++;
398 : #include <synthesis/TransformMachines2/accumulateToGrid.inc>
399 : }
400 0 : sumwt(targetIMPol,targetIMChan) += vbs.imagingWeight_p(ichan, irow)*fabs(norm);
401 : // *(sumWt_ptr+apol+achan*nGridChan)+= *(imgWts_ptr+ichan+irow*nDataChan);
402 0 : }
403 : }
404 : } // End poln-loop
405 : }
406 : }
407 : }
408 : } // End chan-loop
409 : }
410 : } // End row-loop
411 :
412 0 : T *tt=(T *)gridStore;
413 0 : grid.putStorage(tt,gDummy);
414 0 : }
415 : //
416 : //-----------------------------------------------------------------------------------
417 : // Re-sample VisBuffer to a regular grid (griddedData) (a.k.a. de-gridding)
418 : //
419 0 : 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 0 : Vector<Float> sampling(2);//scaledSampling(2);
424 0 : Vector<Int> support(2),loc(3), iloc(4);
425 0 : Vector<Double> pos(3), off(3);
426 :
427 0 : IPosition grdpos(4);
428 :
429 0 : Vector<Complex> norm(4,0.0);
430 0 : Complex phasor, nvalue;
431 0 : CountedPtr<CFBuffer> cfb=(*vb2CFBMap_p)[0];
432 0 : Vector<Int> cfShape=cfb->getStorage()(0,0,0)->getStorage()->shape().asVector();
433 0 : Bool finitePointingOffset=cfb->finitePointingOffsets();
434 :
435 :
436 0 : Vector<Int> convOrigin; // = ((cfShape[0]%2==0) && (cfShape[1]%2==0)) ? (cfShape)/2 : cfShape/2+1;
437 : Double cfRefFreq;//cfScale=1.0
438 :
439 0 : rbeg=0;
440 0 : rend=vbs.nRow_p;
441 0 : rbeg = vbs.beginRow_p;
442 0 : rend = vbs.endRow_p;
443 0 : nx = grid.shape()[0]; ny = grid.shape()[1];
444 0 : nGridPol = grid.shape()[2]; nGridChan = grid.shape()[3];
445 :
446 0 : nDataPol = vbs.flagCube_p.shape()[0];
447 0 : 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 0 : const Complex* __restrict__ gridStore = grid.getStorage(Dummy);
457 : (void)gridStore;
458 0 : Vector<Int> igrdpos(4);
459 0 : Double *freq=vbs.freq_p.getStorage(Dummy);
460 0 : Bool *rowFlag=vbs.rowFlag_p.getStorage(Dummy);
461 :
462 0 : Matrix<Double>& uvw=vbs.uvw_p;
463 0 : Cube<Complex>& visCube=vbs.visCube_p;
464 0 : Cube<Bool>& flagCube=vbs.flagCube_p;
465 :
466 0 : Vector<Int> gridInc, cfInc;
467 :
468 0 : 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 0 : iloc = 0;
475 0 : Int vbSpw = (vbs.vb_p)->spectralWindows()(0);
476 0 : Double vbPA = vbs.paQuant_p.getValue("rad");
477 :
478 0 : Vector<Double> wVals, fVals; PolMapType mVals, mNdx, conjMVals, conjMNdx;
479 : Double fIncr, wIncr;
480 0 : cfb->getCoordList(fVals,wVals,mNdx, mVals, conjMNdx, conjMVals, fIncr, wIncr);
481 0 : nw = wVals.nelements();
482 :
483 0 : for(Int irow=rbeg; irow<rend; irow++) {
484 0 : if(!rowFlag[irow]) {
485 :
486 0 : setFieldPhaseGrad((vb2CFBMap_p->vectorPhaseGradCalculator_p[vb2CFBMap_p->vbRow2BLMap_p[irow]])->field_phaseGrad_p);
487 0 : cfb = (*vb2CFBMap_p)[irow];
488 :
489 0 : for (Int ichan=0; ichan < nDataChan; ichan++) {
490 0 : achan=chanMap_p[ichan];
491 :
492 0 : if((achan>=0) && (achan<nGridChan)) {
493 0 : Double dataWVal = (vbs.vb_p->uvw()(2,irow));
494 0 : Int wndx = cfb->nearestWNdx(abs(dataWVal)*freq[ichan]/C::c);
495 0 : Int fndx = cfb->nearestFreqNdx(vbSpw,ichan);
496 : Float s;
497 :
498 0 : cfb->getParams(cfRefFreq,s,support(0),support(1),fndx,wndx,0);
499 0 : sampling(0) = sampling(1) = SynthesisUtils::nint(s);
500 :
501 0 : sgrid(pos,loc,off,phasor,irow,uvw,dphase_p[irow],freq[ichan],
502 0 : uvwScale_p,offset_p,sampling);
503 :
504 : Bool isOnGrid;
505 :
506 : {
507 0 : for(Int ipol=0; ipol < nDataPol; ipol++)
508 : {
509 0 : if(!flagCube(ipol,ichan,irow))
510 : {
511 0 : apol=polMap_p[ipol];
512 :
513 0 : if((apol>=0) && (apol<nGridPol))
514 : {
515 0 : igrdpos[2]=apol; igrdpos[3]=achan;
516 0 : nvalue=0.0; norm(ipol)=0.0;
517 :
518 : // With VBRow2CFMap in use, CF for each pol. plane is a separate 2D Array.
519 0 : 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 0 : Complex* convFuncV=NULL;
528 : try
529 : {
530 0 : 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 0 : visGridElement=(int)(muellerElement%nDataPol);
541 0 : igrdpos[2]=polMap_p[visGridElement];
542 : //
543 : // Compute the incrmenets and center pixel for the current CF
544 : //
545 0 : if ((isOnGrid=onGrid(nx, ny, nw, loc, support))==false) break;
546 0 : cacheAxisIncrements(cfShape, cfInc_p);
547 :
548 0 : convOrigin = ((cfShape[0]%2==0) && (cfShape[1]%2==0)) ? (cfShape)/2 : cfShape/2+1;
549 :
550 : #include <synthesis/TransformMachines2/accumulateFromGrid.inc>
551 : }
552 0 : 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 0 : }
562 : //
563 : //-----------------------------------------------------------------------------------
564 : //
565 0 : 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 0 : Vector<Double> uvw_l(3,0); // This allows gridding of weights
576 : // centered on the uv-origin
577 0 : if (uvw.nelements() > 0) for(Int i=0;i<3;i++) uvw_l[i]=uvw(i,irow);
578 :
579 0 : pos(2)=sqrt(abs(scale[2]*uvw_l(2)*freq/C::c))+offset[2];
580 : //loc(2)=SynthesisUtils::nint(pos[2]);
581 0 : loc(2)=std::lrint(pos[2]);
582 0 : off(2)=0;
583 :
584 0 : for(Int idim=0;idim<2;idim++)
585 : {
586 0 : pos[idim]=scale[idim]*uvw_l(idim)*freq/C::c+(offset[idim]);
587 :
588 0 : loc[idim]=std::lrint(pos[idim]);
589 :
590 0 : off[idim]=std::lrint((loc[idim]-pos[idim])*sampling[idim]);
591 : }
592 :
593 0 : if (dphase != 0.0)
594 : {
595 0 : phase=-2.0*M_PI*dphase*freq/C::c;
596 : Double sp,cp;
597 0 : SINCOS(phase,sp,cp);
598 0 : phasor=Complex(cp,sp);
599 : }
600 : else
601 0 : phasor=Complex(1.0);
602 0 : }
603 : using namespace casacore;
604 : };// end namespace casa
|
