Line data Source code
1 : //# HetArrayConvFunc.cc: Implementation for HetArrayConvFunc
2 : //# Copyright (C) 2008-2016
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 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 General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU 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 :
29 : #include <casacore/casa/Arrays/ArrayMath.h>
30 : #include <casacore/casa/Arrays/Array.h>
31 : #include <casacore/casa/Arrays/MaskedArray.h>
32 : #include <casacore/casa/Arrays/Vector.h>
33 : #include <casacore/casa/Arrays/Slice.h>
34 : #include <casacore/casa/Arrays/Matrix.h>
35 : #include <casacore/casa/Arrays/Cube.h>
36 : #include <casacore/scimath/Mathematics/FFTServer.h>
37 : #include <casacore/measures/Measures/MeasTable.h>
38 : #include <casacore/scimath/Mathematics/MathFunc.h>
39 : #include <casacore/scimath/Mathematics/ConvolveGridder.h>
40 : #include <casacore/casa/Utilities/Assert.h>
41 : #include <casacore/casa/Utilities/CompositeNumber.h>
42 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
43 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
44 :
45 : #include <casacore/images/Images/ImageInterface.h>
46 : #include <casacore/images/Images/PagedImage.h>
47 : #include <casacore/images/Images/SubImage.h>
48 : #include <casacore/images/Images/TempImage.h>
49 : #include <imageanalysis/Utilities/SpectralImageUtil.h>
50 : #include <casacore/casa/Logging/LogIO.h>
51 : #include <casacore/casa/Logging/LogSink.h>
52 : #include <casacore/casa/Logging/LogMessage.h>
53 :
54 : #include <casacore/lattices/Lattices/ArrayLattice.h>
55 : #include <casacore/lattices/Lattices/SubLattice.h>
56 : #include <casacore/lattices/LRegions/LCBox.h>
57 : #include <casacore/lattices/Lattices/LatticeConcat.h>
58 : #include <casacore/lattices/LEL/LatticeExpr.h>
59 : #include <casacore/lattices/Lattices/LatticeCache.h>
60 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
61 :
62 :
63 : #include <casacore/ms/MeasurementSets/MSColumns.h>
64 :
65 : #include <msvis/MSVis/VisBuffer2.h>
66 :
67 : #include <synthesis/TransformMachines2/Utils.h>
68 : #include <synthesis/TransformMachines/PBMath1DAiry.h>
69 : #include <synthesis/TransformMachines/PBMath1DNumeric.h>
70 : #include <synthesis/TransformMachines/PBMath2DImage.h>
71 : #include <synthesis/TransformMachines/PBMath.h>
72 : #include <synthesis/TransformMachines2/HetArrayConvFunc.h>
73 : #include <synthesis/MeasurementEquations/VPManager.h>
74 :
75 : #include <casacore/casa/OS/Timer.h>
76 :
77 :
78 :
79 : using namespace casacore;
80 : namespace casa { //# NAMESPACE CASA - BEGIN
81 :
82 : namespace refim {
83 :
84 : using namespace casacore;
85 : using namespace casa;
86 : using namespace casacore;
87 : using namespace casa::refim;
88 :
89 :
90 0 : HetArrayConvFunc::HetArrayConvFunc() : convFunctionMap_p(0), nDefined_p(0),msId_p(-1), actualConvIndex_p(-1), vpTable_p("")
91 : {
92 0 : calcFluxScale_p=true;
93 0 : init(PBMathInterface::AIRY);
94 0 : }
95 :
96 182 : HetArrayConvFunc::HetArrayConvFunc(const PBMathInterface::PBClass typeToUse, const String vpTable):
97 182 : convFunctionMap_p(0), nDefined_p(0), msId_p(-1), actualConvIndex_p(-1), vpTable_p(vpTable)
98 : {
99 182 : calcFluxScale_p=true;
100 182 : init(typeToUse);
101 :
102 182 : }
103 :
104 0 : HetArrayConvFunc::HetArrayConvFunc(const RecordInterface& rec, Bool calcFluxneeded):convFunctionMap_p(0), nDefined_p(0), msId_p(-1), actualConvIndex_p(-1) {
105 0 : String err;
106 0 : fromRecord(err, rec, calcFluxneeded);
107 0 : }
108 :
109 364 : HetArrayConvFunc::~HetArrayConvFunc() {
110 : //
111 364 : }
112 :
113 182 : void HetArrayConvFunc::init(const PBMathInterface::PBClass typeTouse) {
114 182 : doneMainConv_p=false;
115 182 : filledFluxScale_p=false;
116 182 : pbClass_p=typeTouse;
117 182 : }
118 :
119 :
120 :
121 73069 : void HetArrayConvFunc::findAntennaSizes(const vi::VisBuffer2& vb) {
122 :
123 73069 : if(msId_p != vb.msId()) {
124 149 : msId_p=vb.msId();
125 : //MSColumns mscol(vb.ms());
126 149 : const MSAntennaColumns& ac=vb.subtableColumns().antenna();
127 149 : antIndexToDiamIndex_p.resize(ac.nrow());
128 149 : antIndexToDiamIndex_p.set(-1);
129 149 : Int diamIndex=antDiam2IndexMap_p.size( );
130 149 : Vector<Double> dishDiam=ac.dishDiameter().getColumn();
131 149 : Vector<String>dishName=ac.name().getColumn();
132 149 : String telescop=vb.subtableColumns().observation().telescopeName()(0);
133 : PBMath::CommonPB whichPB;
134 149 : if(pbClass_p==PBMathInterface::COMMONPB) {
135 97 : String band;
136 97 : String commonPBName;
137 : // This frequency is ONLY required to determine which PB model to use:
138 : // The VLA, the ATNF, and WSRT have frequency - dependent PB models
139 97 : Quantity freq( vb.subtableColumns().spectralWindow().refFrequency()(0), "Hz");
140 :
141 :
142 97 : PBMath::whichCommonPBtoUse( telescop, freq, band, whichPB, commonPBName );
143 : //Revert to using AIRY for unknown common telescope
144 97 : if(whichPB==PBMath::UNKNOWN)
145 0 : pbClass_p=PBMathInterface::AIRY;
146 :
147 97 : }
148 149 : if(pbClass_p== PBMathInterface::AIRY) {
149 52 : LogIO os;
150 52 : os << LogOrigin("HetArrConvFunc", "findAntennaSizes") << LogIO::NORMAL;
151 : ////////We'll be using dish diameter as key
152 1002 : for (uInt k=0; k < dishDiam.nelements(); ++k) {
153 950 : if( (diamIndex !=0) && antDiam2IndexMap_p.find(String::toString(dishDiam(k))) != antDiam2IndexMap_p.end( ) ) {
154 864 : antIndexToDiamIndex_p(k)=antDiam2IndexMap_p[String::toString(dishDiam(k))];
155 : }
156 : else {
157 86 : if(dishDiam[k] > 0.0) { //there may be stations with no dish on
158 86 : antDiam2IndexMap_p.insert(std::pair<casacore::String, casacore::Int>(String::toString(dishDiam(k)), diamIndex));
159 86 : antIndexToDiamIndex_p(k)=diamIndex;
160 86 : antMath_p.resize(diamIndex+1);
161 86 : if(pbClass_p== PBMathInterface::AIRY) {
162 : //ALMA ratio of blockage to dish
163 86 : Quantity qdiam= Quantity (dishDiam(k),"m");
164 86 : Quantity blockDiam= Quantity(dishDiam(k)/12.0*.75, "m");
165 86 : Quantity support=Quantity(150, "arcsec");
166 : ///For ALMA 12m dish it is effectively 10.7 m according to Todd Hunter
167 : ///@ 2011-12-06
168 86 : if((vb.subtableColumns().observation().telescopeName()(0) =="ALMA") || (vb.subtableColumns().observation().telescopeName()(0) =="ACA")){
169 64 : Quantity fov(max(nx_p*dc_p.increment()(0), ny_p*dc_p.increment()(1)), dc_p.worldAxisUnits()(0));
170 64 : if((abs(dishDiam[k] - 12.0) < 0.5)) {
171 35 : qdiam= Quantity(10.7, "m");
172 35 : blockDiam= Quantity(0.75, "m");
173 35 : support=Quantity(max(150.0, fov.getValue("arcsec")/5.0), "arcsec");
174 :
175 : }
176 : else{
177 : //2017 the ACA dishes are best represented by 6.25m:
178 :
179 29 : qdiam= Quantity(6.25,"m");
180 29 : blockDiam = Quantity(0.75,"m");
181 29 : support=Quantity(max(300.0,fov.getValue("arcsec")/3.0) , "arcsec");
182 : }
183 64 : }
184 86 : os << "Overriding PB with Airy of diam,blockage="<<qdiam<<","<<blockDiam<<" starting with antenna "<<k<<LogIO::POST;
185 :
186 :
187 :
188 :
189 86 : antMath_p[diamIndex]=new PBMath1DAiry(qdiam, blockDiam,
190 : support,
191 172 : Quantity(100.0,"GHz"));
192 86 : }
193 :
194 :
195 :
196 : //////Will no longer support this
197 : /*else if(pbClass_p== PBMathInterface::IMAGE){
198 : //Get the image name by calling code for the antenna name and array name
199 : //For now hard wired to ALMA as this part of the code will not be accessed for non-ALMA
200 : //see Imager::setMosaicFTMachine
201 : // When ready to generalize then code that calls with telescope name, antenna name
202 : //(via vb.msColumns) and/or diameter and frequency via vb.frequency (indexing will need to
203 : //be upgraded to account for frequency too) should be done to return the
204 : //right voltage pattern image.
205 : String vpImageName="";
206 : if (abs(dishDiam[k]-7.0) < 1.0)
207 : Aipsrc::find(vpImageName, "alma.vp.7m", "");
208 : else
209 : Aipsrc::find(vpImageName, "alma.vp.12m", "") ;
210 : //cerr << "first vpImagename " << vpImageName << endl;
211 : if(vpImageName==""){
212 : String beamPath;
213 : if(!MeasTable::AntennaResponsesPath(beamPath, "ALMA")){
214 : throw(AipsError("Alma beam images requested cannot be found "));
215 : }
216 : else{
217 : beamPath=beamPath.before(String("AntennaResponses"));
218 : vpImageName= (abs(dishDiam[k]-7.0) < 1.0) ? beamPath
219 : +String("/ALMA_AIRY_7M.VP") :
220 : beamPath+String("/ALMA_AIRY_12M.VP");
221 : }
222 :
223 :
224 : }
225 : //cerr << "Using the image VPs " << vpImageName << endl;
226 : if(Table::isReadable(vpImageName))
227 : antMath_p[diamIndex]=new PBMath2DImage(PagedImage<Complex>(vpImageName));
228 : else
229 : throw(AipsError(String("Cannot find voltage pattern image ") + vpImageName));
230 : }
231 : else{
232 :
233 : throw(AipsError("Do not deal with non airy dishes or images of VP yet "));
234 : }
235 : */
236 86 : ++diamIndex;
237 :
238 : }
239 :
240 : }
241 : }
242 :
243 52 : }
244 97 : else if(pbClass_p== PBMathInterface::IMAGE) {
245 :
246 0 : VPManager *vpman=VPManager::Instance();
247 0 : if(vpTable_p != String(""))
248 0 : vpman->loadfromtable(vpTable_p);
249 : ///else it is already loaded in the static object
250 0 : Vector<Record> recs;
251 0 : Vector<Vector<String> > antnames;
252 :
253 0 : if(vpman->imagepbinfo(antnames, recs)) {
254 0 : Vector<Bool> dishDefined(dishName.nelements(), false);
255 0 : Int nbeams=antnames.nelements();
256 : ///will be keying on file image file name here
257 0 : for (uInt k=0; k < dishDiam.nelements(); ++k) {
258 0 : String key;
259 0 : Bool beamDone=false;
260 0 : Int recordToUse=0;
261 0 : for (Int j =0; j < nbeams; ++j) {
262 0 : key=recs[j].isDefined("realimage") ? recs[j].asString("realimage") : recs[j].asString("compleximage");
263 0 : if(antnames[j][0]=="*" || anyEQ(dishName[k], antnames[j])) {
264 0 : dishDefined[k]=true;
265 0 : recordToUse=j;
266 :
267 0 : if((diamIndex !=0) && antDiam2IndexMap_p.find(key) != antDiam2IndexMap_p.end( )) {
268 0 : antIndexToDiamIndex_p(k)=antDiam2IndexMap_p[key];
269 0 : beamDone=true;
270 : }
271 : }
272 : }
273 0 : if(!beamDone && dishDefined[k]) {
274 0 : key=recs[recordToUse].isDefined("realimage") ? recs[recordToUse].asString("realimage") : recs[recordToUse].asString("compleximage");
275 0 : antDiam2IndexMap_p.insert(std::pair<casacore::String, casacore::Int>(key, diamIndex));
276 0 : antIndexToDiamIndex_p(k)=diamIndex;
277 0 : antMath_p.resize(diamIndex+1);
278 0 : if(recs[recordToUse].isDefined("realimage") && recs[recordToUse].isDefined("imagimage")) {
279 : //PagedImage<Float> realim(recs[recordToUse].asString("realimage"));
280 : // PagedImage<Float> imagim(recs[recordToUse].asString("imagim"));
281 : // antMath_p[diamIndex]=new PBMath2DImage(realim, imagim);
282 :
283 0 : if(!Table::isReadable(recs[recordToUse].asString("realimage")))
284 0 : throw(AipsError("real part of VP "+recs[recordToUse].asString("realimage")+ " is not readable"));
285 0 : PagedImage<Float> realim(recs[recordToUse].asString("realimage"));
286 0 : CountedPtr<ImageInterface<Float> >imagim;
287 0 : if(recs[recordToUse].asString("imagimage").size()==0){
288 0 : imagim=new TempImage<Float>(realim.shape(), realim.coordinates());
289 0 : imagim->set(0.0);
290 : }
291 : else{
292 0 : if(!Table::isReadable(recs[recordToUse].asString("imagimage")))
293 0 : throw(AipsError("Imaginary part of VP "+recs[recordToUse].asString("imagimage")+ " is not readable"));
294 0 : imagim= new PagedImage<Float> (recs[recordToUse].asString("imagimage"));
295 : }
296 0 : antMath_p[diamIndex]=new PBMath2DImage(realim, *imagim);
297 :
298 0 : }
299 : else {
300 : //antMath_p[diamIndex]=new PBMath2DImage(PagedImage<Complex>(recs[recordToUse].asString("compleximage")));
301 :
302 0 : if(!Table::isReadable(recs[recordToUse].asString("compleximage")))
303 0 : throw(AipsError("complex image of VP "+recs[recordToUse].asString("compleximage")+ " is not readable"));
304 0 : antMath_p[diamIndex]=new PBMath2DImage(PagedImage<Complex>(recs[recordToUse].asString("compleximage")));
305 :
306 : }
307 0 : ++diamIndex;
308 : }
309 0 : }
310 0 : if(!allTrue(dishDefined)) {
311 : //cerr << "dishDefined " << dishDefined << endl;
312 0 : throw(AipsError("Some Antennas in the MS did not have a VP defined"));
313 : }
314 0 : }
315 : else {
316 0 : throw(AipsError("Mosaic does not support non-image voltage patterns yet"));
317 : }
318 :
319 : //Get rid of the static class
320 0 : vpman->reset();
321 0 : }
322 97 : else if(vpTable_p != String("")){
323 : ////When we get vpmanager to give beams on antenna name we
324 : //should change this key to antenna name and loop over all antenna names
325 3 : if((diamIndex !=0) && antDiam2IndexMap_p.find(telescop+String("_")+String::toString(dishDiam(0))) != antDiam2IndexMap_p.end( ) ) {
326 0 : antIndexToDiamIndex_p.set(antDiam2IndexMap_p[telescop+String("_")+String::toString(dishDiam(0))]);
327 : }
328 : else{
329 3 : antDiam2IndexMap_p.insert(std::pair<casacore::String, casacore::Int>(telescop+"_"+String::toString(dishDiam(0)), diamIndex));
330 3 : antIndexToDiamIndex_p.set(diamIndex);
331 3 : VPManager *vpman=VPManager::Instance();
332 3 : vpman->loadfromtable(vpTable_p);
333 3 : Record rec;
334 3 : vpman->getvp(rec, telescop);
335 3 : antMath_p.resize(diamIndex+1);
336 3 : antMath_p[diamIndex]=PBMath::pbMathInterfaceFromRecord(rec);
337 3 : vpman->reset();
338 3 : }
339 :
340 : }
341 94 : else if(pbClass_p==PBMathInterface::COMMONPB) {
342 : //cerr << "Doing the commonPB thing" << endl;
343 : ///Have to use telescop part as string as in multims case different
344 : //telescopes may have same dish size but different commonpb
345 : // VLA and EVLA for e.g.
346 94 : if((diamIndex !=0) && antDiam2IndexMap_p.find(telescop+String("_")+String::toString(dishDiam(0))) != antDiam2IndexMap_p.end( )) {
347 0 : antIndexToDiamIndex_p.set(antDiam2IndexMap_p[telescop+String("_")+String::toString(dishDiam(0))]);
348 : }
349 : else {
350 94 : antDiam2IndexMap_p.insert(std::pair<casacore::String, casacore::Int>(telescop+"_"+String::toString(dishDiam(0)), diamIndex));
351 94 : antIndexToDiamIndex_p.set(diamIndex);
352 94 : antMath_p.resize(diamIndex+1);
353 94 : antMath_p[diamIndex]=PBMath::pbMathInterfaceForCommonPB(whichPB, True);
354 : }
355 : }
356 : else {
357 :
358 0 : throw(AipsError("Mosaic supports image based or Airy voltage patterns or known common pb for now"));
359 :
360 : }
361 :
362 : //cerr << "antIndexTodiamIndex " << antIndexToDiamIndex_p << endl;
363 149 : }
364 :
365 :
366 :
367 :
368 :
369 73069 : }
370 :
371 0 : void HetArrayConvFunc::reset() {
372 0 : doneMainConv_p=false;
373 0 : convFunctions_p.resize(0, true);
374 0 : convWeights_p.resize(0, true);
375 0 : convSizes_p.resize(0, true);
376 0 : convSupportBlock_p.resize(0, true);
377 0 : convFunctionMap_p.resize(0);
378 0 : vbConvIndex_p.clear();
379 0 : ft_p=FFT2D(true);
380 0 : }
381 :
382 :
383 :
384 73069 : void HetArrayConvFunc::findConvFunction(const ImageInterface<Complex>& iimage,
385 : const vi::VisBuffer2& vb,
386 : const Int& convSamp, const Vector<Double>& visFreq,
387 : Array<Complex>& convFunc,
388 : Array<Complex>& weightConvFunc,
389 : Vector<Int>& convsize,
390 : Vector<Int>& convSupport,
391 : Vector<Int>& convFuncPolMap,
392 : Vector<Int>& convFuncChanMap,
393 : Vector<Int>& convFuncRowMap, Bool getConjConvFunc,
394 : const MVDirection& extraShift, const Bool useExtraShift)
395 : {
396 :
397 73069 : storeImageParams(iimage,vb);
398 73069 : convFuncChanMap.resize(vb.nChannels());
399 73069 : Vector<Double> beamFreqs;
400 73069 : findUsefulChannels(convFuncChanMap, beamFreqs, vb, visFreq);
401 73069 : Int nBeamChans=beamFreqs.nelements();
402 : /////For now not doing beam rotation or squints but to be enabled easily
403 73069 : convFuncPolMap.resize(vb.nCorrelations());
404 73069 : Int nBeamPols=1;
405 73069 : convFuncPolMap.set(0);
406 73069 : findAntennaSizes(vb);
407 73069 : uInt ndish=antMath_p.nelements();
408 73069 : if(ndish==0)
409 0 : throw(AipsError("Don't have dishsize"));
410 : Int ndishpair;
411 73069 : if(ndish==1)
412 70201 : ndishpair=1;
413 : else
414 2868 : ndishpair=factorial(ndish)/factorial(ndish-2)/2 + ndish;
415 :
416 73069 : convFunc.resize();
417 73069 : weightConvFunc.resize();
418 73069 : convFuncRowMap.resize();
419 73069 : convsize.resize();
420 73069 : convSupport.resize();
421 :
422 73069 : Int isCached=checkPBOfField(vb, convFuncRowMap, extraShift, useExtraShift);
423 : //cout << "isCached " << isCached << endl;
424 73069 : if(isCached==1 && (convFuncRowMap.shape()[0]==(ssize_t)vb.nRows())) {
425 : /*convFunc.reference(convFunc_p);
426 : weightConvFunc.reference(weightConvFunc_p);
427 : convsize=*convSizes_p[actualConvIndex_p];
428 : convSupport=convSupport_p;
429 : return;
430 : */
431 : }
432 73069 : else if(isCached ==2) {
433 0 : convFunc.resize();
434 0 : weightConvFunc.resize();
435 0 : convsize.resize();
436 0 : convSupport.resize();
437 0 : convFuncRowMap.resize();
438 0 : return;
439 :
440 : }
441 : /////TESTOO elkey
442 146138 : String elkey=String::toString(vb.msId())+String("_")+String::toString(vb.spectralWindows()[0])+String("_")+String::toString(visFreq.nelements());
443 :
444 : /////////////////
445 73069 : actualConvIndex_p=convIndex(vb, visFreq.nelements());
446 : //cerr << "actual conv index " << actualConvIndex_p << " doneMainconv " << doneMainConv_p << endl;
447 73069 : if(doneMainConv_p.shape()[0] < (actualConvIndex_p+1)) {
448 : //cerr << "resizing DONEMAIN " << doneMainConv_p.shape()[0] << endl;
449 269 : doneMainConv_p.resize(actualConvIndex_p+1, true);
450 269 : doneMainConv_p[actualConvIndex_p]=false;
451 269 : convFunctions_p.resize(actualConvIndex_p+1);
452 269 : convFunctions_p[actualConvIndex_p]=nullptr;
453 : }
454 : ///// In multi ms mode ndishpair may change when meeting a new ms
455 : //// redo the calculation then
456 73069 : if(msId_p != vb.msId())//doneMainConv_p[actualConvIndex_p] && ((convSupport_p.nelements() != uInt(ndishpair)) || convFunctions_p[actualConvIndex_p]->shape()[3] != nBeamChans))
457 : {
458 0 : doneMainConv_p[actualConvIndex_p]=false;
459 : //cerr << "invalidating doneMainConv " << convFunctions_p[actualConvIndex_p]->shape()[3] << " =? " << nBeamChans << " convsupp " << convSupport_p.nelements() << endl;
460 : }
461 :
462 : ////Trap for cases when the selection seem to have changed
463 73069 : if(doneMainConv_p[actualConvIndex_p]){
464 72800 : if(nBeamChans > (*convFunctions_p[actualConvIndex_p]).shape()[3])
465 0 : doneMainConv_p[actualConvIndex_p]=False;
466 :
467 : }
468 :
469 :
470 :
471 :
472 : // Get the coordinate system
473 73069 : CoordinateSystem coords(iimage.coordinates());
474 73069 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
475 73069 : AlwaysAssert(directionIndex>=0, AipsError);
476 : // Set up the convolution function.
477 73069 : Int nx=nx_p;
478 73069 : Int ny=ny_p;
479 73069 : Int support=max(nx_p, ny_p)/10;
480 73069 : Int convSampling=1;
481 73069 : if(!doneMainConv_p[actualConvIndex_p]) {
482 578 : for (uInt ii=0; ii < ndish; ++ii) {
483 309 : support=max((antMath_p[ii])->support(coords), support);
484 : }
485 :
486 269 : support=Int(min(Float(support), max(Float(nx_p), Float(ny_p)))*2.0)/2;
487 269 : convSize_p=support*convSampling;
488 : // Make this a nice composite number, to speed up FFTs
489 269 : CompositeNumber cn(uInt(convSize_p*2.0));
490 269 : convSize_p = cn.nextLargerEven(Int(convSize_p));
491 269 : convSize_p=(convSize_p/16)*16; // need it to be divisible by 8 in places
492 :
493 269 : }
494 :
495 :
496 73069 : DirectionCoordinate dc=dc_p;
497 : //where in the image in pixels is this pointing
498 73069 : Vector<Double> pixFieldDir(2);
499 73069 : if(doneMainConv_p.shape()[0] < (actualConvIndex_p+1)) {
500 : //cerr << "resizing DONEMAIN " << doneMainConv_p.shape()[0] << endl;
501 0 : doneMainConv_p.resize(actualConvIndex_p+1, true);
502 0 : doneMainConv_p[actualConvIndex_p]=false;
503 : }
504 : //no need to call toPix here as its been done already above in checkOFPB
505 : //thus the values are still current.
506 73069 : pixFieldDir=thePix_p;
507 : //toPix(pixFieldDir, vb);
508 73069 : MDirection fieldDir=direction1_p;
509 : //shift from center
510 73069 : pixFieldDir(0)=pixFieldDir(0)- Double(nx / 2);
511 73069 : pixFieldDir(1)=pixFieldDir(1)- Double(ny / 2);
512 : //phase gradient per pixel to apply
513 73069 : pixFieldDir(0)=-pixFieldDir(0)*2.0*C::pi/Double(nx)/Double(convSamp);
514 73069 : pixFieldDir(1)=-pixFieldDir(1)*2.0*C::pi/Double(ny)/Double(convSamp);
515 :
516 :
517 73069 : if(!doneMainConv_p[actualConvIndex_p]) {
518 : //cerr << "doneMainConv_p " << actualConvIndex_p << endl;
519 :
520 269 : Vector<Double> sampling;
521 :
522 269 : sampling = dc.increment();
523 269 : sampling*=Double(convSampling);
524 269 : sampling(0)*=Double(nx)/Double(convSize_p);
525 269 : sampling(1)*=Double(ny)/Double(convSize_p);
526 269 : dc.setIncrement(sampling);
527 :
528 269 : Vector<Double> unitVec(2);
529 269 : unitVec=convSize_p/2;
530 269 : dc.setReferencePixel(unitVec);
531 : //make sure we are using the same units
532 269 : fieldDir.set(dc.worldAxisUnits()(0));
533 269 : dc.setReferenceValue(fieldDir.getAngle().getValue());
534 269 : coords.replaceCoordinate(dc, directionIndex);
535 269 : Int spind=coords.findCoordinate(Coordinate::SPECTRAL);
536 269 : SpectralCoordinate spCoord=coords.spectralCoordinate(spind);
537 269 : spCoord.setReferencePixel(Vector<Double>(1,0.0));
538 269 : spCoord.setReferenceValue(Vector<Double>(1, beamFreqs(0)));
539 269 : if(beamFreqs.nelements() >1)
540 68 : spCoord.setIncrement(Vector<Double>(1, beamFreqs(1)-beamFreqs(0)));
541 : //cerr << "spcoord " ;
542 : //spCoord.print(std::cerr);
543 269 : coords.replaceCoordinate(spCoord, spind);
544 269 : CoordinateSystem conjCoord=coords;
545 269 : Double centerFreq=SpectralImageUtil::worldFreq(csys_p, 0.0);
546 269 : SpectralCoordinate conjSpCoord=spCoord;
547 : //cerr << "centreFreq " << centerFreq << " beamFreqs " << beamFreqs(0) << " " << beamFreqs(1) << endl;
548 269 : conjSpCoord.setReferenceValue(Vector<Double>(1,SynthesisUtils::conjFreq(beamFreqs[0], centerFreq)));
549 : ///Increment should go in the reverse direction
550 : ////Do a tabular spectral coordinate for more than 1 channel
551 269 : if(beamFreqs.nelements() >1){
552 68 : Vector<Double> conjFreqs(beamFreqs.nelements());
553 272 : for (uInt kk=0; kk< beamFreqs.nelements(); ++kk){
554 : //conjFreqs[kk]=sqrt(2*centerFreq*centerFreq-beamFreqs(kk)*beamFreqs(kk));
555 204 : conjFreqs[kk]=SynthesisUtils::conjFreq(beamFreqs[kk], centerFreq);
556 : }
557 68 : conjSpCoord=SpectralCoordinate(spCoord.frequencySystem(), conjFreqs, spCoord.restFrequency());
558 : //conjSpCoord.setIncrement(Vector<Double>(1, beamFreqs(0)-beamFreqs(1)));
559 68 : }
560 269 : conjCoord.replaceCoordinate(conjSpCoord, spind);
561 269 : IPosition pbShape(4, convSize_p, convSize_p, 1, nBeamChans);
562 : //TempImage<Complex> twoDPB(pbShape, coords);
563 :
564 :
565 538 : TempLattice<Complex> convFuncTemp(TiledShape(IPosition(5, convSize_p/4, convSize_p/4, nBeamPols, nBeamChans, ndishpair), IPosition(5, convSize_p/4, convSize_p/4, 1, 1, 1)), 0);
566 538 : TempLattice<Complex> weightConvFuncTemp(TiledShape(IPosition(5, convSize_p/4, convSize_p/4, nBeamPols, nBeamChans, ndishpair), IPosition(5, convSize_p/4, convSize_p/4, 1, 1, 1)), 0);
567 : //convFunc_p.resize(IPosition(5, convSize_p, convSize_p, nBeamPols, nBeamChans, ndishpair));
568 :
569 : // convFunc_p=0.0;
570 : //weightConvFunc_p.resize(IPosition(5, convSize_p, convSize_p, nBeamPols, nBeamChans, ndishpair));
571 : //weightConvFunc_p=0.0;
572 269 : IPosition begin(5, 0, 0, 0, 0, 0);
573 538 : IPosition end(5, convFuncTemp.shape()[0]-1, convFuncTemp.shape()[1]-1, nBeamPols-1, nBeamChans-1, 0);
574 : //FFTServer<Float, Complex> fft(IPosition(2, convSize_p, convSize_p));
575 : //TempImage<Complex> pBScreen(TiledShape(pbShape, IPosition(4, convSize_p, convSize_p, 1, 1)), coords, 0);
576 : //TempImage<Complex> pB2Screen(TiledShape(pbShape, IPosition(4, convSize_p, convSize_p, 1, 1)), coords, 0);
577 269 : Long memtot=HostInfo::memoryFree();
578 269 : Double memtobeused= Double(memtot)*1024.0;
579 269 : if(memtot <= 2000000)
580 0 : memtobeused=0.0;
581 269 : TempImage<Complex> pBScreen(TiledShape(pbShape), coords, memtobeused/2.2);
582 :
583 269 : TempImage<Complex> pB2Screen(TiledShape(pbShape), ((nchan_p==1) && getConjConvFunc) ?conjCoord : coords , memtobeused/2.2);
584 269 : IPosition start(4, 0, 0, 0, 0);
585 269 : convSupport_p.resize(ndishpair);
586 : //////////////////
587 : /*Double wtime0=0.0;
588 : Double wtime1=0.0;
589 : Double wtime2=0.0;
590 : wtime0=omp_get_wtime()
591 : */;
592 : //////////////
593 578 : for (uInt k=0; k < ndish; ++k) {
594 :
595 658 : for (uInt j =k ; j < ndish; ++j) {
596 : //Timer tim;
597 : //Matrix<Complex> screen(convSize_p, convSize_p);
598 : //screen=1.0;
599 : //pBScreen.putSlice(screen, start);
600 : //cerr << "k " << k << " shape " << pBScreen.shape() << " direction1 " << direction1_p << " direction2 " << direction2_p << endl;
601 :
602 : //pBScreen.set(Complex(1.0, 0.0));
603 : //one antenna
604 349 : antMath_p[k]->setBandOrFeedName(bandName_p);
605 349 : antMath_p[j]->setBandOrFeedName(bandName_p);
606 349 : IPosition blcin(4, 0, 0, 0, 0);
607 349 : IPosition trcin(4, convSize_p-1, convSize_p-1, 0, 0);
608 986 : for (Int kk=0; kk < nBeamChans; ++kk) {
609 637 : blcin[3]=kk;
610 637 : trcin[3]=kk;
611 : //wtime0=omp_get_wtime();
612 637 : Slicer slin(blcin, trcin, Slicer::endIsLast);
613 637 : SubImage<Complex> subim(pBScreen, slin, true);
614 637 : subim.set(Complex(1.0, 0.0));
615 637 : (antMath_p[k])->applyVP(subim, subim, direction1_p);
616 :
617 : //Then the other
618 637 : (antMath_p[j])->applyVP(subim, subim, direction2_p);
619 : //tim.show("After Apply ");
620 : //tim.mark();
621 : //pB2Screen.set(Complex(1.0,0.0));
622 637 : SubImage<Complex> subim2(pB2Screen, slin, true);
623 637 : subim2.set(Complex(1.0,0.0));
624 :
625 637 : if(nchan_p >1 || !getConjConvFunc){
626 : //one antenna
627 616 : (antMath_p[k])->applyPB(subim2, subim2, direction1_p);
628 : //Then the other
629 616 : (antMath_p[j])->applyPB(subim2, subim2, direction2_p);
630 : }
631 : else{
632 : //direct frequency PB
633 : //cerr << "orig coords " << subim.coordinates().toWorld(IPosition(4,0,0,0,0)) << " conj coords " << subim2.coordinates().toWorld(IPosition(4,0,0,0,0)) << endl;
634 : //cerr << "incr " << subim.coordinates().increment() << " " << subim2.coordinates().increment() << endl;
635 21 : subim2.copyData(subim);
636 : //Now do the conjugate freq multiplication
637 21 : (antMath_p[k])->applyVP(subim2, subim2, direction1_p);
638 :
639 : //Then the other
640 21 : (antMath_p[j])->applyVP(subim2, subim2, direction2_p);
641 :
642 : /*
643 : //one antenna
644 : (antMath_p[k])->applyPB(subim2, subim2, direction1_p);
645 : //Then the other
646 : (antMath_p[j])->applyPB(subim2, subim2, direction2_p);
647 : */
648 : }
649 : //tim.show("After Apply2 ");
650 : //tim.mark();
651 : //wtime1+=omp_get_wtime()-wtime0;
652 : //subim.copyData((LatticeExpr<Complex>) (iif(abs(subim)> 5e-2, subim, 0)));
653 : //subim2.copyData((LatticeExpr<Complex>) (iif(abs(subim2)> 25e-4, subim2, 0)));
654 :
655 : //wtime0=omp_get_wtime();
656 :
657 : //make sure fft2d plan shape is the same or else recalculate it
658 637 : auto [ftx, fty] = ft_p.getShape();
659 637 : if(ftx >0 && fty >0 && (ftx != subim.shape()(0)) && (fty != subim.shape()(1))){
660 2 : ft_p = FFT2D(true);
661 : }
662 :
663 637 : ft_p.c2cFFTInDouble(subim);
664 637 : ft_p.c2cFFTInDouble(subim2);
665 : //ft_p.c2cFFT(subim);
666 : //ft_p.c2cFFT(subim2);
667 : //wtime2+=omp_get_wtime()-wtime0;
668 : // tim.show("after ffts ");
669 :
670 :
671 637 : }
672 : //cerr << "Apply " << wtime1 << " fft " << wtime2 << endl;
673 : /*
674 : if(nBeamChans >1){
675 : Slicer slin(blcin, trcin, Slicer::endIsLast);
676 : SubImage<Complex> origPB(pBScreen, slin, false);
677 : IPosition elshape= origPB.shape();
678 : Matrix<Complex> i1=origPB.get(true);
679 : SubImage<Complex> origPB2(pB2Screen, slin, false);
680 : Matrix<Complex> i2=origPB2.get(true);
681 : Int cenX=i1.shape()(0)/2;
682 : Int cenY=i1.shape()(1)/2;
683 :
684 : for (Int kk=0; kk < (nBeamChans-1); ++kk){
685 : Double fratio=beamFreqs(kk)/beamFreqs(nBeamChans-1);
686 : cerr << "fratio " << fratio << endl;
687 : blcin[3]=kk;
688 : trcin[3]=kk;
689 : //Slicer slout(blcin, trcin, Slicer::endIsLast);
690 : Matrix<Complex> o1(i1.shape(), Complex(0.0));
691 : Matrix<Complex> o2(i2.shape(), Complex(0.0));
692 : for (Int yy=0; yy < i1.shape()(1); ++yy){
693 : //Int nyy= (Double(yy-cenY)*fratio) + cenY;
694 : Double nyy= (Double(yy-cenY)/fratio) + cenY;
695 : Double cyy=ceil(nyy);
696 : Double fyy= floor(nyy);
697 : Int iy=nyy > fyy+0.5 ? cyy : fyy;
698 : if(cyy <2*cenY && fyy >=0.0)
699 : for(Int xx=0; xx < i1.shape()(0); ++ xx){
700 : //Int nxx= Int(Double(xx-cenX)*fratio) + cenX;
701 : Double nxx= Int(Double(xx-cenX)/fratio) + cenX;
702 : Double cxx=ceil(nxx);
703 : Double fxx= floor(nxx);
704 : Int ix=nxx > fxx+0.5 ? cxx : fxx ;
705 : if(cxx < 2*cenX && fxx >=0.0 ){
706 : //Double dist=sqrt((nxx-cxx)*(nxx-cxx)+(nyy-cyy)*(nyy-cyy))/sqrt(2.0);
707 : //o1(xx, yy)=float(1-dist)*i1(fxx, fyy)+ dist*i1(cxx,cyy);
708 : o1(xx, yy)=i1( ix, iy);
709 : //o2(xx, yy)=i2(nxx, nyy);
710 : //o2(xx, yy)=float(1-dist)*i2(fxx, fyy)+ dist*i2(cxx,cyy);
711 : o2(xx, yy)=i2(ix, iy);
712 : }
713 : }
714 : }
715 : pBScreen.putSlice(o1.reform(elshape), blcin);
716 : pB2Screen.putSlice(o2.reform(elshape), blcin);
717 : }
718 :
719 : }
720 : */
721 :
722 : //tim.show("after apply+apply2+masking+fft ");
723 : //tim.mark();
724 : //LatticeFFT::cfft2d(pBScreen);
725 : //LatticeFFT::cfft2d(pB2Screen);
726 :
727 : //Matrix<Complex> lala=pBScreen.get(true);
728 : //fft.fft0(lala, true);
729 : //fft.flip(lala, false, false);
730 : // pBScreen.put(lala.reform(IPosition(4, convSize_p, convSize_p, 1, 1)));
731 : //lala=pB2Screen.get(true);
732 : //fft.fft0(lala, true);
733 : //fft.flip(lala, false, false);
734 : //pB2Screen.put(lala.reform(IPosition(4, convSize_p, convSize_p, 1, 1)));
735 :
736 : //////////*****************
737 : /*if(0){
738 : ostringstream os1;
739 : os1 << "PB_field_" << Int(thePix_p[0]) << "_" << Int(thePix_p[1]) << "_antpair_" << k <<"_"<<j ;
740 : PagedImage<Float> thisScreen(pbShape, coords, String(os1));
741 : LatticeExpr<Float> le(abs(pBScreen));
742 : thisScreen.copyData(le);
743 : }*/
744 : ////////*****************/
745 :
746 : //tim.show("after FFT ");
747 : //tim.mark();
748 349 : Int plane=0;
749 389 : for (uInt jj=0; jj < k; ++jj)
750 40 : plane=plane+ndish-jj-1;
751 349 : plane=plane+j;
752 349 : begin[4]=plane;
753 349 : end[4]=plane;
754 349 : Slicer slplane(begin, end, Slicer::endIsLast);
755 : //cerr << "SHAPES " << convFunc_p(begin, end).shape() << " " << pBScreen.get(false).shape() << " begin and end " << begin << " " << end << endl;
756 : //convFunc_p(begin, end).copyMatchingPart(pBScreen.get(false));
757 : //weightConvFunc_p(begin, end).copyMatchingPart(pB2Screen.get(false));
758 349 : IPosition blcQ(4, pbShape(0)/8*3, pbShape(1)/8*3, 0, 0);
759 349 : IPosition trcQ(4, blcQ[0]+ pbShape(0)/4-1, blcQ[1]+pbShape(1)/4-1 , nBeamPols-1, nBeamChans-1);
760 :
761 : //cerr << "blcQ " << blcQ << " trcQ " << trcQ << " pbShape " << pbShape << endl;
762 349 : Slicer slQ(blcQ, trcQ, Slicer::endIsLast);
763 : {
764 349 : SubImage<Complex> pBSSub(pBScreen, slQ, false);
765 349 : SubLattice<Complex> cFTempSub(convFuncTemp, slplane, true);
766 349 : LatticeConcat<Complex> lc(4);
767 349 : lc.setLattice(pBSSub);
768 : //cerr << "SHAPES " << cFTempSub.shape() << " " << lc.shape() << endl;
769 349 : cFTempSub.copyData(lc);
770 : //cFTempSub.copyData(pBScreen);
771 349 : }
772 : {
773 349 : SubImage<Complex> pB2SSub(pB2Screen, slQ, false);
774 349 : SubLattice<Complex> cFTempSub2(weightConvFuncTemp, slplane, true);
775 349 : LatticeConcat<Complex> lc(4);
776 349 : lc.setLattice(pB2SSub);
777 349 : cFTempSub2.copyData(lc);
778 : // cFTempSub2.copyData(pB2Screen);
779 : //weightConvFuncTemp.putSlice(pB2Screen.get(false), begin);
780 :
781 349 : }
782 : // supportAndNormalize(plane, convSampling);
783 349 : supportAndNormalizeLatt( plane, convSampling, convFuncTemp, weightConvFuncTemp);
784 :
785 :
786 :
787 : // tim.show("After search of support ");
788 349 : }
789 :
790 : }
791 :
792 :
793 269 : doneMainConv_p[actualConvIndex_p]=true;
794 269 : convFunctions_p.resize(actualConvIndex_p+1);
795 269 : convWeights_p.resize(actualConvIndex_p+1);
796 269 : convSupportBlock_p.resize(actualConvIndex_p+1);
797 : //Using conjugate change support to be larger of either
798 269 : if((nchan_p == 1) && getConjConvFunc) {
799 21 : Int conjsupp=conjSupport(beamFreqs) ;
800 21 : if(conjsupp > max(convSupport_p)){
801 6 : convSupport_p=conjsupp;
802 : }
803 :
804 : }
805 269 : convFunctions_p[actualConvIndex_p]= new Array<Complex>();
806 269 : convWeights_p[actualConvIndex_p]= new Array<Complex>();
807 269 : convSupportBlock_p[actualConvIndex_p]=new Vector<Int>();
808 269 : Int newConvSize=2*(max(convSupport_p)+2)*convSampling;
809 269 : Int newRealConvSize=newConvSize* Int(Double(convSamp)/Double(convSampling));
810 269 : Int lattSize=convFuncTemp.shape()(0);
811 269 : (*convSupportBlock_p[actualConvIndex_p])=convSupport_p;
812 538 : LogIO os(LogOrigin("HetArrConvFunc", "findConvFunction", WHERE));
813 269 : os << "convolution function support: " << convSupport_p<< "ELKEY " << elkey << " actualConvInd "<< actualConvIndex_p << " pointer " << this << LogIO::POST;
814 :
815 269 : if(newConvSize < lattSize) {
816 269 : IPosition blc(5, (lattSize/2)-(newConvSize/2),
817 269 : (lattSize/2)-(newConvSize/2),0,0,0);
818 269 : IPosition trc(5, (lattSize/2)+(newConvSize/2-1),
819 269 : (lattSize/2)+(newConvSize/2-1), nBeamPols-1, nBeamChans-1,ndishpair-1);
820 269 : IPosition shp(5, newConvSize, newConvSize, nBeamPols, nBeamChans, ndishpair);
821 :
822 269 : convFunctions_p[actualConvIndex_p]= new Array<Complex>(IPosition(5, newRealConvSize, newRealConvSize, nBeamPols, nBeamChans, ndishpair ));
823 269 : convWeights_p[actualConvIndex_p]= new Array<Complex>(IPosition(5, newRealConvSize, newRealConvSize, nBeamPols, nBeamChans, ndishpair ));
824 269 : (*convFunctions_p[actualConvIndex_p])=resample(convFuncTemp.getSlice(blc,shp),Double(convSamp)/Double(convSampling));
825 269 : convSize_p=newRealConvSize;
826 269 : (*convWeights_p[actualConvIndex_p])=resample(weightConvFuncTemp.getSlice(blc, shp),Double(convSamp)/Double(convSampling));
827 : //cerr << "nchan " << nchan_p << " getconj " << getConjConvFunc << endl;
828 :
829 269 : }
830 : else {
831 0 : newRealConvSize=lattSize* Int(Double(convSamp)/Double(convSampling));
832 0 : convFunctions_p[actualConvIndex_p]= new Array<Complex>(IPosition(5, newRealConvSize, newRealConvSize, nBeamPols, nBeamChans, ndishpair ));
833 0 : convWeights_p[actualConvIndex_p]= new Array<Complex>(IPosition(5, newRealConvSize, newRealConvSize, nBeamPols, nBeamChans, ndishpair ));
834 :
835 0 : (*convFunctions_p[actualConvIndex_p])=resample(convFuncTemp.get(), Double(convSamp)/Double(convSampling));
836 0 : (*convWeights_p[actualConvIndex_p])=resample(weightConvFuncTemp.get(), Double(convSamp)/Double(convSampling));
837 0 : convSize_p=newRealConvSize;
838 : }
839 :
840 :
841 269 : if((nchan_p == 1) && getConjConvFunc) {
842 21 : fillConjConvFunc(beamFreqs);
843 : /////////////////////////TESTOOO
844 : /*PagedImage<Complex> SCREEN2(TiledShape(convFunctions_p[actualConvIndex_p]->shape()), TMP, "CONJU"+String::toString(actualConvIndex_p));
845 : SCREEN2.put(*convFunctionsConjFreq_p[actualConvIndex_p] );
846 : */
847 : ////////////////////////
848 : }
849 :
850 269 : convFunc_p.resize();
851 269 : weightConvFunc_p.resize();
852 :
853 269 : }
854 : else {
855 72800 : convSize_p=max(*convSizes_p[actualConvIndex_p]);
856 72800 : convSupport_p.resize();
857 72800 : convSupport_p=*convSupportBlock_p[actualConvIndex_p];
858 : }
859 : /*
860 : rowMap.resize(vb.nRow());
861 : for (Int k=0; k < vb.nRow(); ++k){
862 : //plane of convfunc that match this pair of antennas
863 : rowMap(k)=antIndexToDiamIndex_p(vb.antenna1()(k))*ndish+
864 : antIndexToDiamIndex_p(vb.antenna2()(k));
865 :
866 : }
867 : */
868 : ////////////////TESTOOO
869 : // CoordinateSystem TMP = coords;
870 : // CoordinateUtil::addLinearAxes(TMP, Vector<String>(1,"gulu"), IPosition(1,nBeamChans));
871 : // PagedImage<Complex> SCREEN(TiledShape(convFunctions_p[actualConvIndex_p]->shape()), TMP, "NONCONJUVI2"+String::toString(actualConvIndex_p));
872 : // SCREEN.put(*convFunctions_p[actualConvIndex_p] );
873 : // PagedImage<Complex> SCREEN3(TiledShape(convWeights_p[actualConvIndex_p]->shape()), TMP, "FTWEIGHTVI2"+String::toString(actualConvIndex_p));
874 : // SCREEN3.put(*convWeights_p[actualConvIndex_p] );
875 :
876 : /////////////////
877 :
878 73069 : makerowmap(vb, convFuncRowMap);
879 : ///need to deal with only the maximum of different baselines available in this
880 : ///vb
881 73069 : ndishpair=max(convFuncRowMap)+1;
882 :
883 : //convSupportBlock_p.resize(actualConvIndex_p+1);
884 73069 : convSizes_p.resize(actualConvIndex_p+1);
885 : //convSupportBlock_p[actualConvIndex_p]=new Vector<Int>(ndishpair);
886 : //(*convSupportBlock_p[actualConvIndex_p])=convSupport_p;
887 73069 : convSizes_p[actualConvIndex_p]=new Vector<Int> (ndishpair);
888 :
889 : /* convFunctions_p[actualConvIndex_p]->resize(convSize_p, convSize_p, ndishpair);
890 : *(convFunctions_p[actualConvIndex_p])=convSave_p;
891 : convWeights_p[actualConvIndex_p]->resize(convSize_p, convSize_p, ndishpair);
892 : *(convWeights_p[actualConvIndex_p])=weightSave_p;
893 : */
894 :
895 73069 : convFunc_p.resize();
896 73069 : if((nchan_p == 1) && getConjConvFunc) {
897 : // cerr << this << " recovering " << actualConvIndex_p << " " <<convFunctionsConjFreq_p.size() << endl;
898 7668 : if(Int(convFunctionsConjFreq_p.size()) <= actualConvIndex_p){
899 0 : fillConjConvFunc(beamFreqs);
900 :
901 : }
902 7668 : convFunc_p=(*convFunctionsConjFreq_p[actualConvIndex_p]);
903 : }
904 : else{
905 :
906 65401 : convFunc_p=(*convFunctions_p[actualConvIndex_p]);
907 : }
908 :
909 73069 : weightConvFunc_p.resize();
910 73069 : weightConvFunc_p=(*convWeights_p[actualConvIndex_p]);
911 :
912 :
913 : // cerr << "convfunc shapes " << convFunc_p.shape() << " " << weightConvFunc_p.shape() << " " << convSize_p << " pol " << nBeamPols << " chan " << nBeamChans << " ndishpair " << ndishpair << endl;
914 : /////Due to a bug in buildCoordSysCore...sometimes an image bigger
915 : ///than the spw selection chosen is made
916 73069 : if(nBeamChans > convFunc_p.shape()[3])
917 0 : nBeamChans = convFunc_p.shape()[3];
918 : //convSupport_p.resize();
919 : //convSupport_p=(*convSupportBlock_p[actualConvIndex_p]);
920 : Bool delc;
921 : Bool delw;
922 73069 : Double dirX=pixFieldDir(0);
923 73069 : Double dirY=pixFieldDir(1);
924 73069 : Complex *convstor=convFunc_p.getStorage(delc);
925 73069 : Complex *weightstor=weightConvFunc_p.getStorage(delw);
926 73069 : Int elconvsize=convSize_p;
927 :
928 73069 : #pragma omp parallel default(none) firstprivate(convstor, weightstor, dirX, dirY, elconvsize, ndishpair, nBeamChans, nBeamPols)
929 : {
930 :
931 : #pragma omp for
932 : for(Int iy=0; iy<elconvsize; ++iy) {
933 : applyGradientToYLine(iy, convstor, weightstor, dirX, dirY, elconvsize, ndishpair, nBeamChans, nBeamPols);
934 :
935 : }
936 : }///End of pragma
937 :
938 73069 : convFunc_p.putStorage(convstor, delc);
939 73069 : weightConvFunc_p.putStorage(weightstor, delw);
940 :
941 :
942 :
943 : //For now all have the same size convsize;
944 73069 : convSizes_p[actualConvIndex_p]->set(convSize_p);
945 :
946 : //We have to get the references right now
947 : // convFunc_p.resize();
948 : //convFunc_p.reference(*convFunctions_p[actualConvIndex_p]);
949 : //weightConvFunc_p.resize();
950 : //weightConvFunc_p.reference(*convWeights_p[actualConvIndex_p]);
951 :
952 73069 : convFunc.reference(convFunc_p);
953 73069 : weightConvFunc.reference(weightConvFunc_p);
954 73069 : convsize=*convSizes_p[actualConvIndex_p];
955 73069 : convSupport=convSupport_p;
956 :
957 :
958 73069 : }
959 0 : void HetArrayConvFunc::rephaseConvFunc(const ImageInterface<Complex>& iimage,
960 : const vi::VisBuffer2& vb,const Int& convSampling,Array<Complex>& convFunc,
961 : Array<Complex>& weightConvFunc, const vector<Int>& polmap, const vector<Int>& chanmap, const vector<Int>& rowmap, const MVDirection& extraShift, const Bool useExtraShift){
962 0 : storeImageParams(iimage,vb);
963 0 : toPix(vb, extraShift, useExtraShift);
964 0 : Vector<Double> pixFieldDir(2);
965 0 : pixFieldDir=thePix_p;
966 0 : pixFieldDir(0)=pixFieldDir(0)- Double(nx_p / 2);
967 0 : pixFieldDir(1)=pixFieldDir(1)- Double(ny_p / 2);
968 0 : pixFieldDir(0)=-pixFieldDir(0)*2.0*C::pi/Double(nx_p)/Double(convSampling);
969 0 : pixFieldDir(1)=-pixFieldDir(1)*2.0*C::pi/Double(ny_p)/Double(convSampling);
970 0 : Int nconvrow=convFunc.shape()(4);
971 0 : Int nconvchan=convFunc.shape()(3);
972 0 : Int nconvpol=convFunc.shape()(2);
973 0 : Int convsize=convFunc.shape()(0);
974 : Bool delc;
975 : Bool delw;
976 0 : Double dirX=pixFieldDir(0);
977 0 : Double dirY=pixFieldDir(1);
978 0 : Complex *convstor=convFunc.getStorage(delc);
979 0 : Complex *weightstor=weightConvFunc.getStorage(delw);
980 : //Vector<Int> pmap(polmap);
981 : //Vector<Int> cmap(chanmap);
982 : //Vector<Int> rmap(rowmap);
983 0 : #pragma omp parallel default(none) firstprivate(convstor, weightstor, dirX, dirY, convsize, nconvrow, nconvchan, nconvpol) shared(polmap, chanmap, rowmap)
984 : {
985 :
986 : #pragma omp for
987 : for(Int iy=0; iy<convsize; ++iy) {
988 : applyGradientToYLine(iy, convstor, weightstor, dirX, dirY, convsize, nconvrow, nconvchan, nconvpol, polmap, chanmap, rowmap);
989 :
990 : }
991 : }///End of pragma
992 0 : convFunc.putStorage(convstor, delc);
993 0 : weightConvFunc.putStorage(weightstor, delw);
994 :
995 0 : }
996 :
997 : typedef unsigned long long ooLong;
998 :
999 17694480 : void HetArrayConvFunc::applyGradientToYLine(const Int iy, Complex*& convFunctions, Complex*& convWeights, const Double pixXdir, const Double pixYdir, Int convSize, const Int ndishpair, const Int nChan, const Int nPol) {
1000 : Double cy, sy;
1001 :
1002 17694480 : SINCOS(Double(iy-convSize/2)*pixYdir, sy, cy);
1003 17694480 : Complex phy(cy,sy) ;
1004 4980275280 : for (Int ix=0; ix<convSize; ix++) {
1005 : Double cx, sx;
1006 4962580800 : SINCOS(Double(ix-convSize/2)*pixXdir, sx, cx);
1007 4962580800 : Complex phx(cx,sx) ;
1008 9925161600 : for (Int ipol=0; ipol< nPol; ++ipol) {
1009 : //Int poloffset=ipol*nChan*ndishpair*convSize*convSize;
1010 16443097600 : for (Int ichan=0; ichan < nChan; ++ichan) {
1011 : //Int chanoffset=ichan*ndishpair*convSize*convSize;
1012 23518074800 : for(Int iz=0; iz <ndishpair; ++iz) {
1013 12037558000 : ooLong index=((ooLong(iz*nChan+ichan)*nPol+ipol)*ooLong(convSize)+ooLong(iy))*ooLong(convSize)+ooLong(ix);
1014 12037558000 : convFunctions[index]= convFunctions[index]*phx*phy;
1015 12037558000 : convWeights[index]= convWeights[index]*phx*phy;
1016 : }
1017 : }
1018 : }
1019 :
1020 : }
1021 17694480 : }
1022 0 : void HetArrayConvFunc::applyGradientToYLine(const Int iy, Complex*& convFunctions, Complex*& convWeights, const Double pixXdir, const Double pixYdir, Int convSize, const Int ndishpair, const Int nChan, const Int nPol, const vector<Int>& polmap, const vector<Int>& chanmap, const vector<Int>& rowmap ) {
1023 : Double cy, sy;
1024 :
1025 0 : SINCOS(Double(iy-convSize/2)*pixYdir, sy, cy);
1026 0 : Complex phy(cy,sy) ;
1027 0 : for (Int ix=0; ix<convSize; ix++) {
1028 : Double cx, sx;
1029 0 : SINCOS(Double(ix-convSize/2)*pixXdir, sx, cx);
1030 0 : Complex phx(cx,sx) ;
1031 0 : for (uint p=0; p< polmap.size(); ++p) {
1032 : //for (uint p=0; p < nPol; ++p) {
1033 0 : Int ipol=polmap[p];
1034 : //Int ipol=p;
1035 0 : for (uint c=0; c < chanmap.size(); ++c) {
1036 : //for (uint c=0; c < nChan; ++c) {
1037 0 : Int ichan=chanmap[c];
1038 : //Int ichan=c;
1039 0 : for (uint z=0; z < rowmap.size(); ++z) {
1040 : //for (uint z=0; z < ndishpair; ++z) {
1041 0 : Int iz=rowmap[z];
1042 : //Int iz=z;
1043 0 : ooLong index=((ooLong(iz*nChan+ichan)*nPol+ipol)*ooLong(convSize)+ooLong(iy))*ooLong(convSize)+ooLong(ix);
1044 0 : convFunctions[index]= convFunctions[index]*phx*phy;
1045 0 : convWeights[index]= convWeights[index]*phx*phy;
1046 : }
1047 : }
1048 : }
1049 :
1050 : }
1051 0 : }
1052 :
1053 21 : Int HetArrayConvFunc::conjSupport(const casacore::Vector<casacore::Double>& freqs){
1054 21 : Double centerFreq=SpectralImageUtil::worldFreq(csys_p, 0.0);
1055 21 : Double maxRatio=-1.0;
1056 42 : for (Int k=0; k < freqs.shape()[0]; ++k) {
1057 : //Double conjFreq=(centerFreq-freqs[k])+centerFreq;
1058 21 : Double conjFreq=SynthesisUtils::conjFreq(freqs[k], centerFreq);
1059 21 : if(maxRatio < conjFreq/freqs[k] )
1060 21 : maxRatio=conjFreq/freqs[k];
1061 : }
1062 21 : return Int(max(convSupport_p)*sqrt(maxRatio)/2.0)*2;
1063 : }
1064 21 : void HetArrayConvFunc::fillConjConvFunc(const Vector<Double>& freqs) {
1065 : //cerr << "Actualconv index " << actualConvIndex_p << endl;
1066 21 : convFunctionsConjFreq_p.resize(actualConvIndex_p+1);
1067 21 : Double centerFreq=SpectralImageUtil::worldFreq(csys_p, 0.0);
1068 21 : IPosition shp=convFunctions_p[actualConvIndex_p]->shape();
1069 21 : convFunctionsConjFreq_p[actualConvIndex_p]=new Array<Complex>(shp, Complex(0.0));
1070 : //cerr << "convsize " << convSize_p << " convsupport " << convSupport_p << endl;
1071 : /*
1072 : Double maxRatio=-1.0;
1073 : for (Int k=0; k < freqs.shape()[0]; ++k) {
1074 : Double conjFreq=(centerFreq-freqs[k])+centerFreq;
1075 : if(maxRatio < conjFreq/freqs[k] )
1076 : maxRatio=conjFreq/freqs[k];
1077 : }
1078 : */
1079 21 : IPosition blc(5,0,0,0,0,0);
1080 21 : IPosition trc=shp-1;
1081 : /*
1082 : IPosition trcOut=trc;
1083 : IPosition trcOut(0)= Int(shp(0)*maxRatio/2.0)*2-1;
1084 : IPosition trcOut(1)= Int(shp(1)*maxRatio/2.0)*2-1;
1085 : */
1086 42 : for (Int k=0; k < freqs.shape()[0]; ++k) {
1087 : //Double conjFreq=(centerFreq-freqs[k])+centerFreq;
1088 21 : Double conjFreq=SynthesisUtils::conjFreq(freqs[k], centerFreq);
1089 21 : blc[3]=k;
1090 21 : trc[3]=k;
1091 : //cerr << "blc " << blc << " trc "<< trc << " ratio " << conjFreq/freqs[k] << endl;
1092 : //Matrix<Complex> convSlice((*convFunctions_p[actualConvIndex_p])(blc, trc).reform(IPosition(2, shp[0], shp[1])));
1093 21 : Array<Complex> convSlice((*convFunctions_p[actualConvIndex_p])(blc, trc));
1094 : //Array<Complex> weightSlice((*convWeights_p[actualConvIndex_p])(blc,trc));
1095 21 : Array<Complex> conjFreqSlice(resample(convSlice, conjFreq/freqs[k]));
1096 21 : Double ratio1= Double(Int(Double(convSlice.shape()(0))*conjFreq/freqs[k]/2.0)*2)/Double(convSlice.shape()(0));
1097 21 : Double ratio2= Double(Int(Double(convSlice.shape()(1))*conjFreq/freqs[k]/2.0)*2)/Double(convSlice.shape()(1));
1098 : //cerr << "resample shape " << conjFreqSlice.shape() << " ratio " << ratio1*ratio2 << " trc " << trc << endl;
1099 21 : Array<Complex> conjSlice=(*convFunctionsConjFreq_p[actualConvIndex_p])(blc, trc);
1100 21 : if(conjFreq > freqs[k]) {
1101 14 : IPosition end=shp-1;
1102 14 : IPosition beg(5,0,0,0,0,0);
1103 14 : beg(0)=(conjFreqSlice.shape()(0)-shp(0))/2;
1104 14 : beg(1)=(conjFreqSlice.shape()(1)-shp(1))/2;
1105 14 : end(0)=beg(0)+shp(0)-1;
1106 14 : end(1)=beg(1)+shp(1)-1;
1107 14 : end[3]=0;
1108 14 : conjSlice=conjFreqSlice(beg, end);
1109 14 : }
1110 : else {
1111 7 : IPosition end=conjFreqSlice.shape()-1;
1112 7 : end[3]=0;
1113 7 : IPosition beg(5,0,0,0,0,0);
1114 7 : beg(0)=(shp(0)-conjFreqSlice.shape()(0))/2;
1115 7 : beg(1)=(shp(1)-conjFreqSlice.shape()(1))/2;
1116 7 : end(0)+=beg(0);
1117 7 : end(1)+=beg(1);
1118 7 : conjSlice(beg, end)=conjFreqSlice;
1119 7 : }
1120 : //cerr << "SUMS " << sum(real(convSlice)) << " new " << sum(real(conjSlice))/ratio1/ratio2 << endl; //" weight " << sum(real(weightSlice))/ratio1/ratio2<< endl;
1121 21 : Complex renorm( 1.0/(ratio1*ratio2),0.0);
1122 21 : conjSlice=conjSlice*renorm;
1123 : //weightSlice=weightSlice*Complex(1.0/(ratio1*ratio2),0.0);
1124 :
1125 21 : }
1126 :
1127 :
1128 21 : }
1129 6 : Bool HetArrayConvFunc::toRecord(RecordInterface& rec) {
1130 :
1131 : try {
1132 6 : rec.define("name", "HetArrayConvFunc");
1133 6 : Int numConv=convFunctions_p.nelements();
1134 6 : rec.define("ndefined", numConv);
1135 : //rec.define("convfunctionmap", convFunctionMap_p);
1136 6 : std::map<String, Int>::iterator it=vbConvIndex_p.begin();
1137 7 : for (Int64 k=0; k < numConv; ++k) {
1138 1 : rec.define("convfunctions"+String::toString(k), *(convFunctions_p[k]));
1139 1 : rec.define("convweights"+String::toString(k), *(convWeights_p[k]));
1140 1 : rec.define("convsizes"+String::toString(k), *(convSizes_p[k]));
1141 1 : rec.define("convsupportblock"+String::toString(k), *(convSupportBlock_p[k]));
1142 1 : rec.define(String("key")+String::toString(k), it->first);
1143 1 : rec.define(String("val")+String::toString(k), it->second);
1144 1 : it++;
1145 : }
1146 6 : rec.define("actualconvindex", actualConvIndex_p);
1147 6 : rec.define("donemainconv", doneMainConv_p);
1148 6 : rec.define("vptable", vpTable_p);
1149 6 : rec.define("pbclass", Int(pbClass_p));
1150 :
1151 : }
1152 0 : catch(AipsError& x) {
1153 0 : return false;
1154 0 : }
1155 6 : return true;
1156 :
1157 : }
1158 :
1159 0 : Bool HetArrayConvFunc::fromRecord(String& err, const RecordInterface& rec, Bool calcfluxscale) {
1160 : //Force pbmath stuff and saved image stuff
1161 0 : nchan_p=0;
1162 0 : msId_p=-1;
1163 : try {
1164 0 : if(!rec.isDefined("name") || rec.asString("name") != "HetArrayConvFunc") {
1165 0 : throw(AipsError("Wrong record to recover HetArray from"));
1166 : }
1167 0 : nDefined_p=rec.asInt("ndefined");
1168 : //rec.get("convfunctionmap", convFunctionMap_p);
1169 0 : convFunctions_p.resize(nDefined_p, true, false);
1170 0 : convSupportBlock_p.resize(nDefined_p, true, false);
1171 0 : convWeights_p.resize(nDefined_p, true, false);
1172 0 : convSizes_p.resize(nDefined_p, true, false);
1173 0 : vbConvIndex_p.erase(vbConvIndex_p.begin(), vbConvIndex_p.end());
1174 0 : for (Int64 k=0; k < nDefined_p; ++k) {
1175 0 : convFunctions_p[k]=new Array<Complex>();
1176 0 : convWeights_p[k]=new Array<Complex>();
1177 0 : convSizes_p[k]=new Vector<Int>();
1178 0 : convSupportBlock_p[k]=new Vector<Int>();
1179 0 : rec.get("convfunctions"+String::toString(k), *(convFunctions_p[k]));
1180 0 : rec.get("convweights"+String::toString(k), *(convWeights_p[k]));
1181 0 : rec.get("convsizes"+String::toString(k), *(convSizes_p[k]));
1182 0 : rec.get("convsupportblock"+String::toString(k), *(convSupportBlock_p[k]));
1183 0 : String key;
1184 : Int val;
1185 0 : rec.get(String("key")+String::toString(k), key);
1186 0 : rec.get(String("val")+String::toString(k), val);
1187 0 : vbConvIndex_p[key]=val;
1188 0 : }
1189 : //Now that we are calculating all phase gradients on the fly ..
1190 : ///we should clean up some and get rid of the cached variables
1191 :
1192 0 : convSize_p= nDefined_p > 0 ? (*(convSizes_p[0]))[0] : 0;
1193 : //convSave_p.resize();
1194 : //rec.get("convsave", convSave_p);
1195 : //weightSave_p.resize();
1196 : //rec.get("weightsave", weightSave_p);
1197 0 : rec.get("vptable", vpTable_p);
1198 0 : rec.get("donemainconv", doneMainConv_p);
1199 : //convSupport_p.resize();
1200 : //rec.get("convsupport", convSupport_p);
1201 0 : pbClass_p=static_cast<PBMathInterface::PBClass>(rec.asInt("pbclass"));
1202 0 : calcFluxScale_p=calcfluxscale;
1203 : }
1204 0 : catch(AipsError& x) {
1205 0 : err=x.getMesg();
1206 0 : return false;
1207 0 : }
1208 :
1209 0 : return true;
1210 : }
1211 :
1212 :
1213 0 : void HetArrayConvFunc::supportAndNormalize(Int plane, Int convSampling) {
1214 :
1215 0 : LogIO os;
1216 0 : os << LogOrigin("HetArrConvFunc", "suppAndNorm") << LogIO::NORMAL;
1217 : // Locate support
1218 0 : Int convSupport=-1;
1219 0 : IPosition begin(5, 0, 0, 0, 0, plane);
1220 0 : IPosition end(5, convFunc_p.shape()[0]-1, convFunc_p.shape()[1]-1, 0, 0, plane);
1221 0 : Matrix<Complex> convPlane(convFunc_p(begin, end).reform(IPosition(2,convFunc_p.shape()[0], convFunc_p.shape()[1]))) ;
1222 0 : Float maxAbsConvFunc=max(amplitude(convPlane));
1223 0 : Float minAbsConvFunc=min(amplitude(convPlane));
1224 0 : Bool found=false;
1225 0 : Int trial=0;
1226 0 : for (trial=convSize_p/2-2; trial>0; trial--) {
1227 : //Searching down a diagonal
1228 0 : if(abs(convPlane(convSize_p/2-trial,convSize_p/2-trial)) > (1.0e-2*maxAbsConvFunc) ) {
1229 0 : found=true;
1230 0 : trial=Int(sqrt(2.0*Float(trial*trial)));
1231 :
1232 0 : break;
1233 : }
1234 : }
1235 0 : if(!found) {
1236 0 : if((maxAbsConvFunc-minAbsConvFunc) > (1.0e-2*maxAbsConvFunc))
1237 0 : found=true;
1238 : // if it drops by more than 2 magnitudes per pixel
1239 0 : trial=( (10*convSampling) < convSize_p) ? 5*convSampling : (convSize_p/2 - 4*convSampling);
1240 : }
1241 :
1242 :
1243 0 : if(found) {
1244 0 : if(trial < 5*convSampling)
1245 0 : trial= ( (10*convSampling) < convSize_p) ? 5*convSampling : (convSize_p/2 - 4*convSampling);
1246 0 : convSupport=Int(0.5+Float(trial)/Float(convSampling))+1;
1247 : //support is really over the edge
1248 0 : if( (convSupport*convSampling) >= convSize_p/2) {
1249 0 : convSupport=convSize_p/2/convSampling-1;
1250 : }
1251 : }
1252 : else {
1253 : /*
1254 : os << "Convolution function is misbehaved - support seems to be zero\n"
1255 : << "Reasons can be: \nThe image definition not covering one or more of the pointings selected \n"
1256 : << "Or no unflagged data in a given pointing"
1257 :
1258 : << LogIO::EXCEPTION;
1259 : */
1260 : //OTF may have flagged stuff ...
1261 0 : convSupport=0;
1262 : }
1263 : //cerr << "trial " << trial << " convSupport " << convSupport << " convSize " << convSize_p << endl;
1264 0 : convSupport_p(plane)=convSupport;
1265 0 : Double pbSum=0.0;
1266 : /*
1267 : Double pbSum1=0.0;
1268 :
1269 : for (Int iy=-convSupport;iy<=convSupport;iy++) {
1270 : for (Int ix=-convSupport;ix<=convSupport;ix++) {
1271 : Complex val=convFunc_p.xyPlane(plane)(ix*convSampling+convSize_p/2,
1272 : iy*convSampling+convSize_p/2);
1273 :
1274 : pbSum1+=sqrt(real(val)*real(val)+ imag(val)*imag(val));
1275 : }
1276 : }
1277 :
1278 : */
1279 0 : if(convSupport >0) {
1280 0 : IPosition blc(2, -convSupport*convSampling+convSize_p/2, -convSupport*convSampling+convSize_p/2);
1281 0 : IPosition trc(2, convSupport*convSampling+convSize_p/2, convSupport*convSampling+convSize_p/2);
1282 0 : for (Int chan=0; chan < convFunc_p.shape()[3]; ++chan) {
1283 0 : begin[3]=chan;
1284 0 : end[3]=chan;
1285 0 : convPlane.resize();
1286 0 : convPlane.reference(convFunc_p(begin, end).reform(IPosition(2,convFunc_p.shape()[0], convFunc_p.shape()[1])));
1287 0 : pbSum=real(sum(convPlane(blc,trc)))/Double(convSampling)/Double(convSampling);
1288 0 : if(pbSum>0.0) {
1289 0 : (convPlane)=convPlane*Complex(1.0/pbSum,0.0);
1290 0 : convPlane.resize();
1291 0 : convPlane.reference(weightConvFunc_p(begin, end).reform(IPosition(2,convFunc_p.shape()[0], convFunc_p.shape()[1])));
1292 :
1293 0 : (convPlane) =(convPlane)*Complex(1.0/pbSum,0.0);
1294 : }
1295 : else {
1296 : os << "Convolution function integral is not positive"
1297 0 : << LogIO::EXCEPTION;
1298 : }
1299 : }
1300 0 : }
1301 : else {
1302 : //no valid convolution for this pointing
1303 0 : for (Int chan=0; chan < convFunc_p.shape()[3]; ++chan) {
1304 0 : begin[3]=chan;
1305 0 : end[3]=chan;
1306 0 : convFunc_p(begin, end).set(Complex(0.0));
1307 0 : weightConvFunc_p(begin, end).set(Complex(0.0));
1308 : //convFunc_p.xyPlane(plane).set(0.0);
1309 : //weightConvFunc_p.xyPlane(plane).set(0.0);
1310 : }
1311 : }
1312 :
1313 0 : }
1314 :
1315 349 : void HetArrayConvFunc::supportAndNormalizeLatt(Int plane, Int convSampling, TempLattice<Complex>& convFuncLat,
1316 : TempLattice<Complex>& weightConvFuncLat) {
1317 :
1318 349 : LogIO os;
1319 349 : os << LogOrigin("HetArrConvFunc", "suppAndNorm") << LogIO::NORMAL;
1320 : // Locate support
1321 349 : Int convSupport=-1;
1322 : ///Use largest channel as highest freq thus largest conv func
1323 349 : IPosition begin(5, 0, 0, 0, convFuncLat.shape()(3)-1, plane);
1324 698 : IPosition shape(5, convFuncLat.shape()[0], convFuncLat.shape()[1], 1, 1, 1);
1325 : //Int convSize=convSize_p;
1326 349 : Int convSize=shape(0);
1327 : ///use FT weightconvlat as it is wider
1328 349 : Matrix<Complex> convPlane=weightConvFuncLat.getSlice(begin, shape, true);
1329 : Float maxAbsConvFunc, minAbsConvFunc;
1330 349 : IPosition minpos, maxpos;
1331 349 : minMax(minAbsConvFunc, maxAbsConvFunc, minpos, maxpos, amplitude(convPlane));
1332 349 : Bool found=false;
1333 349 : Int trial=0;
1334 349 : Float cutlevel=2.5e-2;
1335 : //numeric needs a larger ft
1336 818 : for (uInt k=0; k < antMath_p.nelements() ; ++k){
1337 469 : if((antMath_p[k]->whichPBClass()) == PBMathInterface::NUMERIC)
1338 0 : cutlevel=5e-3;
1339 : }
1340 :
1341 3480 : for (trial=0; trial< (convSize-max(maxpos.asVector())-2); ++trial) {
1342 : ///largest along either axis
1343 : //cerr << "rat1 " << abs(convPlane(maxpos[0]-trial,maxpos[1]))/maxAbsConvFunc << " rat2 " << abs(convPlane(maxpos[0],maxpos[1]-trial))/maxAbsConvFunc << endl;
1344 3480 : if((abs(convPlane(maxpos[0]-trial, maxpos[1])) < (cutlevel*maxAbsConvFunc)) &&(abs(convPlane(maxpos[0],maxpos[1]-trial)) < (cutlevel*maxAbsConvFunc)) )
1345 : {
1346 :
1347 349 : found=true;
1348 : //trial=Int(sqrt(2.0*Float(trial*trial)));
1349 :
1350 349 : break;
1351 : }
1352 : }
1353 349 : if(!found) {
1354 0 : if((maxAbsConvFunc-minAbsConvFunc) > (cutlevel*maxAbsConvFunc))
1355 0 : found=true;
1356 : // if it drops by more than 2 magnitudes per pixel
1357 : //trial=( (10*convSampling) < convSize) ? 5*convSampling : (convSize/2 - 4*convSampling);
1358 0 : trial=convSize/2 - 4*convSampling;
1359 : }
1360 :
1361 349 : if(found) {
1362 349 : if(trial < 5*convSampling)
1363 19 : trial= ( (10*convSampling) < convSize) ? 5*convSampling : (convSize/2 - 4*convSampling);
1364 349 : convSupport=(Int(0.5+Float(trial)/Float(convSampling)))+1 ;
1365 : //cerr << "convsupp " << convSupport << endl;
1366 : //support is really over the edge
1367 349 : if( (convSupport*convSampling) >= convSize/2) {
1368 0 : convSupport=convSize/2/convSampling-1;
1369 : }
1370 : }
1371 : else {
1372 : /*
1373 : os << "Convolution function is misbehaved - support seems to be zero\n"
1374 : << "Reasons can be: \nThe image definition not covering one or more of the pointings selected \n"
1375 : << "Or no unflagged data in a given pointing"
1376 :
1377 : << LogIO::EXCEPTION;
1378 : */
1379 : //OTF may have flagged stuff ...
1380 0 : convSupport=0;
1381 : }
1382 349 : convSupport_p(plane)=convSupport;
1383 349 : Double pbSum=0.0;
1384 : /*
1385 : Double pbSum1=0.0;
1386 :
1387 : for (Int iy=-convSupport;iy<=convSupport;iy++) {
1388 : for (Int ix=-convSupport;ix<=convSupport;ix++) {
1389 : Complex val=convFunc_p.xyPlane(plane)(ix*convSampling+convSize_p/2,
1390 : iy*convSampling+convSize_p/2);
1391 :
1392 : pbSum1+=sqrt(real(val)*real(val)+ imag(val)*imag(val));
1393 : }
1394 : }
1395 :
1396 : */
1397 : //cerr << "convSize_p " << convSize_p << " convSize " << convSize << endl;
1398 349 : if(convSupport >0) {
1399 349 : IPosition blc(2, -convSupport*convSampling+convSize/2, -convSupport*convSampling+convSize/2);
1400 349 : IPosition trc(2, convSupport*convSampling+convSize/2, convSupport*convSampling+convSize/2);
1401 986 : for (Int chan=0; chan < convFuncLat.shape()[3]; ++chan) {
1402 637 : begin[3]=chan;
1403 : //end[3]=chan;
1404 637 : convPlane.resize();
1405 637 : convPlane=convFuncLat.getSlice(begin, shape, true);
1406 637 : pbSum=real(sum(convPlane(blc,trc)))/Double(convSampling)/Double(convSampling);
1407 637 : if(pbSum>0.0) {
1408 637 : (convPlane)=convPlane*Complex(1.0/pbSum,0.0);
1409 637 : convFuncLat.putSlice(convPlane, begin);
1410 637 : convPlane.resize();
1411 637 : convPlane=weightConvFuncLat.getSlice(begin, shape, true);
1412 637 : Double pbSum1=0.0;
1413 637 : pbSum1=real(sum(convPlane(blc,trc)))/Double(convSampling)/Double(convSampling);
1414 637 : (convPlane) =(convPlane)*Complex(1.0/pbSum1,0.0);
1415 637 : weightConvFuncLat.putSlice(convPlane, begin);
1416 : }
1417 : else {
1418 : os << "Convolution function integral is not positive"
1419 0 : << LogIO::EXCEPTION;
1420 : }
1421 : }
1422 349 : }
1423 : else {
1424 : //no valid convolution for this pointing
1425 0 : for (Int chan=0; chan < convFuncLat.shape()[3]; ++chan) {
1426 0 : begin[3]=chan;
1427 : //end[3]=chan;
1428 0 : convPlane.resize(shape[0], shape[1]);
1429 0 : convPlane.set(Complex(0.0));
1430 0 : convFuncLat.putSlice(convPlane, begin);
1431 0 : weightConvFuncLat.putSlice(convPlane, begin);
1432 : //convFunc_p.xyPlane(plane).set(0.0);
1433 : //weightConvFunc_p.xyPlane(plane).set(0.0);
1434 : }
1435 : }
1436 :
1437 349 : }
1438 :
1439 5736 : Int HetArrayConvFunc::factorial(Int n) {
1440 5736 : Int fact=1;
1441 11472 : for (Int k=1; k<=n; ++k)
1442 5736 : fact *=k;
1443 5736 : return fact;
1444 : }
1445 :
1446 :
1447 73069 : Int HetArrayConvFunc::checkPBOfField(const vi::VisBuffer2& vb,
1448 : Vector<Int>& /*rowMap*/, const MVDirection& extraShift, const Bool useExtraShift) {
1449 :
1450 73069 : toPix(vb, extraShift, useExtraShift);
1451 73069 : Vector<Int> pixdepoint(2);
1452 73069 : convertArray(pixdepoint, thePix_p);
1453 146138 : if((pixdepoint(0) < 0) || pixdepoint(0) >= nx_p || pixdepoint(1) < 0 ||
1454 73069 : pixdepoint(1) >=ny_p) {
1455 : //cout << "in pix de point off " << pixdepoint << endl;
1456 0 : return 2;
1457 : }
1458 146138 : String pointingid=String::toString(pixdepoint(0))+"_"+String::toString(pixdepoint(1));
1459 73069 : String msid=vb.msName(true);
1460 :
1461 :
1462 73069 : if(convFunctionMap_p.nelements() == 0) {
1463 143 : convFunctionMap_p.resize(nx_p*ny_p);
1464 143 : convFunctionMap_p.set(-1);
1465 143 : convFunctionMap_p[pixdepoint[1]*nx_p+pixdepoint[0]]=0;
1466 143 : nDefined_p=1;
1467 143 : actualConvIndex_p=0;
1468 143 : return -1;
1469 : }
1470 :
1471 72926 : if(convFunctionMap_p[pixdepoint[1]*nx_p+pixdepoint[0]] <0) {
1472 887 : actualConvIndex_p=nDefined_p;
1473 887 : convFunctionMap_p[pixdepoint[1]*nx_p+pixdepoint[0]]=nDefined_p;
1474 : // ++nDefined_p;
1475 887 : nDefined_p=1;
1476 887 : return -1;
1477 : }
1478 : else {
1479 72039 : actualConvIndex_p=0;
1480 72039 : return -1;
1481 : }
1482 :
1483 : return 1;
1484 :
1485 :
1486 73069 : }
1487 :
1488 73069 : void HetArrayConvFunc::makerowmap(const vi::VisBuffer2& vb,
1489 : Vector<Int>& rowMap) {
1490 :
1491 73069 : uInt ndish=antMath_p.nelements();
1492 73069 : rowMap.resize(vb.nRows());
1493 14034896 : for (rownr_t k=0; k < vb.nRows(); ++k) {
1494 13961827 : Int index1=antIndexToDiamIndex_p(vb.antenna1()(k));
1495 13961827 : Int index2=antIndexToDiamIndex_p(vb.antenna2()(k));
1496 13961827 : if(index2 < index1) {
1497 0 : index1=index2;
1498 0 : index2=antIndexToDiamIndex_p(vb.antenna1()(k));
1499 : }
1500 13961827 : Int plane=0;
1501 14030711 : for (Int jj=0; jj < index1; ++jj)
1502 68884 : plane=plane+ndish-jj-1;
1503 13961827 : plane=plane+index2;
1504 : //plane of convfunc that match this pair of antennas
1505 13961827 : rowMap(k)=plane;
1506 :
1507 : }
1508 :
1509 73069 : }
1510 :
1511 559 : Array<Complex> HetArrayConvFunc::resample(const Array<Complex>& inarray, const Double factor) {
1512 :
1513 559 : Double nx=Double(inarray.shape()(0));
1514 559 : Double ny=Double(inarray.shape()(1));
1515 559 : IPosition shp=inarray.shape();
1516 559 : shp(0)=Int(nx*factor/2.0)*2;
1517 559 : shp(1)=Int(ny*factor/2.0)*2;
1518 559 : Int newNx=shp(0);
1519 559 : Int newNy=shp(1);
1520 :
1521 559 : Array<Complex> out(shp, Complex(0.0));
1522 : // cerr << "SHP " << shp << endl;
1523 :
1524 559 : IPosition incursor=IPosition(inarray.shape().nelements(),1);
1525 559 : incursor[0]=nx;
1526 559 : incursor[1]=ny;
1527 559 : IPosition outcursor=IPosition(inarray.shape().nelements(),1);
1528 559 : outcursor[0]=shp[0];
1529 559 : outcursor[1]=shp[1];
1530 559 : ArrayIterator<Complex> inIt(inarray, IPosition(2,0,1), True);
1531 559 : ArrayIterator<Complex> outIt(out, IPosition(2,0,1),True);
1532 559 : inIt.origin();
1533 559 : outIt.origin();
1534 : //for (zzz=0; zzz< shp.(4); ++zzz){
1535 : // for(yyy=0; yyy< shp.(3); ++yyy){
1536 : // for(xxx=0; xxx< shp.(2); ++xxx){
1537 1854 : while(!inIt.pastEnd()) {
1538 : // cerr << "Iter shape " << inIt.array().shape() << endl;
1539 1295 : Matrix<Complex> inmat;
1540 1295 : inmat=inIt.array();
1541 : //Matrix<Float> leReal=real(Matrix<Complex>(inIt.array()));
1542 : //Matrix<Float> leImag=imag(Matrix<Complex>(inIt.array()));
1543 1295 : Matrix<Float> leReal=real(inmat);
1544 1295 : Matrix<Float> leImag=imag(inmat);
1545 : Bool leRealCopy, leImagCopy;
1546 1295 : Float *realptr=leReal.getStorage(leRealCopy);
1547 1295 : Float *imagptr=leImag.getStorage(leImagCopy);
1548 : Bool isCopy;
1549 1295 : Matrix<Complex> outMat(outIt.array());
1550 1295 : Complex *intPtr=outMat.getStorage(isCopy);
1551 : Float realval, imagval;
1552 : #ifdef _OPENMP
1553 1295 : omp_set_nested(0);
1554 : #endif
1555 1295 : #pragma omp parallel for default(none) private(realval, imagval) firstprivate(intPtr, realptr, imagptr, nx, ny, newNx, newNy) shared(leReal, leImag)
1556 :
1557 : for (Int k =0; k < newNy; ++k) {
1558 : Double y =Double(k)/Double(newNy)*Double(ny);
1559 :
1560 : for (Int j=0; j < newNx; ++j) {
1561 : // Interpolate2D interp(Interpolate2D::LANCZOS);
1562 : Double x=Double(j)/Double(newNx)*Double(nx);
1563 : //interp.interp(realval, where, leReal);
1564 : realval=interpLanczos(x , y, nx, ny,
1565 : realptr);
1566 : imagval=interpLanczos(x , y, nx, ny,
1567 : imagptr);
1568 : //interp.interp(imagval, where, leImag);
1569 : intPtr[k*Int(newNx)+j]=Complex(realval, imagval);
1570 : }
1571 :
1572 : }
1573 1295 : outMat.putStorage(intPtr, isCopy);
1574 1295 : leReal.putStorage(realptr, leRealCopy);
1575 1295 : leImag.putStorage(imagptr, leImagCopy);
1576 1295 : inIt.next();
1577 1295 : outIt.next();
1578 1295 : }
1579 1118 : return out;
1580 559 : }
1581 0 : Matrix <Complex> HetArrayConvFunc::resample2(const Matrix<Complex>& inarray, const Double factor) {
1582 :
1583 0 : Double nx=Double(inarray.shape()(0));
1584 0 : Double ny=Double(inarray.shape()(1));
1585 0 : IPosition shp=inarray.shape();
1586 0 : shp(0)=Int(nx*factor/2.0)*2;
1587 0 : shp(1)=Int(ny*factor/2.0)*2;
1588 :
1589 :
1590 0 : Matrix<Complex> outMat(shp, Complex(0.0));
1591 :
1592 :
1593 :
1594 : {
1595 : //cerr << "Iter shape " << inarray.shape() << endl;
1596 :
1597 0 : Matrix<Float> leReal=real(inarray);
1598 0 : Matrix<Float> leImag=imag(inarray);
1599 : Bool leRealCopy, leImagCopy;
1600 0 : Float *realptr=leReal.getStorage(leRealCopy);
1601 0 : Float *imagptr=leImag.getStorage(leImagCopy);
1602 : Bool isCopy;
1603 0 : Complex *intPtr=outMat.getStorage(isCopy);
1604 : Float realval, imagval;
1605 : #ifdef _OPENMP
1606 : // omp_set_nested(0);
1607 : #endif
1608 : // #pragma omp parallel for default(none) private(realval, imagval) firstprivate(intPtr, realptr, imagptr, nx, ny) shared(leReal, leImag)
1609 :
1610 0 : for (Int k =0; k < shp(1); ++k) {
1611 0 : Double y =Double(k)/Double(shp(1))*Double(ny);
1612 :
1613 0 : for (Int j=0; j < Int(nx*factor); ++j) {
1614 : // Interpolate2D interp(Interpolate2D::LANCZOS);
1615 0 : Double x=Double(j)/Double(factor);
1616 : //interp.interp(realval, where, leReal);
1617 0 : realval=interpLanczos(x , y, nx, ny,
1618 : realptr);
1619 0 : imagval=interpLanczos(x , y, nx, ny,
1620 : imagptr);
1621 : //interp.interp(imagval, where, leImag);
1622 0 : intPtr[k*Int(nx*factor)+j]=Complex(realval, imagval);
1623 : }
1624 :
1625 : }
1626 0 : outMat.putStorage(intPtr, isCopy);
1627 0 : leReal.putStorage(realptr, leRealCopy);
1628 0 : leImag.putStorage(imagptr, leImagCopy);
1629 :
1630 :
1631 0 : }
1632 0 : return outMat;
1633 0 : }
1634 21918286656 : Float HetArrayConvFunc::sinc(const Float x) {
1635 21918286656 : if (x == 0) {
1636 361448448 : return 1;
1637 : }
1638 21556838208 : return sin(C::pi * x) / (C::pi * x);
1639 : }
1640 231377672 : Float HetArrayConvFunc::interpLanczos( const Double& x , const Double& y, const Double& nx, const Double& ny, const Float* data, const Float a) {
1641 231377672 : Double floorx = floor(x);
1642 231377672 : Double floory = floor(y);
1643 231377672 : Float result=0.0;
1644 231377672 : if (floorx < a || floorx >= nx - a || floory < a || floory >= ny - a) {
1645 79167348 : result = 0;
1646 79167348 : return result;
1647 : }
1648 1065472268 : for (Float i = floorx - a + 1; i <= floorx + a; ++i) {
1649 6392833608 : for (Float j = floory - a + 1; j <= floory + a; ++j) {
1650 5479571664 : result += Float(Double(data[Int(j*nx+i)]) * sinc(x - i)*sinc((x-i)/ a) * sinc(y - j)*sinc((y-j)/ a));
1651 : }
1652 : }
1653 152210324 : return result;
1654 : }
1655 :
1656 0 : ImageInterface<Float>& HetArrayConvFunc::getFluxScaleImage() {
1657 0 : if(!calcFluxScale_p)
1658 0 : throw(AipsError("Programmer Error: flux image cannot be retrieved"));
1659 0 : if(!filledFluxScale_p) {
1660 : //The best flux image for a heterogenous array is the weighted coverage
1661 0 : fluxScale_p=TempImage<Float>(IPosition(4, nx_p, ny_p, npol_p, nchan_p), csys_p);
1662 0 : fluxScale_p.copyData(*(convWeightImage_p));
1663 0 : IPosition blc(4,nx_p, ny_p, npol_p, nchan_p);
1664 0 : IPosition trc(4, ny_p, ny_p, npol_p, nchan_p);
1665 0 : blc(0)=0;
1666 0 : blc(1)=0;
1667 0 : trc(0)=nx_p-1;
1668 0 : trc(1)=ny_p-1;
1669 :
1670 0 : for (Int j=0; j < npol_p; ++j) {
1671 0 : for (Int k=0; k < nchan_p ; ++k) {
1672 :
1673 0 : blc(2)=j;
1674 0 : trc(2)=j;
1675 0 : blc(3)=k;
1676 0 : trc(3)=k;
1677 0 : Slicer sl(blc, trc, Slicer::endIsLast);
1678 0 : SubImage<Float> fscalesub(fluxScale_p, sl, true);
1679 : Float planeMax;
1680 0 : LatticeExprNode LEN = max( fscalesub );
1681 0 : planeMax = LEN.getFloat();
1682 0 : if(planeMax !=0) {
1683 0 : fscalesub.copyData( (LatticeExpr<Float>) (fscalesub/planeMax));
1684 :
1685 : }
1686 0 : }
1687 : }
1688 0 : filledFluxScale_p=true;
1689 0 : }
1690 :
1691 :
1692 0 : return fluxScale_p;
1693 :
1694 : }
1695 :
1696 0 : void HetArrayConvFunc::sliceFluxScale(Int npol) {
1697 0 : IPosition fshp=fluxScale_p.shape();
1698 0 : if (fshp(2)>npol) {
1699 0 : npol_p=npol;
1700 : // use first npol planes...
1701 0 : IPosition blc(4,0,0,0,0);
1702 0 : IPosition trc(4,fluxScale_p.shape()(0)-1, fluxScale_p.shape()(1)-1,npol-1,fluxScale_p.shape()(3)-1);
1703 0 : Slicer sl=Slicer(blc, trc, Slicer::endIsLast);
1704 : //writeable if possible
1705 0 : SubImage<Float> fluxScaleSub = SubImage<Float> (fluxScale_p, sl, true);
1706 0 : SubImage<Float> convWeightImageSub = SubImage<Float> (*convWeightImage_p, sl, true);
1707 0 : fluxScale_p = TempImage<Float>(fluxScaleSub.shape(),fluxScaleSub.coordinates());
1708 0 : convWeightImage_p = new TempImage<Float> (convWeightImageSub.shape(),convWeightImageSub.coordinates());
1709 0 : LatticeExpr<Float> le(fluxScaleSub);
1710 0 : fluxScale_p.copyData(le);
1711 0 : LatticeExpr<Float> le2(convWeightImageSub);
1712 0 : convWeightImage_p->copyData(le2);
1713 0 : }
1714 0 : }
1715 : } // namespace refim end
1716 : } //# NAMESPACE CASA - END
1717 :
1718 :
1719 :
1720 :
|
