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