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 776 : SimplePBConvFunc::SimplePBConvFunc(): nchan_p(-1),
81 776 : npol_p(-1), pointToPix_p(), directionIndex_p(-1), thePix_p(0),
82 388 : filledFluxScale_p(false),doneMainConv_p(0),
83 :
84 1164 : calcFluxScale_p(true), usePointingTable_p(False), actualConvIndex_p(-1), convSize_p(0), convSupport_p(0), pointingPix_p() {
85 : //
86 :
87 388 : pbClass_p=PBMathInterface::COMMONPB;
88 388 : ft_p=FFT2D(true);
89 388 : usePointingTable_p=False;
90 388 : }
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 595 : SimplePBConvFunc::~SimplePBConvFunc(){
112 : //
113 :
114 595 : }
115 :
116 85077 : void SimplePBConvFunc::storeImageParams(const ImageInterface<Complex>& iimage,
117 : const vi::VisBuffer2& vb){
118 : //image signature changed...rather simplistic for now
119 85077 : if((iimage.shape().product() != nx_p*ny_p*nchan_p*npol_p) || nchan_p < 1){
120 162 : csys_p=iimage.coordinates();
121 162 : Int coordIndex=csys_p.findCoordinate(Coordinate::DIRECTION);
122 162 : AlwaysAssert(coordIndex>=0, AipsError);
123 162 : directionIndex_p=coordIndex;
124 162 : dc_p=csys_p.directionCoordinate(directionIndex_p);
125 162 : ObsInfo imInfo=csys_p.obsInfo();
126 162 : String tel= imInfo.telescope();
127 162 : MPosition pos;
128 162 : MSColumns mscol(vb.ms());
129 162 : if (vb.subtableColumns().observation().nrow() > 0) {
130 162 : tel =vb.subtableColumns().observation().telescopeName()(mscol.observationId()(0));
131 : }
132 290 : if (tel.length() == 0 || !tel.contains("VLA") ||
133 128 : !MeasTable::Observatory(pos,tel)) {
134 : // unknown observatory, use first antenna
135 34 : Int ant1=vb.antenna1()(0);
136 34 : pos=vb.subtableColumns().antenna().positionMeas()(ant1);
137 : }
138 162 : imInfo.setTelescope(tel);
139 162 : csys_p.setObsInfo(imInfo);
140 : //Store this to build epochs via the time access of visbuffer later
141 162 : timeMType_p=MEpoch::castType(mscol.timeMeas()(0).getRef().getType());
142 162 : 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 162 : pointFrame_p=MeasFrame(imInfo.obsDate(), pos);
146 162 : MDirection::Ref elRef(dc_p.directionType(), pointFrame_p);
147 : //For now we set the conversion from this direction
148 162 : pointToPix_p=MDirection::Convert( MDirection(), elRef);
149 162 : nx_p=iimage.shape()(coordIndex);
150 162 : ny_p=iimage.shape()(coordIndex+1);
151 162 : pointingPix_p.resize(nx_p, ny_p);
152 162 : pointingPix_p.set(false);
153 162 : coordIndex=csys_p.findCoordinate(Coordinate::SPECTRAL);
154 162 : Int pixAxis=csys_p.pixelAxes(coordIndex)[0];
155 162 : nchan_p=iimage.shape()(pixAxis);
156 162 : coordIndex=csys_p.findCoordinate(Coordinate::STOKES);
157 162 : pixAxis=csys_p.pixelAxes(coordIndex)[0];
158 162 : npol_p=iimage.shape()(pixAxis);
159 162 : if(calcFluxScale_p){
160 162 : if(fluxScale_p.shape().nelements()==0){
161 162 : fluxScale_p=TempImage<Float>(IPosition(4,nx_p,ny_p,npol_p,nchan_p), csys_p);
162 162 : fluxScale_p.set(0.0);
163 : }
164 162 : filledFluxScale_p=false;
165 : }
166 :
167 162 : }
168 :
169 85077 : }
170 :
171 85077 : void SimplePBConvFunc::toPix(const vi::VisBuffer2& vb, const MVDirection& extraShift, const Bool useExtraShift){
172 85077 : thePix_p.resize(2);
173 :
174 85077 : const MDirection& p1=pointingDirAnt1(vb);
175 85077 : 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 85077 : direction1_p=p1;
209 :
210 : //direction2_p=vb.direction2()(0);
211 : //For now
212 85077 : direction2_p=direction1_p;
213 :
214 : }
215 : //Should return both antennas direction in the future
216 :
217 85077 : if(useExtraShift){
218 1120 : direction1_p.shift(extraShift, False);
219 1120 : direction2_p.shift(extraShift, False);
220 : }
221 85077 : dc_p.toPixel(thePix_p, direction1_p);
222 :
223 85077 : }
224 :
225 204 : void SimplePBConvFunc::setWeightImage(CountedPtr<TempImage<Float> >& wgtimage){
226 204 : convWeightImage_p=wgtimage;
227 204 : filledFluxScale_p=false;
228 204 : calcFluxScale_p=true;
229 :
230 204 : }
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 85077 : Int SimplePBConvFunc::convIndex(const vi::VisBuffer2& vb, const uInt nchan){
246 170154 : String elkey=String::toString(vb.msId())+String("_")+String::toString(vb.spectralWindows()[0])+String("_")+String::toString(nchan);
247 85077 : if(vbConvIndex_p.count(elkey) > 0){
248 84790 : return vbConvIndex_p[elkey];
249 : }
250 287 : Int val=vbConvIndex_p.size();
251 287 : vbConvIndex_p[elkey]=val;
252 287 : return val;
253 85077 : }
254 :
255 85077 : const MDirection& SimplePBConvFunc::pointingDirAnt1(const vi::VisBuffer2& vb){
256 :
257 :
258 85077 : std::ostringstream oss;
259 :
260 85077 : oss << vb.msId() << "_" << vb.antenna1()(0) << "_";
261 85077 : oss.precision(13);
262 85077 : oss << vb.time()(0);
263 85077 : 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 85077 : if(ant1PointVal_p.count(elkey) > 0){
269 55686 : return ant1PointingCache_p[ant1PointVal_p[elkey]];
270 :
271 : }
272 29391 : Bool hasValidPointing=False;
273 29391 : if(Table::isReadable(vb.ms().pointingTableName())){
274 29391 : hasValidPointing=usePointingTable_p && (vb.ms().pointing().nrow() >0);
275 : }
276 :
277 29391 : Int val=ant1PointingCache_p.nelements();
278 29391 : ant1PointingCache_p.resize(val+1, true);
279 29391 : 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 29391 : ant1PointingCache_p[val]=vbutil_p->getPhaseCenter(vb);
285 : //ant1PointingCache_p[val]=vbUtil_p.getPointingDir(vb, vb.antenna1()(0), 0);
286 29391 : ant1PointVal_p[elkey]=val;
287 29391 : return ant1PointingCache_p[val];
288 :
289 85077 : }
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 13448 : 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 13448 : Int convSamp=2*convSampling;
346 : /////////////////////////
347 13448 : LogIO os1;
348 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 ;
349 : //////////////////////
350 13448 : storeImageParams(iimage, vb);
351 13448 : convFuncChanMap.resize(vb.nChannels());
352 13448 : Vector<Double> beamFreqs;
353 13448 : findUsefulChannels(convFuncChanMap, beamFreqs, vb, visFreq);
354 : //cerr << "CHANMAP " << convFuncChanMap << " beamFreqs " << beamFreqs << endl;
355 13448 : Int nBeamChans=beamFreqs.nelements();
356 : //indgen(convFuncChanMap);
357 13448 : convFuncPolMap.resize(vb.nCorrelations());
358 13448 : convFuncPolMap.set(0);
359 : //Only one plane in this version
360 13448 : convFuncRowMap.resize();
361 13448 : convFuncRowMap=Vector<Int>(vb.nRows(),0);
362 : //break reference
363 13448 : convFunc.resize();
364 13448 : weightConvFunc.resize();
365 13448 : LogIO os;
366 13448 : os << LogOrigin("SimplePBConv", "findConvFunction") << LogIO::NORMAL;
367 :
368 :
369 : // Get the coordinate system
370 13448 : CoordinateSystem coords(iimage.coordinates());
371 :
372 :
373 13448 : 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 13448 : Int directionIndex=directionIndex_p;
379 13448 : AlwaysAssert(directionIndex>=0, AipsError);
380 :
381 : // Set up the convolution function.
382 13448 : Int nx=nx_p;
383 13448 : Int ny=ny_p;
384 : // convSize_p=max(nx,ny)*convSampling;
385 : //cerr << "size " << nx << " " << ny << endl;
386 : //3 times the support size
387 13448 : if(doneMainConv_p.shape()[0] < (actualConvIndex_p+1)){
388 : // cerr << "resizing DONEMAIN " << doneMainConv_p.shape()[0] << endl;
389 20 : doneMainConv_p.resize(actualConvIndex_p+1, true);
390 20 : doneMainConv_p[actualConvIndex_p]=false;
391 : }
392 :
393 13448 : if(!(doneMainConv_p[actualConvIndex_p])){
394 :
395 : //convSize_p=4*(sj_p->support(vb, coords));
396 20 : 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 20 : CompositeNumber cn(uInt(convSize_p*2.0));
400 :
401 20 : convSize_p = cn.nextLargerEven(Int(convSize_p));
402 : //cerr << "convSize : " << convSize_p << endl;
403 :
404 20 : }
405 :
406 :
407 13448 : toPix(vb, extraShift, useExtraShift);
408 : //Timer tim;
409 : //tim.mark();
410 13448 : addPBToFlux(vb);
411 : //tim.show("After addPBToFlux");
412 13448 : DirectionCoordinate dc=dc_p;
413 :
414 : //where in the image in pixels is this pointing
415 13448 : Vector<Double> pixFieldDir(2);
416 13448 : pixFieldDir=thePix_p;
417 :
418 : //cerr << "pix of pointing " << pixFieldDir << endl;
419 13448 : MDirection fieldDir=direction1_p;
420 : //shift from center
421 13448 : pixFieldDir(0) = pixFieldDir(0) - Double(nx / 2);
422 13448 : pixFieldDir(1) = pixFieldDir(1) - Double(ny / 2);
423 :
424 : //phase gradient per pixel to apply
425 13448 : pixFieldDir(0) = -pixFieldDir(0)*2.0*C::pi/Double(nx)/Double(convSampling);
426 13448 : pixFieldDir(1) = -pixFieldDir(1)*2.0*C::pi/Double(ny)/Double(convSampling);
427 :
428 : //cerr << "DonemainConv " << doneMainConv_p[actualConvIndex_p] << endl;
429 13448 : if(!doneMainConv_p[actualConvIndex_p]){
430 20 : Vector<Double> sampling;
431 20 : sampling = dc.increment();
432 20 : sampling*=Double(convSamp);
433 20 : sampling(0)*=Double(nx)/Double(convSize_p);
434 20 : sampling(1)*=Double(ny)/Double(convSize_p);
435 20 : dc.setIncrement(sampling);
436 :
437 :
438 20 : Vector<Double> unitVec(2);
439 20 : unitVec=convSize_p/2;
440 20 : dc.setReferencePixel(unitVec);
441 :
442 :
443 : //make sure we are using the same units
444 20 : fieldDir.set(dc.worldAxisUnits()(0));
445 :
446 20 : dc.setReferenceValue(fieldDir.getAngle().getValue());
447 :
448 20 : coords.replaceCoordinate(dc, directionIndex);
449 20 : Int spind=coords.findCoordinate(Coordinate::SPECTRAL);
450 20 : SpectralCoordinate spCoord=coords.spectralCoordinate(spind);
451 20 : spCoord.setReferencePixel(Vector<Double>(1,0.0));
452 20 : spCoord.setReferenceValue(Vector<Double>(1, beamFreqs(0)));
453 20 : if(beamFreqs.nelements() >1)
454 20 : spCoord.setIncrement(Vector<Double>(1, beamFreqs(1)-beamFreqs(0)));
455 20 : coords.replaceCoordinate(spCoord, spind);
456 :
457 :
458 20 : CoordinateSystem coordLastPlane= coords;
459 20 : spCoord.setReferenceValue(Vector<Double>(1, beamFreqs(nBeamChans-1)));
460 20 : coordLastPlane.replaceCoordinate(spCoord, spind);
461 : //cerr << "BEAM freqs " << beamFreqs << endl;
462 :
463 : // coords.list(logIO(), MDoppler::RADIO, IPosition(), IPosition());
464 :
465 20 : Int tempConvSize=((convSize_p/4/(convSamp/convSampling))/2)*2;
466 20 : IPosition pbShape(4, tempConvSize, tempConvSize, 1, nBeamChans);
467 :
468 20 : Long memtot=HostInfo::memoryFree();
469 20 : 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 20 : if(memtot <= 2000000)
476 0 : memtobeused=0.0;
477 : //cerr << "mem to be used " << memtobeused << endl;
478 : //tim.mark();
479 20 : IPosition start(4, 0, 0, 0, 0);
480 : //IPosition pbSlice(4, convSize_p, convSize_p, 1, 1);
481 : //cerr << "pbshape " << pbShape << endl;
482 40 : 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 20 : convFunc_p.resize(tempConvSize, tempConvSize);
486 20 : 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 40 : TempImage<Complex> subim(IPosition(4, convSize_p, convSize_p, 1, 1), coordLastPlane, memtobeused/2.2);
502 20 : subim.set(Complex(1.0,0.0));
503 : //twoDPB.putSlice(screen, start);
504 20 : sj_p->apply(subim, subim, vb, 0);
505 : //LatticeFFT::cfft2d(subim);
506 20 : ft_p.c2cFFT(subim);
507 : // }
508 : //tim.show("after an apply" );
509 : //tim.mark();
510 40 : TempImage<Float> screen2(TiledShape(IPosition(4, convSize_p, convSize_p, 1, 1)), coordLastPlane, memtobeused/2.2);
511 20 : screen2.set(1.0);
512 40 : 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 40 : 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 20 : sj_p->applySquare(screen2, screen2, vb, 0);
517 20 : 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 20 : IPosition iblc(4, 0, 3*subout.shape()(1)/8, 0, 0);
522 20 : IPosition itrc(4, 0, 5*subout.shape()(1)/8, 0, 0);
523 16332 : for(Int x=subout.shape()(0)/2; x <(5*subout.shape()(0)/8); ++x){
524 :
525 16312 : iblc[0]=x-subout.shape()(0)/2;
526 16312 : itrc[0]=x-subout.shape()(0)/2;
527 16312 : Slicer isl(iblc, itrc, Slicer::endIsLast);
528 16312 : iblc[0]=x;
529 16312 : subout.putSlice(subout.getSlice(isl), iblc);
530 16312 : }
531 16312 : for(Int x=subout.shape()(0)/2+1; x <(5*subout.shape()(0)/8); ++x){
532 :
533 16292 : iblc[0]=x;
534 16292 : itrc[0]=x;
535 16292 : Slicer isl(iblc, itrc, Slicer::endIsLast);
536 16292 : iblc[0]=subout.shape()(0)-x;
537 16292 : subout.putSlice(subout.getSlice(isl), iblc);
538 16292 : if(x==(subout.shape()(0)-1)){
539 0 : iblc[0]=0;
540 0 : subout.putSlice(subout.getSlice(isl), iblc);
541 : }
542 16292 : }
543 : //End of FFT's
544 : //tim.show("After apply2 ");
545 40 : TempImage<Complex> twoDPB2(TiledShape(pbShape, IPosition(4, pbShape(0), pbShape(1), 1, 1)), coords, memtobeused/10.0);
546 :
547 20 : IPosition blcout(4, 0, 0, 0, nBeamChans-1);
548 20 : IPosition trcout(4, pbShape(0)-1, pbShape(1)-1, 0,nBeamChans-1);
549 20 : Slicer outsl(blcout, trcout, Slicer::endIsLast);
550 :
551 20 : IPosition blcin(4, convSize_p/2-pbShape(0)/2, convSize_p/2-pbShape(1)/2, 0, 0);
552 20 : IPosition trcin(4, convSize_p/2+pbShape(0)/2-1, convSize_p/2+pbShape(1)/2-1, 0, 0);
553 20 : Slicer insl(blcin, trcin, Slicer::endIsLast);
554 : {
555 20 : SubImage<Complex> subtwoDPB(twoDPB, outsl, true);
556 20 : SubImage<Complex> intwoDPB(subim, insl, false);
557 20 : subtwoDPB.copyData(intwoDPB);
558 20 : }
559 : {
560 20 : SubImage<Complex> subtwoDPB2(twoDPB2, outsl, true);
561 20 : SubImage<Complex> intwoDPB2(subout, insl, false);
562 20 : subtwoDPB2.copyData(intwoDPB2);
563 20 : }
564 :
565 20 : if(nBeamChans > 0){
566 20 : blcin=IPosition(4,0,0,0, nBeamChans-1);
567 20 : trcin=IPosition(4, pbShape(0)-1, pbShape(1)-1, 0, nBeamChans-1);
568 20 : Slicer slin(blcin, trcin, Slicer::endIsLast);
569 20 : SubImage<Complex> origPB(twoDPB, slin, false);
570 20 : IPosition elshape= origPB.shape();
571 20 : Matrix<Complex> i1=origPB.get(true);
572 20 : SubImage<Complex> origPB2(twoDPB2, slin, false);
573 20 : Matrix<Complex> i2=origPB2.get(true);
574 20 : Int cenX=i1.shape()(0)/2;
575 20 : Int cenY=i1.shape()(1)/2;
576 :
577 :
578 68 : for (Int kk=0; kk < nBeamChans; ++kk){
579 48 : Double fratio=beamFreqs(kk)/beamFreqs(nBeamChans-1);
580 : //cerr << "fratio " << fratio << endl;
581 48 : Float convRatio=convSamp/convSampling;
582 48 : blcin[3]=kk;
583 48 : trcin[3]=kk;
584 : //Slicer slout(blcin, trcin, Slicer::endIsLast);
585 48 : Matrix<Complex> o1(i1.shape(), Complex(0.0));
586 48 : Matrix<Complex> o2(i2.shape(), Complex(0.0));
587 42912 : for (Int yy=0; yy < i1.shape()(1); ++yy){
588 : //Int nyy= (Double(yy-cenY)*fratio) + cenY;
589 42864 : Double nyy= (Double((yy-cenY)*convRatio)/fratio) + cenY;
590 42864 : Double cyy=ceil(nyy);
591 42864 : Double fyy= floor(nyy);
592 42864 : Int iy=nyy > fyy+0.5 ? Int(cyy) : Int(fyy);
593 42864 : if(cyy <2*cenY && fyy >=0.0)
594 21076008 : for(Int xx=0; xx < i1.shape()(0); ++ xx){
595 : //Int nxx= Int(Double(xx-cenX)*fratio) + cenX;
596 21056800 : Double nxx= Int(Double((xx-cenX)*convRatio)/fratio) + cenX;
597 21056800 : Double cxx=ceil(nxx);
598 21056800 : Double fxx= floor(nxx);
599 21056800 : Int ix=nxx > fxx+0.5 ? Int(cxx) : Int(fxx) ;
600 21056800 : 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 9925272 : 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 9925272 : o2(xx, yy)=i2(ix, iy);
607 : }
608 : }
609 : }
610 48 : twoDPB.putSlice(o1.reform(elshape), blcin);
611 48 : twoDPB2.putSlice(o2.reform(elshape), blcin);
612 48 : }
613 :
614 20 : }
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 20 : convFunc_p=twoDPB.getSlice(IPosition(4,0,0,0,0), IPosition(4, tempConvSize, tempConvSize, 1, 1), true);
716 20 : weightConvFunc_p.resize();
717 20 : weightConvFunc_p=twoDPB2.getSlice(IPosition(4,0,0,0,nBeamChans-1), IPosition(4, tempConvSize, tempConvSize, 1, 1), true);
718 : //convFunc/=max(abs(convFunc));
719 20 : Float maxAbsConvFunc=max(amplitude(weightConvFunc_p));
720 :
721 20 : Float minAbsConvFunc=min(amplitude(weightConvFunc_p));
722 : //cerr << "min max " << minAbsConvFunc << " " << maxAbsConvFunc << endl;
723 20 : convSupport_p=-1;
724 20 : Bool found=false;
725 : //Bool found2=true;
726 : //Int trial2=0;
727 20 : Int trial=0;
728 7748 : for (trial=tempConvSize/2-2;trial>0;trial--) {
729 : //Searching down a diagonal
730 7748 : if(abs(weightConvFunc_p(tempConvSize/2-trial, tempConvSize/2-trial)) > (1.0e-3*maxAbsConvFunc)) {
731 20 : found=true;
732 20 : trial=Int(sqrt(2.0*Float(trial*trial)));
733 20 : break;
734 : }
735 : }
736 20 : 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 20 : if(trial < 5*convSampling)
744 20 : trial=( tempConvSize > (10*convSampling)) ? 5*convSampling : (tempConvSize/2 - 4*convSampling);
745 :
746 20 : if(found) {
747 20 : 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 20 : Double pbSum=0.0;
762 : //Double pbSum2=0.0;
763 :
764 :
765 :
766 160 : for (Int iy=-convSupport_p;iy<=convSupport_p;iy++) {
767 1600 : for (Int ix=-convSupport_p;ix<=convSupport_p;ix++) {
768 2920 : Complex val=convFunc_p(ix*convSampling+tempConvSize/2,
769 1460 : iy*convSampling+tempConvSize/2);
770 1460 : 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 20 : if(pbSum>0.0) {
780 20 : 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 20 : << LogIO::POST;
791 :
792 20 : convSupportBlock_p.resize(actualConvIndex_p+1);
793 20 : 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 20 : convSupportBlock_p[actualConvIndex_p]=new Vector<Int>(1);
797 20 : convSizes_p[actualConvIndex_p]=new Vector<Int> (1);
798 20 : (*(convSupportBlock_p[actualConvIndex_p]))[0]=convSupport_p;
799 20 : convFunctions_p.resize(actualConvIndex_p+1);
800 20 : convWeights_p.resize(actualConvIndex_p+1);
801 20 : convFunctions_p[actualConvIndex_p]= new Array<Complex>();
802 20 : convWeights_p[actualConvIndex_p]= new Array<Complex>();
803 20 : Int newConvSize=2*(convSupport_p+2)*convSampling;
804 : //NEED to chop this right ...and in the centre
805 20 : if(newConvSize >=tempConvSize)
806 12 : newConvSize=tempConvSize;
807 :
808 20 : IPosition blc(4, (tempConvSize/2)-(newConvSize/2),
809 20 : (tempConvSize/2)-(newConvSize/2), 0, 0);
810 20 : IPosition trc(4, (tempConvSize/2)+(newConvSize/2-1),
811 20 : (tempConvSize/2)+(newConvSize/2-1), 0, nBeamChans-1);
812 20 : 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 20 : convFunctions_p[actualConvIndex_p]->copyMatchingPart(twoDPB.get(false)(blc,trc));//*Complex(1.0/pbSum,0.0));
816 20 : convSize_p=newConvSize;
817 20 : convWeights_p[actualConvIndex_p]->resize(IPosition(5, newConvSize, newConvSize, 1, nBeamChans,1));
818 20 : convWeights_p[actualConvIndex_p]->copyMatchingPart(twoDPB2.get(false)(blc,trc));//*Complex(1.0/pbSum2,0.0));
819 20 : blc.resize(5, false);
820 20 : trc.resize(5,false);
821 20 : blc=IPosition(5, 0, 0, 0, 0,0);
822 20 : trc=IPosition(5, newConvSize-1, newConvSize-1, 0, 0,0);
823 68 : for (Int bc=0; bc< nBeamChans; ++bc){
824 48 : blc[3]=bc;
825 48 : trc[3]=bc;
826 48 : pbSum=real(sum((*convFunctions_p[actualConvIndex_p])(blc,trc)))/Double(convSampling)/Double(convSampling);
827 48 : (*convFunctions_p[actualConvIndex_p])(blc,trc)= (*convFunctions_p[actualConvIndex_p])(blc,trc)/pbSum;
828 48 : (*convWeights_p[actualConvIndex_p])(blc,trc)= (*convWeights_p[actualConvIndex_p])(blc,trc)/pbSum;
829 : }
830 :
831 :
832 :
833 20 : convFunc_p.resize();//break any reference
834 20 : (*convSizes_p[actualConvIndex_p])[0]=convSize_p;
835 20 : doneMainConv_p[actualConvIndex_p]=true;
836 :
837 20 : }
838 : else{
839 13428 : convSize_p=(*convSizes_p[actualConvIndex_p])[0];
840 :
841 : }
842 :
843 : //Apply the shift phase gradient
844 13448 : convFunc.resize();
845 13448 : weightConvFunc.resize();
846 13448 : convFunc.assign(*(convFunctions_p[actualConvIndex_p]));
847 13448 : weightConvFunc.assign(*(convWeights_p[actualConvIndex_p]));
848 : Bool copyconv, copywgt;
849 13448 : Complex *cv=convFunc.getStorage(copyconv);
850 13448 : Complex *wcv=weightConvFunc.getStorage(copywgt);
851 : //cerr << "Field " << vb.fieldId() << " spw " << vb.spectralWindow() << " phase grad: " << pixFieldDir << endl;
852 :
853 45146 : for (Int nc=0; nc < nBeamChans; ++nc){
854 31698 : Int planeoffset=nc*convSize_p*convSize_p;
855 14543850 : for (Int iy=0;iy<convSize_p;iy++) {
856 : Double cy, sy;
857 : Int offset;
858 :
859 14512152 : SINCOS(Double(iy-convSize_p/2)*pixFieldDir(1), sy, cy);
860 14512152 : Complex phy(cy,sy) ;
861 14512152 : offset = iy*convSize_p+planeoffset;
862 6746846664 : for (Int ix=0;ix<convSize_p;ix++) {
863 : Double cx, sx;
864 6732334512 : SINCOS(Double(ix-convSize_p/2)*pixFieldDir(0), sx, cx);
865 6732334512 : Complex phx(cx,sx) ;
866 6732334512 : cv[ix+offset]= cv[ix+offset]*phx*phy;
867 6732334512 : wcv[ix+offset]= wcv[ix+offset]*phx*phy;
868 :
869 : }
870 : }
871 : }
872 13448 : convFunc.putStorage(cv, copyconv);
873 13448 : weightConvFunc.putStorage(wcv, copywgt);
874 13448 : convsize.resize();
875 13448 : convsize=*(convSizes_p[actualConvIndex_p]);
876 13448 : convSupport.resize();
877 13448 : convSupport=(*(convSupportBlock_p[actualConvIndex_p]));
878 :
879 :
880 13448 : }
881 :
882 85077 : void SimplePBConvFunc::setSkyJones(SkyJones* sj){
883 85077 : sj_p=sj;
884 85077 : }
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 85077 : void SimplePBConvFunc::findUsefulChannels(Vector<Int>& chanMap, Vector<Double>& chanFreqs, const vi::VisBuffer2& vb, const Vector<Double>& freq){
911 :
912 :
913 85077 : Int spw=vb.spectralWindows()(0);
914 85077 : bandName_p=vb.subtableColumns().spectralWindow().name()(spw);
915 85077 : Vector<Double> spwfreq=vb.subtableColumns().spectralWindow().chanFreq()(spw);
916 85077 : Double spwfreqwidth=abs(Vector<Double>(vb.subtableColumns().spectralWindow().chanWidth()(spw))(0));
917 85077 : chanMap.resize(freq.nelements());
918 85077 : Vector<Double> localfreq=vb.getFrequencies(0, MFrequency::TOPO);
919 85077 : Double minfreq=min(freq);
920 85077 : Double maxfreq=max(freq);
921 85077 : Double origwidth=freq.nelements()==1 ? 1e12 : (maxfreq-minfreq)/(freq.nelements()-1);
922 : ///Fractional bandwidth which will trigger mutiple PB in one spw
923 :
924 85077 : Double tol=(max(spwfreq))*0.5/100;
925 85077 : if(tol < origwidth/2.0) tol=origwidth/2.0;
926 85077 : Double topFreq=max(spwfreq);
927 333787 : while (topFreq > maxfreq){
928 248710 : topFreq -= tol;
929 : }
930 : // For large tol
931 85077 : if(topFreq < minfreq)
932 47425 : topFreq += tol;
933 : //Int nchan=Int(lround((max(freq)-min(freq))/tol));
934 85077 : Double bottomFreq=topFreq;
935 85077 : Int nchan=0;
936 : //cerr << std::setprecision(12) << "top " << bottomFreq << " minFreq " << minfreq << " maxFreq " << maxfreq << endl;
937 260058 : while(bottomFreq > minfreq){
938 174981 : ++nchan;
939 174981 : bottomFreq -= tol;
940 : }
941 85077 : if(nchan > 1){
942 36212 : nchan-=1;
943 36212 : 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 85077 : if(nchan > vb.nChannels()){
949 0 : nchan=vb.nChannels();
950 0 : tol=spwfreqwidth;
951 0 : bottomFreq=min(localfreq);
952 : }
953 :
954 85077 : chanFreqs.resize();
955 85077 : if(nchan >= (Int)(freq.nelements()-1)) { indgen(chanMap); chanFreqs=freq; return;}
956 280 : if((nchan==0) || (freq.nelements()==1)) { chanFreqs=Vector<Double>(1, bottomFreq);chanMap.set(0); return;}
957 :
958 280 : chanFreqs.resize(nchan);
959 560 : for(Int k=0; k < nchan; ++k){
960 280 : chanFreqs[k]=k*tol+bottomFreq;
961 : }
962 :
963 :
964 280 : Int activechan=0;
965 280 : chanMap.set(0);
966 3376 : for (uInt k=0; k < chanMap.nelements(); ++k){
967 3096 : Double mindiff=DBL_MAX;
968 3096 : Int nearestchan=0;
969 3096 : 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 3096 : if(activechan != nchan){
980 3096 : 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 3096 : activechan=0;
989 : }
990 : //cerr << "chanfreqs " << chanFreqs << endl;
991 : //cerr << "USEFULchan " << chanMap << endl;
992 280 : return;
993 85077 : }
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 13448 : void SimplePBConvFunc::addPBToFlux(const vi::VisBuffer2& vb){
1233 13448 : if(calcFluxScale_p){
1234 13448 : if(fluxScale_p.shape().nelements()==0)
1235 : { //cerr << "nx_p " << nx_p << endl;
1236 0 : calcFluxScale_p=False;
1237 0 : return;
1238 : }
1239 13448 : Vector<Int> pixdepoint(2, -100000);
1240 13448 : convertArray(pixdepoint, thePix_p);
1241 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))){
1242 40 : TempImage<Float> thispb(fluxScale_p.shape(), fluxScale_p.coordinates());
1243 20 : thispb.set(1.0);
1244 20 : sj_p->applySquare(thispb, thispb, vb, 0);
1245 40 : LatticeExpr<Float> le(fluxScale_p+thispb);
1246 20 : fluxScale_p.copyData(le);
1247 20 : 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 20 : }
1270 13448 : }
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 :
|
