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