Line data Source code
1 : //# SimplePBConvFunc.cc: implementation of SimplePBConvFunc
2 : //# Copyright (C) 2007
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be adressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //#
27 : //# $Id$
28 : #include <casacore/casa/Arrays/ArrayMath.h>
29 : #include <casacore/casa/Arrays/Array.h>
30 : #include <casacore/casa/Arrays/MaskedArray.h>
31 : #include <casacore/casa/Arrays/Vector.h>
32 : #include <casacore/casa/Arrays/Slice.h>
33 : #include <casacore/casa/Arrays/Slicer.h>
34 : #include <casacore/casa/Arrays/Matrix.h>
35 : #include <casacore/casa/Arrays/Cube.h>
36 : #include <casacore/casa/OS/Timer.h>
37 : #include <casacore/casa/Utilities/Assert.h>
38 : #include <casacore/casa/Quanta/MVTime.h>
39 : #include <casacore/casa/Quanta/MVAngle.h>
40 : #include <casacore/measures/Measures/MeasTable.h>
41 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
42 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
43 :
44 : #include <casacore/images/Images/ImageInterface.h>
45 : #include <casacore/images/Images/PagedImage.h>
46 : #include <casacore/images/Images/TempImage.h>
47 : #include <casacore/images/Images/SubImage.h>
48 : #include <casacore/casa/Logging/LogIO.h>
49 : #include <casacore/casa/Logging/LogSink.h>
50 : #include <casacore/casa/Logging/LogMessage.h>
51 :
52 : #include <casacore/ms/MeasurementSets/MSColumns.h>
53 : #include <casacore/lattices/Lattices/ArrayLattice.h>
54 : #include <casacore/lattices/Lattices/SubLattice.h>
55 : #include <casacore/lattices/LRegions/LCBox.h>
56 : #include <casacore/lattices/LEL/LatticeExpr.h>
57 : #include <casacore/lattices/Lattices/LatticeCache.h>
58 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
59 :
60 : #include <casacore/scimath/Mathematics/ConvolveGridder.h>
61 : #include <casacore/scimath/Mathematics/MathFunc.h>
62 : #include <msvis/MSVis/VisBuffer.h>
63 : #include <msvis/MSVis/VisibilityIterator.h>
64 :
65 : #include <synthesis/TransformMachines/SimplePBConvFunc.h>
66 : #include <synthesis/TransformMachines/SkyJones.h>
67 : #include <synthesis/TransformMachines/PBMath1DNumeric.h>
68 : #include <casacore/scimath/Mathematics/FFTPack.h>
69 : #include <casacore/scimath/Mathematics/FFTW.h>
70 : #include <casacore/scimath/Mathematics/FFTServer.h>
71 : #include <casacore/casa/Utilities/CompositeNumber.h>
72 : #include <math.h>
73 : #ifdef _OPENMP
74 : #include <omp.h>
75 : #endif
76 : using namespace casacore;
77 : namespace casa { //# NAMESPACE CASA - BEGIN
78 :
79 0 : SimplePBConvFunc::SimplePBConvFunc(): nchan_p(-1),
80 : npol_p(-1), pointToPix_p(), directionIndex_p(-1), thePix_p(0),
81 : filledFluxScale_p(false),doneMainConv_p(0),
82 :
83 0 : calcFluxScale_p(true), actualConvIndex_p(-1), convSize_p(0), convSupport_p(0), pointingPix_p() {
84 : //
85 :
86 0 : pbClass_p=PBMathInterface::COMMONPB;
87 0 : ft_p=FFT2D(true);
88 0 : }
89 :
90 0 : SimplePBConvFunc::SimplePBConvFunc(const PBMathInterface::PBClass typeToUse):
91 : nchan_p(-1),npol_p(-1),pointToPix_p(),
92 : directionIndex_p(-1), thePix_p(0), filledFluxScale_p(false),doneMainConv_p(0),
93 0 : calcFluxScale_p(true), actualConvIndex_p(-1), convSize_p(0), convSupport_p(0), pointingPix_p() {
94 : //
95 0 : pbClass_p=typeToUse;
96 0 : ft_p=FFT2D(true);
97 0 : }
98 0 : SimplePBConvFunc::SimplePBConvFunc(const RecordInterface& rec, const Bool calcfluxneeded)
99 : : nchan_p(-1),npol_p(-1),pointToPix_p(), directionIndex_p(-1), thePix_p(0), filledFluxScale_p(false),
100 : doneMainConv_p(0),
101 0 : calcFluxScale_p(calcfluxneeded), actualConvIndex_p(-1), convSize_p(0), convSupport_p(0), pointingPix_p()
102 : {
103 0 : String err;
104 0 : fromRecord(err, rec, calcfluxneeded);
105 0 : ft_p=FFT2D(true);
106 0 : }
107 0 : SimplePBConvFunc::~SimplePBConvFunc(){
108 : //
109 :
110 :
111 0 : }
112 :
113 0 : void SimplePBConvFunc::storeImageParams(const ImageInterface<Complex>& iimage,
114 : const VisBuffer& vb){
115 : //image signature changed...rather simplistic for now
116 0 : if((iimage.shape().product() != nx_p*ny_p*nchan_p*npol_p) || nchan_p < 1){
117 0 : csys_p=iimage.coordinates();
118 0 : Int coordIndex=csys_p.findCoordinate(Coordinate::DIRECTION);
119 0 : AlwaysAssert(coordIndex>=0, AipsError);
120 0 : directionIndex_p=coordIndex;
121 0 : dc_p=csys_p.directionCoordinate(directionIndex_p);
122 0 : ObsInfo imInfo=csys_p.obsInfo();
123 0 : String tel= imInfo.telescope();
124 0 : MPosition pos;
125 0 : if (vb.msColumns().observation().nrow() > 0) {
126 0 : tel = vb.msColumns().observation().telescopeName()(vb.msColumns().observationId()(0));
127 : }
128 0 : if (tel.length() == 0 || !tel.contains("VLA") ||
129 0 : !MeasTable::Observatory(pos,tel)) {
130 : // unknown observatory, use first antenna
131 0 : Int ant1=vb.antenna1()(0);
132 0 : pos=vb.msColumns().antenna().positionMeas()(ant1);
133 : }
134 : //cout << "TELESCOPE " << tel << endl;
135 : //Store this to build epochs via the time access of visbuffer later
136 0 : timeMType_p=MEpoch::castType(vb.msColumns().timeMeas()(0).getRef().getType());
137 0 : timeUnit_p=Unit(vb.msColumns().timeMeas().measDesc().getUnits()(0).getName());
138 : // timeUnit_p=Unit("s");
139 : //cout << "UNIT " << timeUnit_p.getValue() << " name " << timeUnit_p.getName() << endl;
140 0 : pointFrame_p=MeasFrame(imInfo.obsDate(), pos);
141 0 : MDirection::Ref elRef(dc_p.directionType(), pointFrame_p);
142 : //For now we set the conversion from this direction
143 0 : pointToPix_p=MDirection::Convert( MDirection(), elRef);
144 0 : nx_p=iimage.shape()(coordIndex);
145 0 : ny_p=iimage.shape()(coordIndex+1);
146 0 : pointingPix_p.resize(nx_p, ny_p);
147 0 : pointingPix_p.set(false);
148 0 : coordIndex=csys_p.findCoordinate(Coordinate::SPECTRAL);
149 0 : Int pixAxis=csys_p.pixelAxes(coordIndex)[0];
150 0 : nchan_p=iimage.shape()(pixAxis);
151 0 : coordIndex=csys_p.findCoordinate(Coordinate::STOKES);
152 0 : pixAxis=csys_p.pixelAxes(coordIndex)[0];
153 0 : npol_p=iimage.shape()(pixAxis);
154 0 : if(calcFluxScale_p){
155 0 : if(fluxScale_p.shape().nelements()==0){
156 0 : fluxScale_p=TempImage<Float>(IPosition(4,nx_p,ny_p,npol_p,nchan_p), csys_p);
157 0 : fluxScale_p.set(0.0);
158 : }
159 0 : filledFluxScale_p=false;
160 : }
161 :
162 : }
163 :
164 0 : }
165 :
166 0 : void SimplePBConvFunc::toPix(const VisBuffer& vb){
167 0 : thePix_p.resize(2);
168 :
169 0 : const MDirection& p1=pointingDirAnt1(vb);
170 0 : if(dc_p.directionType() != MDirection::castType(p1.getRef().getType())){
171 : //pointToPix_p.setModel(theDir);
172 :
173 0 : String tel= csys_p.obsInfo().telescope();
174 0 : if(!tel.contains("VLA")) {
175 : //use first antenna as direction1_p is used to calculate pointing
176 : // as only VLA uses observatory pos for calculations
177 0 : Int ant1=vb.antenna1()(0);
178 0 : MPosition pos=vb.msColumns().antenna().positionMeas()(ant1);
179 0 : pointFrame_p.resetPosition(pos);
180 : }
181 0 : MEpoch timenow(Quantity(vb.time()(0), timeUnit_p), timeMType_p);
182 : //cerr << "Ref " << p1.toString() << " ep " << timenow.getRefString() << " time " << MVTime(timenow.getValue().getTime()).string(MVTime::YMD) << endl;
183 0 : pointFrame_p.resetEpoch(timenow);
184 : ///////////////////////////
185 : //MDirection some=pointToPix_p(vb.direction1()(0));
186 : //MVAngle mvRa=some.getAngle().getValue()(0);
187 : //MVAngle mvDec=some.getAngle().getValue()(1);
188 :
189 : //cout << mvRa(0.0).string(MVAngle::TIME,8) << " ";
190 : // cout << mvDec.string(MVAngle::DIG2,8) << " ";
191 : //cout << MDirection::showType(some.getRefPtr()->getType()) << endl;
192 :
193 : //////////////////////////
194 : //pointToPix holds pointFrame_p by reference...
195 : //thus good to go for conversion
196 0 : direction1_p=pointToPix_p(p1);
197 : //cerr << "converted " << direction1_p.toString() << endl;
198 : //direction2_p=pointToPix_p(vb.direction2()(0));
199 0 : direction2_p=direction1_p;
200 0 : dc_p.toPixel(thePix_p, direction1_p);
201 :
202 : }
203 : else{
204 0 : direction1_p=p1;
205 : //direction2_p=vb.direction2()(0);
206 : //For now
207 0 : direction2_p=direction1_p;
208 0 : dc_p.toPixel(thePix_p, direction1_p);
209 : }
210 0 : }
211 :
212 0 : void SimplePBConvFunc::setWeightImage(CountedPtr<TempImage<Float> >& wgtimage){
213 0 : convWeightImage_p=wgtimage;
214 0 : calcFluxScale_p=true;
215 :
216 0 : }
217 :
218 0 : void SimplePBConvFunc::reset(){
219 0 : doneMainConv_p.resize();
220 0 : convFunctions_p.resize(0, true);
221 0 : convWeights_p.resize(0, true);
222 0 : convSizes_p.resize(0, true);
223 0 : convSupportBlock_p.resize(0, true);
224 0 : convFunctionMap_p.clear();
225 0 : vbConvIndex_p.clear();
226 0 : }
227 :
228 :
229 :
230 0 : Int SimplePBConvFunc::convIndex(const VisBuffer& vb){
231 0 : String elkey=String::toString(vb.msId())+String("_")+String::toString(vb.spectralWindow());
232 0 : if(vbConvIndex_p.count(elkey) > 0){
233 0 : return vbConvIndex_p[elkey];
234 : }
235 0 : Int val=vbConvIndex_p.size();
236 0 : vbConvIndex_p[elkey]=val;
237 0 : return val;
238 : }
239 :
240 0 : const MDirection& SimplePBConvFunc::pointingDirAnt1(const VisBuffer& vb){
241 0 : std::ostringstream oss;
242 0 : oss << vb.msId() << "_" << vb.antenna1()(0) << "_";
243 0 : oss.precision(13);
244 0 : oss << vb.time()(0);
245 0 : String elkey=oss.str();
246 : // String elkey=String::toString(vb.msId())+String("_")+String::toString(vb.antenna1()(0))+String("_")
247 : // +String::toString(vb.time()(0));
248 :
249 : //cerr << "key " << elkey << " count " << ant1PointVal_p.count(elkey) << " size " << ant1PointVal_p.size() << " " << ant1PointingCache_p.nelements() << endl;
250 0 : if(ant1PointVal_p.count(elkey) > 0){
251 0 : return ant1PointingCache_p[ant1PointVal_p[elkey]];
252 :
253 : }
254 0 : Int val=ant1PointingCache_p.nelements();
255 0 : ant1PointingCache_p.resize(val+1, true);
256 0 : ant1PointingCache_p[val]=vb.firstDirection1();
257 : //cerr << "POINTINGDIRANT1 " << vb.firstDirection1() << " " << vb.direction1()(0) << endl;
258 : //ant1PointingCache_p[val]=vbUtil_p.getPointingDir(vb, vb.antenna1()(0), 0);
259 0 : ant1PointVal_p[elkey]=val;
260 0 : return ant1PointingCache_p[val];
261 :
262 : }
263 0 : void SimplePBConvFunc::findConvFunction(const ImageInterface<Complex>& iimage,
264 : const VisBuffer& vb,
265 : const Int& convSampling,
266 : const Vector<Double>& visFreq,
267 : Array<Complex>& convFunc,
268 : Array<Complex>& weightConvFunc,
269 : Vector<Int>& convsize,
270 : Vector<Int>& convSupport,
271 : Vector<Int>& convFuncPolMap,
272 : Vector<Int>& convFuncChanMap,
273 : Vector<Int>& convFuncRowMap,
274 : const Bool /*getConjFunc*/
275 : ){
276 :
277 :
278 :
279 0 : LogIO os;
280 0 : os << LogOrigin("SimplePBConv", "findConvFunction") << LogIO::NORMAL;
281 : /////////////////////////
282 : os<< LogIO::DEBUG1
283 0 : << "msID " << vb.msId() << " ANT1 id" << vb.antenna1()(0)
284 0 : << " direction " << vb.firstDirection1().toString() << " ANT2 id"
285 0 : << vb.antenna2()(0) << " direction " << vb.direction2()(0).toString()
286 0 : << LogIO::POST ;
287 : //////////////////////
288 0 : Int convSamp=2*convSampling;
289 0 : storeImageParams(iimage, vb);
290 0 : convFuncChanMap.resize(vb.nChannel());
291 0 : Vector<Double> beamFreqs;
292 0 : findUsefulChannels(convFuncChanMap, beamFreqs, vb, visFreq);
293 : //cerr << "CHANMAP " << convFuncChanMap << " beamFreqs " << beamFreqs << endl;
294 0 : Int nBeamChans=beamFreqs.nelements();
295 : //indgen(convFuncChanMap);
296 0 : convFuncPolMap.resize(vb.nCorr());
297 0 : convFuncPolMap.set(0);
298 : //Only one plane in this version
299 0 : convFuncRowMap.resize();
300 0 : convFuncRowMap=Vector<Int>(vb.nRow(),0);
301 : //break reference
302 0 : convFunc.resize();
303 0 : weightConvFunc.resize();
304 :
305 :
306 : // Get the coordinate system
307 0 : CoordinateSystem coords(iimage.coordinates());
308 :
309 :
310 0 : actualConvIndex_p=convIndex(vb);
311 : //cerr << "In findConv " << actualConvIndex_p << endl;
312 : // Make a two dimensional image to calculate the
313 : // primary beam. We want this on a fine grid in the
314 : // UV plane
315 0 : Int directionIndex=directionIndex_p;
316 0 : AlwaysAssert(directionIndex>=0, AipsError);
317 :
318 : // Set up the convolution function.
319 0 : Int nx=nx_p;
320 0 : Int ny=ny_p;
321 : // convSize_p=max(nx,ny)*convSampling;
322 : //cerr << "size " << nx << " " << ny << endl;
323 : //3 times the support size
324 0 : if(doneMainConv_p.shape()[0] < (actualConvIndex_p+1)){
325 : // cerr << "resizing DONEMAIN " << doneMainConv_p.shape()[0] << endl;
326 0 : doneMainConv_p.resize(actualConvIndex_p+1, true);
327 0 : doneMainConv_p[actualConvIndex_p]=false;
328 : }
329 :
330 0 : if(!(doneMainConv_p[actualConvIndex_p])){
331 :
332 : //convSize_p=4*(sj_p->support(vb, coords));
333 0 : convSize_p=Int(max(nx_p, ny_p)/2)*2*convSamp;
334 : //cerr << "convSize_p " << convSize_p << " support " << sj_p->support(vb, coords) << endl;
335 : // Make this a nice composite number, to speed up FFTs
336 : //cerr << "convSize_p 0 " << convSize_p << " convSamp " << convSamp<< endl;
337 0 : CompositeNumber cn(uInt(convSize_p*2.0));
338 :
339 0 : convSize_p = cn.nextLargerEven(Int(convSize_p));
340 : //cerr << "convSize : " << convSize_p << endl;
341 :
342 : }
343 :
344 :
345 0 : toPix(vb);
346 : //Timer tim;
347 : //tim.mark();
348 0 : addPBToFlux(vb);
349 : //tim.show("After addPBToFlux");
350 0 : DirectionCoordinate dc=dc_p;
351 :
352 : //where in the image in pixels is this pointing
353 0 : Vector<Double> pixFieldDir(2);
354 0 : pixFieldDir=thePix_p;
355 :
356 : //cerr << "pix of pointing " << pixFieldDir << endl;
357 0 : MDirection fieldDir=direction1_p;
358 : //shift from center
359 0 : pixFieldDir(0) = pixFieldDir(0) - Double(nx / 2);
360 0 : pixFieldDir(1) = pixFieldDir(1) - Double(ny / 2);
361 :
362 : //phase gradient per pixel to apply
363 0 : pixFieldDir(0) = -pixFieldDir(0)*2.0*C::pi/Double(nx)/Double(convSampling);
364 0 : pixFieldDir(1) = -pixFieldDir(1)*2.0*C::pi/Double(ny)/Double(convSampling);
365 :
366 : //cerr << "DonemainConv " << doneMainConv_p[actualConvIndex_p] << endl;
367 0 : if(!doneMainConv_p[actualConvIndex_p]){
368 0 : Vector<Double> sampling;
369 0 : sampling = dc.increment();
370 0 : sampling*=Double(convSamp);
371 0 : sampling(0)*=Double(nx)/Double(convSize_p);
372 0 : sampling(1)*=Double(ny)/Double(convSize_p);
373 0 : dc.setIncrement(sampling);
374 :
375 :
376 0 : Vector<Double> unitVec(2);
377 0 : unitVec=convSize_p/2;
378 0 : dc.setReferencePixel(unitVec);
379 :
380 :
381 : //make sure we are using the same units
382 0 : fieldDir.set(dc.worldAxisUnits()(0));
383 : //cerr << "SPB direction " << fieldDir << endl;
384 :
385 0 : dc.setReferenceValue(fieldDir.getAngle().getValue());
386 :
387 0 : coords.replaceCoordinate(dc, directionIndex);
388 0 : Int spind=coords.findCoordinate(Coordinate::SPECTRAL);
389 0 : SpectralCoordinate spCoord=coords.spectralCoordinate(spind);
390 0 : spCoord.setReferencePixel(Vector<Double>(1,0.0));
391 0 : spCoord.setReferenceValue(Vector<Double>(1, beamFreqs(0)));
392 0 : if(beamFreqs.nelements() >1)
393 0 : spCoord.setIncrement(Vector<Double>(1, beamFreqs(1)-beamFreqs(0)));
394 0 : coords.replaceCoordinate(spCoord, spind);
395 :
396 :
397 0 : CoordinateSystem coordLastPlane= coords;
398 0 : spCoord.setReferenceValue(Vector<Double>(1, beamFreqs(nBeamChans-1)));
399 0 : coordLastPlane.replaceCoordinate(spCoord, spind);
400 : //cerr << "BEAM freqs " << beamFreqs << endl;
401 :
402 : // coords.list(logIO(), MDoppler::RADIO, IPosition(), IPosition());
403 :
404 0 : Int tempConvSize=((convSize_p/4/(convSamp/convSampling))/2)*2;
405 0 : IPosition pbShape(4, tempConvSize, tempConvSize, 1, nBeamChans);
406 :
407 0 : Long memtot=HostInfo::memoryFree();
408 0 : Double memtobeused= Double(memtot)*1024.0;
409 : //cerr << "Mem to be used " << memtobeused << " arraysize " << Double(convSize_p*convSize_p)*8.0 << endl;;
410 : //check for 32 bit OS and limit it to 2Gbyte
411 : if( sizeof(void*) == 4){
412 : if(memtot > 2000000)
413 : memtot=2000000;
414 : }
415 0 : if(memtot <= 2000000)
416 0 : memtobeused=0.0;
417 : //cerr << "mem to be used " << memtobeused << endl;
418 : //tim.mark();
419 0 : IPosition start(4, 0, 0, 0, 0);
420 : //IPosition pbSlice(4, convSize_p, convSize_p, 1, 1);
421 : //cerr << "pbshape " << pbShape << endl;
422 0 : TempImage<Complex> twoDPB(TiledShape(pbShape, IPosition(4, pbShape(0), pbShape(1), 1, 1)), coords, memtobeused/10.0);
423 :
424 : //tim.show("after making one image");
425 0 : convFunc_p.resize(tempConvSize, tempConvSize);
426 0 : convFunc_p=0.0;
427 :
428 :
429 :
430 : // Accumulate terms
431 : //Matrix<Complex> screen(convSize_p, convSize_p);
432 : //screen=1.0;
433 : // Either the SkyJones
434 : //tim.mark();
435 : //twoDPB.set(Complex(1.0,0.0));
436 : //for (Int k=0; k < nBeamChans; ++k){
437 : //blcin[3]=k;
438 : //trcin[3]=k;
439 : //Slicer slin(blcin, trcin, Slicer::endIsLast);
440 : //SubImage<Complex> subim(twoDPB, slin, true);
441 0 : TempImage<Complex> subim(IPosition(4, convSize_p, convSize_p, 1, 1), coordLastPlane, memtobeused/2.2);
442 0 : subim.set(Complex(1.0,0.0));
443 :
444 :
445 : /////Let's make a spheroid taper...
446 : /* Vector<Float> valsph(300);
447 : valsph.set(1.0);
448 : Sph_Conv<Float> d(1.0, 1.0);
449 : for(Int k=200; k <300; ++k){
450 : valsph(k)=d.value((Float)(k-200)/100.0);
451 : }
452 : Quantity fulrad(27,"arcmin");
453 : Quantity lefreq(1.4, "GHz");
454 : PBMath1DNumeric taper(valsph,fulrad,lefreq);
455 : taper.applyPB(subim, subim, fieldDir);
456 : */
457 : ///////////////
458 : /*{
459 : PagedImage<Complex> thisScreen(subim.shape(), subim.coordinates(), "ELTAPER");
460 : //LatticeExpr<Float> le(abs(twoDPB2));
461 : thisScreen.copyData(subim);
462 : }*/
463 : //twoDPB.putSlice(screen, start);
464 : /////////////
465 : //Double wtime0=omp_get_wtime();
466 : //////////
467 0 : sj_p->apply(subim, subim, vb, 0);
468 : ///////////////
469 : //Double wtime1=omp_get_wtime();
470 : /////////////////
471 : /*{
472 : PagedImage<Complex> thisScreen(subim.shape(), subim.coordinates(), "BEAM-ELTAPER");
473 : //LatticeExpr<Float> le(abs(twoDPB2));
474 : thisScreen.copyData(subim);
475 : }*/
476 : //LatticeFFT::cfft2d(subim);
477 : //////
478 :
479 :
480 0 : ft_p.c2cFFT(subim);
481 :
482 : // cerr << "make pb " << wtime1-wtime0 << " fft " << omp_get_wtime()-wtime1 << endl;
483 : // }
484 : //tim.show("after an apply" );
485 : //tim.mark();
486 0 : TempImage<Float> screen2(TiledShape(IPosition(4, convSize_p, convSize_p, 1, 1)), coordLastPlane, memtobeused/2.2);
487 0 : screen2.set(1.0);
488 0 : TempImage<Complex> subout(TiledShape(IPosition(4, convSize_p, convSize_p, 1, 1)), coordLastPlane, memtobeused/2.2);
489 :
490 : //////Making a reference on half of the lattice as on the Mac rcfft is failing for some
491 : //////reason
492 0 : SubImage<Complex> halfsubout(subout, Slicer(IPosition(4, 0, 0, 0, 0),
493 0 : IPosition(4, convSize_p/2, convSize_p-1, 0, 0),
494 0 : Slicer::endIsLast), true);
495 : /////////////////////
496 : //wtime0=omp_get_wtime();
497 : ///////////////////
498 0 : sj_p->applySquare(screen2, screen2, vb, 0);
499 : ///////////////
500 : //wtime1=omp_get_wtime();
501 : /////////////////
502 : /////////
503 : /* arr.resize();
504 : Array<Float> arrf;
505 : screen2.get(arrf, false);
506 : isRef=halfsubout.get (arr,false);
507 : Bool delf;
508 : Float* scrf=arrf.getStorage(delf);
509 : scr= arr.getStorage(del);
510 :
511 : ft_p.r2cFFT(scr, scrf, Long(convSize_p), Long(convSize_p));
512 : cerr << "isRef " << isRef << " iscopy " << del << " delf " << delf << endl;
513 : if(del)
514 : arr.putStorage(scr, del);
515 : if(delf)
516 : arrf.putStorage(scrf, delf);
517 : if(!isRef)
518 : halfsubout.put(arr.reform(IPosition(4, convSize_p/2+1, convSize_p, 1, 1)));
519 : */
520 : ////////////////////
521 0 : ft_p.r2cFFT(halfsubout, screen2);
522 : //LatticeFFT::rcfft(halfsubout, screen2, true, false);
523 : //cerr << "make pb2 " << wtime1-wtime0 << " fft " << omp_get_wtime()-wtime1 << endl;
524 : //Real FFT fills only first half of the array
525 : //making it look like a Complex to Complex FFT
526 0 : IPosition iblc(4, 0, 3*subout.shape()(1)/8, 0, 0);
527 0 : IPosition itrc(4, 0, 5*subout.shape()(1)/8, 0, 0);
528 0 : for(Int x=subout.shape()(0)/2; x <(5*subout.shape()(0)/8); ++x){
529 :
530 0 : iblc[0]=x-subout.shape()(0)/2;
531 0 : itrc[0]=x-subout.shape()(0)/2;
532 0 : Slicer isl(iblc, itrc, Slicer::endIsLast);
533 0 : iblc[0]=x;
534 0 : subout.putSlice(subout.getSlice(isl), iblc);
535 : }
536 0 : for(Int x=subout.shape()(0)/2+1; x <(5*subout.shape()(0)/8); ++x){
537 :
538 0 : iblc[0]=x;
539 0 : itrc[0]=x;
540 0 : Slicer isl(iblc, itrc, Slicer::endIsLast);
541 0 : iblc[0]=subout.shape()(0)-x;
542 0 : subout.putSlice(subout.getSlice(isl), iblc);
543 0 : if(x==(subout.shape()(0)-1)){
544 0 : iblc[0]=0;
545 0 : subout.putSlice(subout.getSlice(isl), iblc);
546 : }
547 : }
548 : //End of FFT's
549 : //tim.show("After apply2 ");
550 : /////////////////////
551 : //wtime0=omp_get_wtime();
552 : ///////////////////
553 0 : TempImage<Complex> twoDPB2(TiledShape(pbShape, IPosition(4, pbShape(0), pbShape(1), 1, 1)), coords, memtobeused/10.0);
554 :
555 0 : IPosition blcout(4, 0, 0, 0, nBeamChans-1);
556 0 : IPosition trcout(4, pbShape(0)-1, pbShape(1)-1, 0,nBeamChans-1);
557 0 : Slicer outsl(blcout, trcout, Slicer::endIsLast);
558 :
559 0 : IPosition blcin(4, convSize_p/2-pbShape(0)/2, convSize_p/2-pbShape(1)/2, 0, 0);
560 0 : IPosition trcin(4, convSize_p/2+pbShape(0)/2-1, convSize_p/2+pbShape(1)/2-1, 0, 0);
561 0 : Slicer insl(blcin, trcin, Slicer::endIsLast);
562 : {
563 0 : SubImage<Complex> subtwoDPB(twoDPB, outsl, true);
564 0 : SubImage<Complex> intwoDPB(subim, insl, false);
565 0 : subtwoDPB.copyData(intwoDPB);
566 : }
567 : {
568 0 : SubImage<Complex> subtwoDPB2(twoDPB2, outsl, true);
569 0 : SubImage<Complex> intwoDPB2(subout, insl, false);
570 0 : subtwoDPB2.copyData(intwoDPB2);
571 : }
572 :
573 0 : if(nBeamChans > 0){
574 0 : blcin=IPosition(4,0,0,0, nBeamChans-1);
575 0 : trcin=IPosition(4, pbShape(0)-1, pbShape(1)-1, 0, nBeamChans-1);
576 0 : Slicer slin(blcin, trcin, Slicer::endIsLast);
577 0 : SubImage<Complex> origPB(twoDPB, slin, false);
578 0 : IPosition elshape= origPB.shape();
579 0 : Matrix<Complex> i1=origPB.get(true);
580 0 : SubImage<Complex> origPB2(twoDPB2, slin, false);
581 0 : Matrix<Complex> i2=origPB2.get(true);
582 0 : Int cenX=i1.shape()(0)/2;
583 0 : Int cenY=i1.shape()(1)/2;
584 :
585 :
586 0 : for (Int kk=0; kk < nBeamChans; ++kk){
587 0 : Double fratio=beamFreqs(kk)/beamFreqs(nBeamChans-1);
588 : //cerr << "fratio " << fratio << endl;
589 0 : Float convRatio=convSamp/convSampling;
590 0 : blcin[3]=kk;
591 0 : trcin[3]=kk;
592 : //Slicer slout(blcin, trcin, Slicer::endIsLast);
593 0 : Matrix<Complex> o1(i1.shape(), Complex(0.0));
594 0 : Matrix<Complex> o2(i2.shape(), Complex(0.0));
595 0 : for (Int yy=0; yy < i1.shape()(1); ++yy){
596 : //Int nyy= (Double(yy-cenY)*fratio) + cenY;
597 0 : Double nyy= (Double((yy-cenY)*convRatio)/fratio) + cenY;
598 0 : Double cyy=ceil(nyy);
599 0 : Double fyy= floor(nyy);
600 0 : Int iy=nyy > fyy+0.5 ? Int(cyy) : Int(fyy);
601 0 : if(cyy <2*cenY && fyy >=0.0)
602 0 : for(Int xx=0; xx < i1.shape()(0); ++ xx){
603 : //Int nxx= Int(Double(xx-cenX)*fratio) + cenX;
604 0 : Double nxx= Int(Double((xx-cenX)*convRatio)/fratio) + cenX;
605 0 : Double cxx=ceil(nxx);
606 0 : Double fxx= floor(nxx);
607 0 : Int ix=nxx > fxx+0.5 ? Int(cxx) : Int(fxx) ;
608 0 : if(cxx < 2*cenX && fxx >=0.0 ){
609 : //Double dist=sqrt((nxx-cxx)*(nxx-cxx)+(nyy-cyy)*(nyy-cyy))/sqrt(2.0);
610 : //o1(xx, yy)=float(1-dist)*i1(fxx, fyy)+ dist*i1(cxx,cyy);
611 0 : o1(xx, yy)=i1( ix, iy);
612 : //o2(xx, yy)=i2(nxx, nyy);
613 : //o2(xx, yy)=float(1-dist)*i2(fxx, fyy)+ dist*i2(cxx,cyy);
614 0 : o2(xx, yy)=i2(ix, iy);
615 : }
616 : }
617 : }
618 0 : twoDPB.putSlice(o1.reform(elshape), blcin);
619 0 : twoDPB2.putSlice(o2.reform(elshape), blcin);
620 : }
621 :
622 : }
623 : //cerr << "make multfreq beam " << omp_get_wtime()-wtime0 << endl;
624 : /*
625 : {
626 : TempImage<Float> screen2(TiledShape(pbShape, IPosition(4, pbShape(0), pbShape(1), 1, 1)), coords, memtobeused);
627 : // Matrix<Float> screenoo(convSize_p, convSize_p);
628 : //screenoo.set(1.0);
629 : //screen2.putSlice(screenoo,start);
630 : //screen2.set(1.0);
631 : for (Int k=0; k < nBeamChans; ++k){
632 : blcin[3]=k;
633 : trcin[3]=k;
634 : Slicer slin(blcin, trcin, Slicer::endIsLast);
635 : SubImage<Float> subim(screen2, slin, true);
636 : SubImage<Complex> subout(twoDPB2, slin, true);
637 : subim.set(1.0);
638 : //twoDPB.putSlice(screen, start);
639 : sj_p->applySquare(subim, subim, vb, 0);
640 : //// LatticeExpr<Complex> le(subim);
641 : //// subout.copyData(le);
642 : ///// LatticeFFT::cfft2d(subout);
643 :
644 : LatticeFFT::rcfft(subout, subim, true, false);
645 : IPosition iblc(4, 0, 3*subout.shape()(1)/8, 0, 0);
646 : IPosition itrc(4, 0, 5*subout.shape()(1)/8, 0, 0);
647 : for(Int x=subout.shape()(0)/2; x <(5*subout.shape()(0)/8); ++x){
648 :
649 : iblc[0]=x-subout.shape()(0)/2;
650 : itrc[0]=x-subout.shape()(0)/2;
651 : Slicer isl(iblc, itrc, Slicer::endIsLast);
652 : iblc[0]=x;
653 : subout.putSlice(subout.getSlice(isl), iblc);
654 : }
655 : for(Int x=subout.shape()(0)/2+1; x <(5*subout.shape()(0)/8); ++x){
656 :
657 : iblc[0]=x;
658 : itrc[0]=x;
659 : Slicer isl(iblc, itrc, Slicer::endIsLast);
660 : iblc[0]=subout.shape()(0)-x;
661 : subout.putSlice(subout.getSlice(isl), iblc);
662 : if(x==(subout.shape()(0)-1)){
663 : iblc[0]=0;
664 : subout.putSlice(subout.getSlice(isl), iblc);
665 : }
666 : }
667 :
668 : }
669 :
670 : //sj_p->applySquare(screen2, screen2, vb, 0);
671 : //LatticeExpr<Complex> le(screen2);
672 : //twoDPB2.copyData(le);
673 : }
674 :
675 : */
676 :
677 : /*
678 : if(0) {
679 : CoordinateSystem ftCoords(coords);
680 : //directionIndex=ftCoords.findCoordinate(Coordinate::DIRECTION);
681 : //AlwaysAssert(directionIndex>=0, AipsError);
682 : dc=coords.directionCoordinate(directionIndex);
683 : Vector<Bool> axes(2); axes(0)=true;axes(1)=true;
684 : Vector<Int> shape(2); shape(0)=convSize_p;shape(1)=convSize_p;
685 : Coordinate* ftdc=dc.makeFourierCoordinate(axes,shape);
686 : //ftCoords.replaceCoordinate(*ftdc, directionIndex);
687 : delete ftdc; ftdc=0;
688 : ostringstream os1;
689 : os1 << "Screen_" << vb.fieldId() ;
690 : PagedImage<Complex> thisScreen(twoDPB2.shape(), ftCoords, String(os1));
691 : //LatticeExpr<Float> le(abs(twoDPB2));
692 : thisScreen.copyData(twoDPB2);
693 : LatticeFFT::cfft2d(thisScreen, false);
694 : PagedImage<Complex> thisScreen0(twoDPB.shape(), ftCoords, String("PB_")+String(os1));
695 : thisScreen0.copyData(twoDPB);
696 : LatticeFFT::cfft2d(thisScreen0, false);
697 : }
698 : */
699 : /*
700 : // Now FFT and get the result back
701 : //LatticeFFT::cfft2d(twoDPB);
702 : //LatticeFFT::cfft2d(twoDPB2);
703 :
704 : // Write out FT of screen as an image
705 : if(0) {
706 : CoordinateSystem ftCoords(coords);
707 : directionIndex=ftCoords.findCoordinate(Coordinate::DIRECTION);
708 : AlwaysAssert(directionIndex>=0, AipsError);
709 : dc=coords.directionCoordinate(directionIndex);
710 : Vector<Bool> axes(2); axes(0)=true;axes(1)=true;
711 : Vector<Int> shape(2); shape(0)=convSize_p;shape(1)=convSize_p;
712 : Coordinate* ftdc=dc.makeFourierCoordinate(axes,shape);
713 : ftCoords.replaceCoordinate(*ftdc, directionIndex);
714 : delete ftdc; ftdc=0;
715 : ostringstream os1;
716 : os1 << "FTScreen_" << vb.fieldId() ;
717 : PagedImage<Float> thisScreen(pbShape, ftCoords, String(os1));
718 : LatticeExpr<Float> le(abs(twoDPB2));
719 : thisScreen.copyData(le);
720 : }
721 : */
722 : //cerr << "twoDPB shape " << twoDPB.shape() << " slice shape " << IPosition(4, tempConvSize, tempConvSize, 1, 1) << endl;
723 0 : convFunc_p=twoDPB.getSlice(IPosition(4,0,0,0,0), IPosition(4, tempConvSize, tempConvSize, 1, 1), true);
724 0 : weightConvFunc_p.resize();
725 0 : weightConvFunc_p=twoDPB2.getSlice(IPosition(4,0,0,0,nBeamChans-1), IPosition(4, tempConvSize, tempConvSize, 1, 1), true);
726 : //convFunc/=max(abs(convFunc));
727 0 : Float maxAbsConvFunc=max(amplitude(weightConvFunc_p));
728 :
729 0 : Float minAbsConvFunc=min(amplitude(weightConvFunc_p));
730 : //cerr << "min max " << minAbsConvFunc << " " << maxAbsConvFunc << endl;
731 0 : convSupport_p=-1;
732 0 : Bool found=false;
733 : //Bool found2=true;
734 : //Int trial2=0;
735 0 : Int trial=0;
736 0 : for (trial=tempConvSize/2-2;trial>0;trial--) {
737 : //Searching down a diagonal
738 0 : if(abs(weightConvFunc_p(tempConvSize/2-trial, tempConvSize/2-trial)) > (1.0e-3*maxAbsConvFunc)) {
739 0 : found=true;
740 0 : trial=Int(sqrt(2.0*Float(trial*trial)));
741 0 : break;
742 : }
743 : }
744 0 : if(!found){
745 0 : if((maxAbsConvFunc-minAbsConvFunc) > (1.0e-3*maxAbsConvFunc))
746 0 : found=true;
747 : // if it drops by more than 2 magnitudes per pixel
748 0 : trial=( tempConvSize > (10*convSampling)) ? 5*convSampling : (tempConvSize/2 - 4*convSampling);
749 : }
750 :
751 0 : if(trial < 5*convSampling)
752 0 : trial=( tempConvSize > (10*convSampling)) ? 5*convSampling : (tempConvSize/2 - 4*convSampling);
753 :
754 0 : if(found) {
755 0 : convSupport_p=Int(0.5+Float(trial)/Float(convSampling))+1;
756 : }
757 : else {
758 : os << "Convolution function is misbehaved - support seems to be zero\n"
759 : << "Reasons can be: \n(1)The image definition not covering one or more of the pointings selected\n"
760 : << "(2) No unflagged data in a given pointing\n"
761 : << "(3) The entries in the POINTING subtable do not match the field being imaged."
762 : << "Please check, and try again with an empty POINTING subtable.)\n"
763 0 : << LogIO::EXCEPTION;
764 : }
765 :
766 : // Normalize such that plane 0 sums to 1 (when jumping in
767 : // steps of convSampling)
768 :
769 0 : Double pbSum=0.0;
770 : //Double pbSum2=0.0;
771 :
772 :
773 0 : for (Int iy=-convSupport_p;iy<=convSupport_p;iy++) {
774 0 : for (Int ix=-convSupport_p;ix<=convSupport_p;ix++) {
775 0 : Complex val=convFunc_p(ix*convSampling+tempConvSize/2,
776 0 : iy*convSampling+tempConvSize/2);
777 0 : pbSum+=real(val);
778 :
779 : //pbSum+=sqrt(real(val)*real(val)+ imag(val)*imag(val));
780 : }
781 : }
782 :
783 :
784 : //pbSum=sum(amplitude(convFunc_p))/Double(convSampling)/Double(convSampling);
785 :
786 0 : if(pbSum>0.0) {
787 0 : convFunc_p*=Complex(1.0/pbSum,0.0);
788 : }
789 : else {
790 : os << "Convolution function integral is not positive"
791 0 : << LogIO::EXCEPTION;
792 : }
793 :
794 : //##########################################
795 : os << "Convolution support = " << convSupport_p
796 : << " pixels in Fourier plane"
797 0 : << LogIO::POST;
798 :
799 0 : convSupportBlock_p.resize(actualConvIndex_p+1);
800 0 : convSizes_p.resize(actualConvIndex_p+1);
801 : //Only one beam for now...but later this should be able to
802 : // take all the beams for the different antennas.
803 0 : convSupportBlock_p[actualConvIndex_p]=new Vector<Int>(1);
804 0 : convSizes_p[actualConvIndex_p]=new Vector<Int> (1);
805 0 : (*(convSupportBlock_p[actualConvIndex_p]))[0]=convSupport_p;
806 0 : convFunctions_p.resize(actualConvIndex_p+1);
807 0 : convWeights_p.resize(actualConvIndex_p+1);
808 0 : convFunctions_p[actualConvIndex_p]= new Array<Complex>();
809 0 : convWeights_p[actualConvIndex_p]= new Array<Complex>();
810 0 : Int newConvSize=2*(convSupport_p+2)*convSampling;
811 : //NEED to chop this right ...and in the centre
812 0 : if(newConvSize >=tempConvSize)
813 0 : newConvSize=tempConvSize;
814 :
815 0 : IPosition blc(4, (tempConvSize/2)-(newConvSize/2),
816 0 : (tempConvSize/2)-(newConvSize/2), 0, 0);
817 0 : IPosition trc(4, (tempConvSize/2)+(newConvSize/2-1),
818 0 : (tempConvSize/2)+(newConvSize/2-1), 0, nBeamChans-1);
819 0 : convFunctions_p[actualConvIndex_p]->resize(IPosition(5, newConvSize, newConvSize, 1, nBeamChans,1));
820 : //cerr << "convFunc shape " << (convFunctions_p[actualConvIndex_p])->shape() <<
821 : //" " << " twoDPB shape " <<twoDPB.get(false)(blc,trc).shape() << endl;
822 0 : convFunctions_p[actualConvIndex_p]->copyMatchingPart(twoDPB.get(false)(blc,trc));//*Complex(1.0/pbSum,0.0));
823 0 : convSize_p=newConvSize;
824 0 : convWeights_p[actualConvIndex_p]->resize(IPosition(5, newConvSize, newConvSize, 1, nBeamChans,1));
825 0 : convWeights_p[actualConvIndex_p]->copyMatchingPart(twoDPB2.get(false)(blc,trc));//*Complex(1.0/pbSum2,0.0));
826 0 : blc.resize(5, false);
827 0 : trc.resize(5,false);
828 0 : blc=IPosition(5, 0, 0, 0, 0,0);
829 0 : trc=IPosition(5, newConvSize-1, newConvSize-1, 0, 0,0);
830 0 : for (Int bc=0; bc< nBeamChans; ++bc){
831 0 : blc[3]=bc;
832 0 : trc[3]=bc;
833 0 : pbSum=real(sum((*convFunctions_p[actualConvIndex_p])(blc,trc)))/Double(convSampling)/Double(convSampling);
834 0 : (*convFunctions_p[actualConvIndex_p])(blc,trc)= (*convFunctions_p[actualConvIndex_p])(blc,trc)/pbSum;
835 0 : (*convWeights_p[actualConvIndex_p])(blc,trc)= (*convWeights_p[actualConvIndex_p])(blc,trc)/pbSum;
836 : }
837 :
838 :
839 :
840 0 : convFunc_p.resize();//break any reference
841 0 : (*convSizes_p[actualConvIndex_p])[0]=convSize_p;
842 0 : doneMainConv_p[actualConvIndex_p]=true;
843 :
844 : }
845 : else{
846 0 : convSize_p=(*convSizes_p[actualConvIndex_p])[0];
847 :
848 : }
849 :
850 : //Apply the shift phase gradient
851 0 : convFunc.resize();
852 0 : weightConvFunc.resize();
853 0 : convFunc.assign(*(convFunctions_p[actualConvIndex_p]));
854 0 : weightConvFunc.assign(*(convWeights_p[actualConvIndex_p]));
855 : Bool copyconv, copywgt;
856 0 : Complex *cv=convFunc.getStorage(copyconv);
857 0 : Complex *wcv=weightConvFunc.getStorage(copywgt);
858 : //cerr << "Field " << vb.fieldId() << " spw " << vb.spectralWindow() << " phase grad: " << pixFieldDir << endl;
859 :
860 0 : for (Int nc=0; nc < nBeamChans; ++nc){
861 0 : Int planeoffset=nc*convSize_p*convSize_p;
862 0 : for (Int iy=0;iy<convSize_p;iy++) {
863 : Double cy, sy;
864 : Int offset;
865 :
866 0 : SINCOS(Double(iy-convSize_p/2)*pixFieldDir(1), sy, cy);
867 0 : Complex phy(cy,sy) ;
868 0 : offset = iy*convSize_p+planeoffset;
869 0 : for (Int ix=0;ix<convSize_p;ix++) {
870 : Double cx, sx;
871 0 : SINCOS(Double(ix-convSize_p/2)*pixFieldDir(0), sx, cx);
872 0 : Complex phx(cx,sx) ;
873 0 : cv[ix+offset]= cv[ix+offset]*phx*phy;
874 0 : wcv[ix+offset]= wcv[ix+offset]*phx*phy;
875 :
876 : }
877 : }
878 : }
879 0 : convFunc.putStorage(cv, copyconv);
880 0 : weightConvFunc.putStorage(wcv, copywgt);
881 0 : convsize.resize();
882 0 : convsize=*(convSizes_p[actualConvIndex_p]);
883 0 : convSupport.resize();
884 0 : convSupport=(*(convSupportBlock_p[actualConvIndex_p]));
885 :
886 :
887 0 : }
888 :
889 0 : void SimplePBConvFunc::setSkyJones(SkyJones* sj){
890 0 : sj_p=sj;
891 0 : }
892 :
893 0 : void SimplePBConvFunc::findUsefulChannels(Vector<Int>& chanMap, Vector<Double>& chanFreqs, const VisBuffer& vb, const Vector<Double>& freq){
894 0 : chanMap.resize(freq.nelements());
895 0 : Vector<Double> localfreq=vb.frequency();
896 0 : Double minfreq=min(freq);
897 :
898 0 : Double origwidth=freq.nelements()==1 ? 1e12 : (max(freq)-min(freq))/(freq.nelements()-1);
899 : ///Fractional bandwidth which will trigger mutiple PB in one spw
900 0 : Double tol=(max(freq))*0.5/100;
901 0 : Int nchan=Int(lround((max(freq)-min(freq))/tol));
902 :
903 : //cerr << "TOLERA " << tol << " nchan " << nchan << " vb.nchan " << vb.nChannel() << endl;
904 : //Number of beams that matters are the ones in the data
905 0 : if(nchan > vb.nChannel())
906 0 : nchan=vb.nChannel();
907 :
908 0 : if(tol < origwidth) tol=origwidth;
909 0 : chanFreqs.resize();
910 0 : if(nchan >= (Int)(freq.nelements()-1)) { indgen(chanMap); chanFreqs=freq; return;}
911 0 : if((nchan==0) || (freq.nelements()==1)) { chanFreqs=Vector<Double>(1, freq[0]);chanMap.set(0); return;}
912 :
913 : //readjust the tolerance...
914 0 : tol=(max(freq)-min(freq)+origwidth)/Double(nchan);
915 0 : chanFreqs.resize(nchan);
916 0 : for (Int k=0; k < nchan; ++k)
917 0 : chanFreqs[k]=minfreq-origwidth+tol/2.0+tol*Double(k);
918 0 : Int activechan=0;
919 0 : chanMap.set(-1);
920 0 : for (uInt k=0; k < chanMap.nelements(); ++k){
921 :
922 0 : while((activechan< nchan) && Float(fabs(freq[k]-chanFreqs[activechan])) > Float(tol/2.0)){
923 : // cerr << "k " << k << " atcivechan " << activechan << " comparison "
924 : // << freq[k] << " " << chanFreqs[activechan] << endl;
925 0 : ++activechan;
926 : }
927 0 : if(activechan != nchan)
928 0 : chanMap[k]=activechan;
929 : //////////////////
930 : //if(chanMap[k] < 0)
931 : //cerr << "freq diffs " << freq[k]-chanFreqs << " TOL " << tol/2.0 << endl;
932 :
933 : ///////////////////////////
934 0 : activechan=0;
935 : }
936 :
937 0 : return;
938 : }
939 :
940 :
941 0 : Bool SimplePBConvFunc::checkPBOfField(const VisBuffer& vb){
942 : //Int fieldid=vb.fieldId();
943 0 : String msid=vb.msName(true);
944 : /*
945 : if(convFunctionMap_p.ndefined() > 0){
946 : if (((fluxScale_p.shape()[3] != nchan_p) || (fluxScale_p.shape()[2] != npol_p)) && calcFluxScale_p){
947 : convFunctionMap_p.clear();
948 : }
949 : }
950 : // if you rename the ms might be a problem
951 : String mapid=msid+String("_")+String::toString(fieldid);
952 : if(convFunctionMap_p.ndefined() == 0){
953 : convFunctionMap_p.define(mapid, 0);
954 : actualConvIndex_p=0;
955 : if(calcFluxScale_p){
956 : // 0ne channel only is needed to keep track of pb coverage
957 : if(fluxScale_p.shape().nelements()==0){
958 : fluxScale_p=TempImage<Float>(IPosition(4,nx_p,ny_p,npol_p,1), csys_p);
959 : fluxScale_p.set(0.0);
960 : }
961 : filledFluxScale_p=false;
962 : }
963 : return false;
964 : }
965 :
966 : if(!convFunctionMap_p.isDefined(mapid)){
967 : actualConvIndex_p=convFunctionMap_p.ndefined();
968 : convFunctionMap_p.define(mapid, actualConvIndex_p);
969 : return false;
970 : }
971 : else{
972 : actualConvIndex_p=convFunctionMap_p(mapid);
973 : convFunc_p.resize(); // break any reference
974 : weightConvFunc_p.resize();
975 : //Here we will need to use the right xyPlane for different PA range.
976 : //convFunc_p.reference(convFunctions_p[actualConvIndex_p]->xyPlane(0));
977 : //weightConvFunc_p.reference(convWeights_p[actualConvIndex_p]->xyPlane(0));
978 : //Again this for one time of antenna only later should be fixed for all
979 : // antennas independently
980 : convSupport_p=(*convSupportBlock_p[actualConvIndex_p])[0];
981 : convSize_p=(*convSizes_p[actualConvIndex_p])[0];
982 :
983 : }
984 : */
985 :
986 0 : return true;
987 :
988 :
989 :
990 : }
991 0 : ImageInterface<Float>& SimplePBConvFunc::getFluxScaleImage(){
992 :
993 0 : if(!calcFluxScale_p)
994 0 : throw(AipsError("Programmer error: Cannot get flux scale"));
995 0 : if(!filledFluxScale_p){
996 0 : IPosition blc=fluxScale_p.shape();
997 0 : IPosition trc=fluxScale_p.shape();
998 0 : blc(0)=0; blc(1)=0; trc(0)=nx_p-1; trc(1)=ny_p-1;
999 :
1000 0 : for (Int j=0; j < fluxScale_p.shape()(2); ++j){
1001 0 : for (Int k=0; k < fluxScale_p.shape()(3) ; ++k){
1002 :
1003 0 : blc(2)=j; trc(2)=j;
1004 0 : blc(3)=k; trc(3)=k;
1005 0 : Slicer sl(blc, trc, Slicer::endIsLast);
1006 0 : SubImage<Float> fscalesub(fluxScale_p, sl, true);
1007 : Float planeMax;
1008 0 : LatticeExprNode LEN = max( fscalesub );
1009 0 : planeMax = LEN.getFloat();
1010 0 : if(planeMax !=0){
1011 0 : fscalesub.copyData( (LatticeExpr<Float>) (fscalesub/planeMax));
1012 :
1013 : }
1014 : }
1015 : }
1016 :
1017 0 : filledFluxScale_p=true;
1018 : }
1019 :
1020 :
1021 0 : return fluxScale_p;
1022 :
1023 : }
1024 :
1025 :
1026 0 : Bool SimplePBConvFunc::toRecord(RecordInterface& rec){
1027 0 : Int numConv=convFunctions_p.nelements();
1028 : // not saving the protected variables as they are generated by
1029 : // the first call to storeImageParams
1030 : try{
1031 0 : rec.define("name", "SimplePBConvFunc");
1032 0 : rec.define("numconv", numConv);
1033 : //cerr << "num of conv " << numConv << " " << convFunctionMap_p.ndefined() << " " <<convFunctions_p.nelements() << endl;
1034 0 : std::map<String, Int>::iterator it=vbConvIndex_p.begin();
1035 0 : for (Int k=0; k < numConv; ++k){
1036 0 : rec.define("convfunctions"+String::toString(k), *(convFunctions_p[k]));
1037 0 : rec.define("convweights"+String::toString(k), *(convWeights_p[k]));
1038 0 : rec.define("convsizes"+String::toString(k), *(convSizes_p[k]));
1039 0 : rec.define("convsupportblock"+String::toString(k), *(convSupportBlock_p[k]));
1040 : //cerr << "k " << k << " key " << convFunctionMap_p.getKey(k) << " val " << convFunctionMap_p.getVal(k) << endl;
1041 0 : rec.define(String("key")+String::toString(k), it->first);
1042 0 : rec.define(String("val")+String::toString(k), it->second);
1043 0 : it++;
1044 : }
1045 0 : rec.define("pbclass", Int(pbClass_p));
1046 0 : rec.define("actualconvindex", actualConvIndex_p);
1047 0 : rec.define("donemainconv", doneMainConv_p);
1048 : //The following is not needed ..can be regenerated
1049 : //rec.define("pointingpix", pointingPix_p);
1050 : }
1051 0 : catch(AipsError &x) {
1052 0 : return false;
1053 : }
1054 0 : return true;
1055 : }
1056 :
1057 0 : Bool SimplePBConvFunc::fromRecord(String& err, const RecordInterface& rec, Bool calcFluxneeded){
1058 0 : Int numConv=0;
1059 : //make sure storeImageParams is triggered
1060 0 : nchan_p=0;
1061 :
1062 : try{
1063 0 : if(!rec.isDefined("name") || rec.asString("name") != "SimplePBConvFunc"){
1064 0 : throw(AipsError("Wrong record to recover SimplePBConvFunc from"));
1065 : }
1066 0 : rec.get("numconv", numConv);
1067 0 : convFunctions_p.resize(numConv, true, false);
1068 0 : convSupportBlock_p.resize(numConv, true, false);
1069 0 : convWeights_p.resize(numConv, true, false);
1070 0 : convSizes_p.resize(numConv, true, false);
1071 0 : convFunctionMap_p.clear( );
1072 0 : vbConvIndex_p.erase(vbConvIndex_p.begin(), vbConvIndex_p.end());
1073 0 : for (Int k=0; k < numConv; ++k){
1074 0 : convFunctions_p[k]=new Array<Complex>();
1075 0 : convWeights_p[k]=new Array<Complex>();
1076 0 : convSizes_p[k]=new Vector<Int>();
1077 0 : convSupportBlock_p[k]=new Vector<Int>();
1078 0 : rec.get("convfunctions"+String::toString(k), *(convFunctions_p[k]));
1079 0 : rec.get("convsupportblock"+String::toString(k), *(convSupportBlock_p[k]));
1080 0 : rec.get("convweights"+String::toString(k), *(convWeights_p[k]));
1081 0 : rec.get("convsizes"+String::toString(k), *(convSizes_p[k]));
1082 0 : String key;
1083 : Int val;
1084 0 : rec.get(String("key")+String::toString(k), key);
1085 0 : rec.get(String("val")+String::toString(k), val);
1086 0 : vbConvIndex_p[key]=val;
1087 0 : ant1PointVal_p.clear();
1088 0 : ant1PointingCache_p.resize();
1089 : //convFunctionMap_p.define(key,val);
1090 : }
1091 0 : pbClass_p=static_cast<PBMathInterface::PBClass>(rec.asInt("pbclass"));
1092 0 : rec.get("actualconvindex", actualConvIndex_p);
1093 0 : pointingPix_p.resize();
1094 : //rec.get("pointingpix", pointingPix_p);
1095 0 : calcFluxScale_p=calcFluxneeded;
1096 :
1097 : }
1098 0 : catch(AipsError & x) {
1099 0 : err=x.getMesg();
1100 0 : return false;
1101 : }
1102 0 : return true;
1103 :
1104 : }
1105 0 : void SimplePBConvFunc::addPBToFlux(const VisBuffer& vb){
1106 :
1107 0 : if(calcFluxScale_p){
1108 0 : if(fluxScale_p.shape().nelements()==0)
1109 : { //cerr << "nx_p " << nx_p << endl;
1110 0 : calcFluxScale_p=False;
1111 0 : return;
1112 : }
1113 0 : Vector<Int> pixdepoint(2, -100000);
1114 0 : convertArray(pixdepoint, thePix_p);
1115 0 : if((pixdepoint(0) >=0) && (pixdepoint(0) < pointingPix_p.shape()[0]) && (pixdepoint(1) >=0) && (pixdepoint(1) < pointingPix_p.shape()[1]) && !pointingPix_p(pixdepoint(0), pixdepoint(1))){
1116 0 : TempImage<Float> thispb(fluxScale_p.shape(), fluxScale_p.coordinates());
1117 0 : thispb.set(1.0);
1118 0 : sj_p->applySquare(thispb, thispb, vb, 0);
1119 0 : LatticeExpr<Float> le(fluxScale_p+thispb);
1120 0 : fluxScale_p.copyData(le);
1121 0 : pointingPix_p(pixdepoint(0), pixdepoint(1))=true;
1122 : //LatticeExprNode LEN = max(fluxScale_p);
1123 : //Float maxsca=LEN.getFloat();
1124 : //Tempporary fix when cubesky is chunking...do not add on
1125 : //already defined position
1126 : //if(maxsca > 1.98){
1127 : // cerr << "avoiding subtract " << endl;
1128 : //fluxScale_p.copyData(LatticeExpr<Float>(fluxScale_p-thispb));
1129 :
1130 : //}
1131 : /*
1132 : if(0) {
1133 : ostringstream os1;
1134 : os1 << "SINGLE_" << vb.fieldId() ;
1135 : PagedImage<Float> thisScreen(fluxScale_p.shape(), fluxScale_p.coordinates(), String(os1));
1136 : thisScreen.copyData(thispb);
1137 : ostringstream os2;
1138 : os2 << "ALL_" << vb.fieldId() ;
1139 : PagedImage<Float> thisScreen2(fluxScale_p.shape(), fluxScale_p.coordinates(), String(os2));
1140 : thisScreen2.copyData(fluxScale_p);
1141 : }
1142 : */
1143 : }
1144 : }
1145 :
1146 : }
1147 :
1148 0 : void SimplePBConvFunc::sliceFluxScale(Int npol) {
1149 0 : IPosition fshp=fluxScale_p.shape();
1150 0 : if (fshp(2)>npol){
1151 0 : npol_p=npol;
1152 : // use first npol planes...
1153 0 : IPosition blc(4,0,0,0,0);
1154 0 : IPosition trc(4,fluxScale_p.shape()(0)-1, fluxScale_p.shape()(1)-1,npol-1,fluxScale_p.shape()(3)-1);
1155 0 : Slicer sl=Slicer(blc, trc, Slicer::endIsLast);
1156 : //writeable if possible
1157 0 : SubImage<Float> fluxScaleSub = SubImage<Float> (fluxScale_p, sl, true);
1158 0 : fluxScale_p = TempImage<Float>(fluxScaleSub.shape(),fluxScaleSub.coordinates());
1159 0 : LatticeExpr<Float> le(fluxScaleSub);
1160 0 : fluxScale_p.copyData(le);
1161 : }
1162 0 : }
1163 :
1164 : } //# NAMESPACE CASA - END
1165 :
1166 :
1167 :
1168 :
1169 :
1170 :
1171 :
|
