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 0 : HetArrayConvFunc::HetArrayConvFunc(const PBMathInterface::PBClass typeToUse, const String vpTable):
97 0 : convFunctionMap_p(0), nDefined_p(0), msId_p(-1), actualConvIndex_p(-1), vpTable_p(vpTable)
98 : {
99 0 : calcFluxScale_p=true;
100 0 : init(typeToUse);
101 :
102 0 : }
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 0 : HetArrayConvFunc::~HetArrayConvFunc() {
110 : //
111 0 : }
112 :
113 0 : void HetArrayConvFunc::init(const PBMathInterface::PBClass typeTouse) {
114 0 : doneMainConv_p=false;
115 0 : filledFluxScale_p=false;
116 0 : pbClass_p=typeTouse;
117 0 : }
118 :
119 :
120 :
121 0 : void HetArrayConvFunc::findAntennaSizes(const vi::VisBuffer2& vb) {
122 :
123 0 : if(msId_p != vb.msId()) {
124 0 : msId_p=vb.msId();
125 : //MSColumns mscol(vb.ms());
126 0 : const MSAntennaColumns& ac=vb.subtableColumns().antenna();
127 0 : antIndexToDiamIndex_p.resize(ac.nrow());
128 0 : antIndexToDiamIndex_p.set(-1);
129 0 : Int diamIndex=antDiam2IndexMap_p.size( );
130 0 : Vector<Double> dishDiam=ac.dishDiameter().getColumn();
131 0 : Vector<String>dishName=ac.name().getColumn();
132 0 : String telescop=vb.subtableColumns().observation().telescopeName()(0);
133 : PBMath::CommonPB whichPB;
134 0 : if(pbClass_p==PBMathInterface::COMMONPB) {
135 0 : String band;
136 0 : 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 0 : Quantity freq( vb.subtableColumns().spectralWindow().refFrequency()(0), "Hz");
140 :
141 :
142 0 : PBMath::whichCommonPBtoUse( telescop, freq, band, whichPB, commonPBName );
143 : //Revert to using AIRY for unknown common telescope
144 0 : if(whichPB==PBMath::UNKNOWN)
145 0 : pbClass_p=PBMathInterface::AIRY;
146 :
147 0 : }
148 0 : if(pbClass_p== PBMathInterface::AIRY) {
149 0 : LogIO os;
150 0 : os << LogOrigin("HetArrConvFunc", "findAntennaSizes") << LogIO::NORMAL;
151 : ////////We'll be using dish diameter as key
152 0 : for (uInt k=0; k < dishDiam.nelements(); ++k) {
153 0 : if( (diamIndex !=0) && antDiam2IndexMap_p.find(String::toString(dishDiam(k))) != antDiam2IndexMap_p.end( ) ) {
154 0 : antIndexToDiamIndex_p(k)=antDiam2IndexMap_p[String::toString(dishDiam(k))];
155 : }
156 : else {
157 0 : if(dishDiam[k] > 0.0) { //there may be stations with no dish on
158 0 : antDiam2IndexMap_p.insert(std::pair<casacore::String, casacore::Int>(String::toString(dishDiam(k)), diamIndex));
159 0 : antIndexToDiamIndex_p(k)=diamIndex;
160 0 : antMath_p.resize(diamIndex+1);
161 0 : if(pbClass_p== PBMathInterface::AIRY) {
162 : //ALMA ratio of blockage to dish
163 0 : Quantity qdiam= Quantity (dishDiam(k),"m");
164 0 : Quantity blockDiam= Quantity(dishDiam(k)/12.0*.75, "m");
165 0 : 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 0 : if((vb.subtableColumns().observation().telescopeName()(0) =="ALMA") || (vb.subtableColumns().observation().telescopeName()(0) =="ACA")){
169 0 : Quantity fov(max(nx_p*dc_p.increment()(0), ny_p*dc_p.increment()(1)), dc_p.worldAxisUnits()(0));
170 0 : if((abs(dishDiam[k] - 12.0) < 0.5)) {
171 0 : qdiam= Quantity(10.7, "m");
172 0 : blockDiam= Quantity(0.75, "m");
173 0 : 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 0 : qdiam= Quantity(6.25,"m");
180 0 : blockDiam = Quantity(0.75,"m");
181 0 : support=Quantity(max(300.0,fov.getValue("arcsec")/3.0) , "arcsec");
182 : }
183 0 : }
184 0 : os << "Overriding PB with Airy of diam,blockage="<<qdiam<<","<<blockDiam<<" starting with antenna "<<k<<LogIO::POST;
185 :
186 :
187 :
188 :
189 0 : antMath_p[diamIndex]=new PBMath1DAiry(qdiam, blockDiam,
190 : support,
191 0 : Quantity(100.0,"GHz"));
192 0 : }
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 0 : ++diamIndex;
237 :
238 : }
239 :
240 : }
241 : }
242 :
243 0 : }
244 0 : 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 0 : 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 0 : 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 0 : antDiam2IndexMap_p.insert(std::pair<casacore::String, casacore::Int>(telescop+"_"+String::toString(dishDiam(0)), diamIndex));
330 0 : antIndexToDiamIndex_p.set(diamIndex);
331 0 : VPManager *vpman=VPManager::Instance();
332 0 : vpman->loadfromtable(vpTable_p);
333 0 : Record rec;
334 0 : vpman->getvp(rec, telescop);
335 0 : antMath_p.resize(diamIndex+1);
336 0 : antMath_p[diamIndex]=PBMath::pbMathInterfaceFromRecord(rec);
337 0 : vpman->reset();
338 0 : }
339 :
340 : }
341 0 : 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 0 : 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 0 : antDiam2IndexMap_p.insert(std::pair<casacore::String, casacore::Int>(telescop+"_"+String::toString(dishDiam(0)), diamIndex));
351 0 : antIndexToDiamIndex_p.set(diamIndex);
352 0 : antMath_p.resize(diamIndex+1);
353 0 : 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 0 : }
364 :
365 :
366 :
367 :
368 :
369 0 : }
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 0 : 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 0 : storeImageParams(iimage,vb);
398 0 : convFuncChanMap.resize(vb.nChannels());
399 0 : Vector<Double> beamFreqs;
400 0 : findUsefulChannels(convFuncChanMap, beamFreqs, vb, visFreq);
401 0 : Int nBeamChans=beamFreqs.nelements();
402 : /////For now not doing beam rotation or squints but to be enabled easily
403 0 : convFuncPolMap.resize(vb.nCorrelations());
404 0 : Int nBeamPols=1;
405 0 : convFuncPolMap.set(0);
406 0 : findAntennaSizes(vb);
407 0 : uInt ndish=antMath_p.nelements();
408 0 : if(ndish==0)
409 0 : throw(AipsError("Don't have dishsize"));
410 : Int ndishpair;
411 0 : if(ndish==1)
412 0 : ndishpair=1;
413 : else
414 0 : ndishpair=factorial(ndish)/factorial(ndish-2)/2 + ndish;
415 :
416 0 : convFunc.resize();
417 0 : weightConvFunc.resize();
418 0 : convFuncRowMap.resize();
419 0 : convsize.resize();
420 0 : convSupport.resize();
421 :
422 0 : Int isCached=checkPBOfField(vb, convFuncRowMap, extraShift, useExtraShift);
423 : //cout << "isCached " << isCached << endl;
424 0 : 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 0 : 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 0 : String elkey=String::toString(vb.msId())+String("_")+String::toString(vb.spectralWindows()[0])+String("_")+String::toString(visFreq.nelements());
443 :
444 : /////////////////
445 0 : actualConvIndex_p=convIndex(vb, visFreq.nelements());
446 : //cerr << "actual conv index " << actualConvIndex_p << " doneMainconv " << doneMainConv_p << endl;
447 0 : if(doneMainConv_p.shape()[0] < (actualConvIndex_p+1)) {
448 : //cerr << "resizing DONEMAIN " << doneMainConv_p.shape()[0] << endl;
449 0 : doneMainConv_p.resize(actualConvIndex_p+1, true);
450 0 : doneMainConv_p[actualConvIndex_p]=false;
451 0 : convFunctions_p.resize(actualConvIndex_p+1);
452 0 : 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 0 : 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 0 : if(doneMainConv_p[actualConvIndex_p]){
464 0 : 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 0 : CoordinateSystem coords(iimage.coordinates());
474 0 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
475 0 : AlwaysAssert(directionIndex>=0, AipsError);
476 : // Set up the convolution function.
477 0 : Int nx=nx_p;
478 0 : Int ny=ny_p;
479 0 : Int support=max(nx_p, ny_p)/10;
480 0 : Int convSampling=1;
481 0 : if(!doneMainConv_p[actualConvIndex_p]) {
482 0 : for (uInt ii=0; ii < ndish; ++ii) {
483 0 : support=max((antMath_p[ii])->support(coords), support);
484 : }
485 :
486 0 : support=Int(min(Float(support), max(Float(nx_p), Float(ny_p)))*2.0)/2;
487 0 : convSize_p=support*convSampling;
488 : // Make this a nice composite number, to speed up FFTs
489 0 : CompositeNumber cn(uInt(convSize_p*2.0));
490 0 : convSize_p = cn.nextLargerEven(Int(convSize_p));
491 0 : convSize_p=(convSize_p/16)*16; // need it to be divisible by 8 in places
492 :
493 0 : }
494 :
495 :
496 0 : DirectionCoordinate dc=dc_p;
497 : //where in the image in pixels is this pointing
498 0 : Vector<Double> pixFieldDir(2);
499 0 : 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 0 : pixFieldDir=thePix_p;
507 : //toPix(pixFieldDir, vb);
508 0 : MDirection fieldDir=direction1_p;
509 : //shift from center
510 0 : pixFieldDir(0)=pixFieldDir(0)- Double(nx / 2);
511 0 : pixFieldDir(1)=pixFieldDir(1)- Double(ny / 2);
512 : //phase gradient per pixel to apply
513 0 : pixFieldDir(0)=-pixFieldDir(0)*2.0*C::pi/Double(nx)/Double(convSamp);
514 0 : pixFieldDir(1)=-pixFieldDir(1)*2.0*C::pi/Double(ny)/Double(convSamp);
515 :
516 :
517 0 : if(!doneMainConv_p[actualConvIndex_p]) {
518 : //cerr << "doneMainConv_p " << actualConvIndex_p << endl;
519 :
520 0 : Vector<Double> sampling;
521 :
522 0 : sampling = dc.increment();
523 0 : sampling*=Double(convSampling);
524 0 : sampling(0)*=Double(nx)/Double(convSize_p);
525 0 : sampling(1)*=Double(ny)/Double(convSize_p);
526 0 : dc.setIncrement(sampling);
527 :
528 0 : Vector<Double> unitVec(2);
529 0 : unitVec=convSize_p/2;
530 0 : dc.setReferencePixel(unitVec);
531 : //make sure we are using the same units
532 0 : fieldDir.set(dc.worldAxisUnits()(0));
533 0 : dc.setReferenceValue(fieldDir.getAngle().getValue());
534 0 : coords.replaceCoordinate(dc, directionIndex);
535 0 : Int spind=coords.findCoordinate(Coordinate::SPECTRAL);
536 0 : SpectralCoordinate spCoord=coords.spectralCoordinate(spind);
537 0 : spCoord.setReferencePixel(Vector<Double>(1,0.0));
538 0 : spCoord.setReferenceValue(Vector<Double>(1, beamFreqs(0)));
539 0 : if(beamFreqs.nelements() >1)
540 0 : spCoord.setIncrement(Vector<Double>(1, beamFreqs(1)-beamFreqs(0)));
541 : //cerr << "spcoord " ;
542 : //spCoord.print(std::cerr);
543 0 : coords.replaceCoordinate(spCoord, spind);
544 0 : CoordinateSystem conjCoord=coords;
545 0 : Double centerFreq=SpectralImageUtil::worldFreq(csys_p, 0.0);
546 0 : SpectralCoordinate conjSpCoord=spCoord;
547 : //cerr << "centreFreq " << centerFreq << " beamFreqs " << beamFreqs(0) << " " << beamFreqs(1) << endl;
548 0 : 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 0 : if(beamFreqs.nelements() >1){
552 0 : Vector<Double> conjFreqs(beamFreqs.nelements());
553 0 : for (uInt kk=0; kk< beamFreqs.nelements(); ++kk){
554 : //conjFreqs[kk]=sqrt(2*centerFreq*centerFreq-beamFreqs(kk)*beamFreqs(kk));
555 0 : conjFreqs[kk]=SynthesisUtils::conjFreq(beamFreqs[kk], centerFreq);
556 : }
557 0 : conjSpCoord=SpectralCoordinate(spCoord.frequencySystem(), conjFreqs, spCoord.restFrequency());
558 : //conjSpCoord.setIncrement(Vector<Double>(1, beamFreqs(0)-beamFreqs(1)));
559 0 : }
560 0 : conjCoord.replaceCoordinate(conjSpCoord, spind);
561 0 : IPosition pbShape(4, convSize_p, convSize_p, 1, nBeamChans);
562 : //TempImage<Complex> twoDPB(pbShape, coords);
563 :
564 :
565 0 : 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 0 : 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 0 : IPosition begin(5, 0, 0, 0, 0, 0);
573 0 : 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 0 : Long memtot=HostInfo::memoryFree();
578 0 : Double memtobeused= Double(memtot)*1024.0;
579 0 : if(memtot <= 2000000)
580 0 : memtobeused=0.0;
581 0 : TempImage<Complex> pBScreen(TiledShape(pbShape), coords, memtobeused/2.2);
582 :
583 0 : TempImage<Complex> pB2Screen(TiledShape(pbShape), ((nchan_p==1) && getConjConvFunc) ?conjCoord : coords , memtobeused/2.2);
584 0 : IPosition start(4, 0, 0, 0, 0);
585 0 : 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 0 : for (uInt k=0; k < ndish; ++k) {
594 :
595 0 : 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 0 : antMath_p[k]->setBandOrFeedName(bandName_p);
605 0 : antMath_p[j]->setBandOrFeedName(bandName_p);
606 0 : IPosition blcin(4, 0, 0, 0, 0);
607 0 : IPosition trcin(4, convSize_p-1, convSize_p-1, 0, 0);
608 0 : for (Int kk=0; kk < nBeamChans; ++kk) {
609 0 : blcin[3]=kk;
610 0 : trcin[3]=kk;
611 : //wtime0=omp_get_wtime();
612 0 : Slicer slin(blcin, trcin, Slicer::endIsLast);
613 0 : SubImage<Complex> subim(pBScreen, slin, true);
614 0 : subim.set(Complex(1.0, 0.0));
615 0 : (antMath_p[k])->applyVP(subim, subim, direction1_p);
616 :
617 : //Then the other
618 0 : (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 0 : SubImage<Complex> subim2(pB2Screen, slin, true);
623 0 : subim2.set(Complex(1.0,0.0));
624 :
625 0 : if(nchan_p >1 || !getConjConvFunc){
626 : //one antenna
627 0 : (antMath_p[k])->applyPB(subim2, subim2, direction1_p);
628 : //Then the other
629 0 : (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 0 : subim2.copyData(subim);
636 : //Now do the conjugate freq multiplication
637 0 : (antMath_p[k])->applyVP(subim2, subim2, direction1_p);
638 :
639 : //Then the other
640 0 : (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 0 : ft_p.c2cFFTInDouble(subim);
657 0 : ft_p.c2cFFTInDouble(subim2);
658 : //ft_p.c2cFFT(subim);
659 : //ft_p.c2cFFT(subim2);
660 : //wtime2+=omp_get_wtime()-wtime0;
661 : // tim.show("after ffts ");
662 :
663 :
664 0 : }
665 : //cerr << "Apply " << wtime1 << " fft " << wtime2 << endl;
666 : /*
667 : if(nBeamChans >1){
668 : Slicer slin(blcin, trcin, Slicer::endIsLast);
669 : SubImage<Complex> origPB(pBScreen, slin, false);
670 : IPosition elshape= origPB.shape();
671 : Matrix<Complex> i1=origPB.get(true);
672 : SubImage<Complex> origPB2(pB2Screen, slin, false);
673 : Matrix<Complex> i2=origPB2.get(true);
674 : Int cenX=i1.shape()(0)/2;
675 : Int cenY=i1.shape()(1)/2;
676 :
677 : for (Int kk=0; kk < (nBeamChans-1); ++kk){
678 : Double fratio=beamFreqs(kk)/beamFreqs(nBeamChans-1);
679 : cerr << "fratio " << fratio << endl;
680 : blcin[3]=kk;
681 : trcin[3]=kk;
682 : //Slicer slout(blcin, trcin, Slicer::endIsLast);
683 : Matrix<Complex> o1(i1.shape(), Complex(0.0));
684 : Matrix<Complex> o2(i2.shape(), Complex(0.0));
685 : for (Int yy=0; yy < i1.shape()(1); ++yy){
686 : //Int nyy= (Double(yy-cenY)*fratio) + cenY;
687 : Double nyy= (Double(yy-cenY)/fratio) + cenY;
688 : Double cyy=ceil(nyy);
689 : Double fyy= floor(nyy);
690 : Int iy=nyy > fyy+0.5 ? cyy : fyy;
691 : if(cyy <2*cenY && fyy >=0.0)
692 : for(Int xx=0; xx < i1.shape()(0); ++ xx){
693 : //Int nxx= Int(Double(xx-cenX)*fratio) + cenX;
694 : Double nxx= Int(Double(xx-cenX)/fratio) + cenX;
695 : Double cxx=ceil(nxx);
696 : Double fxx= floor(nxx);
697 : Int ix=nxx > fxx+0.5 ? cxx : fxx ;
698 : if(cxx < 2*cenX && fxx >=0.0 ){
699 : //Double dist=sqrt((nxx-cxx)*(nxx-cxx)+(nyy-cyy)*(nyy-cyy))/sqrt(2.0);
700 : //o1(xx, yy)=float(1-dist)*i1(fxx, fyy)+ dist*i1(cxx,cyy);
701 : o1(xx, yy)=i1( ix, iy);
702 : //o2(xx, yy)=i2(nxx, nyy);
703 : //o2(xx, yy)=float(1-dist)*i2(fxx, fyy)+ dist*i2(cxx,cyy);
704 : o2(xx, yy)=i2(ix, iy);
705 : }
706 : }
707 : }
708 : pBScreen.putSlice(o1.reform(elshape), blcin);
709 : pB2Screen.putSlice(o2.reform(elshape), blcin);
710 : }
711 :
712 : }
713 : */
714 :
715 : //tim.show("after apply+apply2+masking+fft ");
716 : //tim.mark();
717 : //LatticeFFT::cfft2d(pBScreen);
718 : //LatticeFFT::cfft2d(pB2Screen);
719 :
720 : //Matrix<Complex> lala=pBScreen.get(true);
721 : //fft.fft0(lala, true);
722 : //fft.flip(lala, false, false);
723 : // pBScreen.put(lala.reform(IPosition(4, convSize_p, convSize_p, 1, 1)));
724 : //lala=pB2Screen.get(true);
725 : //fft.fft0(lala, true);
726 : //fft.flip(lala, false, false);
727 : //pB2Screen.put(lala.reform(IPosition(4, convSize_p, convSize_p, 1, 1)));
728 :
729 : //////////*****************
730 : /*if(0){
731 : ostringstream os1;
732 : os1 << "PB_field_" << Int(thePix_p[0]) << "_" << Int(thePix_p[1]) << "_antpair_" << k <<"_"<<j ;
733 : PagedImage<Float> thisScreen(pbShape, coords, String(os1));
734 : LatticeExpr<Float> le(abs(pBScreen));
735 : thisScreen.copyData(le);
736 : }*/
737 : ////////*****************/
738 :
739 : //tim.show("after FFT ");
740 : //tim.mark();
741 0 : Int plane=0;
742 0 : for (uInt jj=0; jj < k; ++jj)
743 0 : plane=plane+ndish-jj-1;
744 0 : plane=plane+j;
745 0 : begin[4]=plane;
746 0 : end[4]=plane;
747 0 : Slicer slplane(begin, end, Slicer::endIsLast);
748 : //cerr << "SHAPES " << convFunc_p(begin, end).shape() << " " << pBScreen.get(false).shape() << " begin and end " << begin << " " << end << endl;
749 : //convFunc_p(begin, end).copyMatchingPart(pBScreen.get(false));
750 : //weightConvFunc_p(begin, end).copyMatchingPart(pB2Screen.get(false));
751 0 : IPosition blcQ(4, pbShape(0)/8*3, pbShape(1)/8*3, 0, 0);
752 0 : IPosition trcQ(4, blcQ[0]+ pbShape(0)/4-1, blcQ[1]+pbShape(1)/4-1 , nBeamPols-1, nBeamChans-1);
753 :
754 : //cerr << "blcQ " << blcQ << " trcQ " << trcQ << " pbShape " << pbShape << endl;
755 0 : Slicer slQ(blcQ, trcQ, Slicer::endIsLast);
756 : {
757 0 : SubImage<Complex> pBSSub(pBScreen, slQ, false);
758 0 : SubLattice<Complex> cFTempSub(convFuncTemp, slplane, true);
759 0 : LatticeConcat<Complex> lc(4);
760 0 : lc.setLattice(pBSSub);
761 : //cerr << "SHAPES " << cFTempSub.shape() << " " << lc.shape() << endl;
762 0 : cFTempSub.copyData(lc);
763 : //cFTempSub.copyData(pBScreen);
764 0 : }
765 : {
766 0 : SubImage<Complex> pB2SSub(pB2Screen, slQ, false);
767 0 : SubLattice<Complex> cFTempSub2(weightConvFuncTemp, slplane, true);
768 0 : LatticeConcat<Complex> lc(4);
769 0 : lc.setLattice(pB2SSub);
770 0 : cFTempSub2.copyData(lc);
771 : // cFTempSub2.copyData(pB2Screen);
772 : //weightConvFuncTemp.putSlice(pB2Screen.get(false), begin);
773 :
774 0 : }
775 : // supportAndNormalize(plane, convSampling);
776 0 : supportAndNormalizeLatt( plane, convSampling, convFuncTemp, weightConvFuncTemp);
777 :
778 :
779 :
780 : // tim.show("After search of support ");
781 0 : }
782 :
783 : }
784 :
785 :
786 0 : doneMainConv_p[actualConvIndex_p]=true;
787 0 : convFunctions_p.resize(actualConvIndex_p+1);
788 0 : convWeights_p.resize(actualConvIndex_p+1);
789 0 : convSupportBlock_p.resize(actualConvIndex_p+1);
790 : //Using conjugate change support to be larger of either
791 0 : if((nchan_p == 1) && getConjConvFunc) {
792 0 : Int conjsupp=conjSupport(beamFreqs) ;
793 0 : if(conjsupp > max(convSupport_p)){
794 0 : convSupport_p=conjsupp;
795 : }
796 :
797 : }
798 0 : convFunctions_p[actualConvIndex_p]= new Array<Complex>();
799 0 : convWeights_p[actualConvIndex_p]= new Array<Complex>();
800 0 : convSupportBlock_p[actualConvIndex_p]=new Vector<Int>();
801 0 : Int newConvSize=2*(max(convSupport_p)+2)*convSampling;
802 0 : Int newRealConvSize=newConvSize* Int(Double(convSamp)/Double(convSampling));
803 0 : Int lattSize=convFuncTemp.shape()(0);
804 0 : (*convSupportBlock_p[actualConvIndex_p])=convSupport_p;
805 0 : LogIO os(LogOrigin("HetArrConvFunc", "findConvFunction", WHERE));
806 0 : os << "convolution function support: " << convSupport_p<< "ELKEY " << elkey << " actualConvInd "<< actualConvIndex_p << " pointer " << this << LogIO::POST;
807 :
808 0 : if(newConvSize < lattSize) {
809 0 : IPosition blc(5, (lattSize/2)-(newConvSize/2),
810 0 : (lattSize/2)-(newConvSize/2),0,0,0);
811 0 : IPosition trc(5, (lattSize/2)+(newConvSize/2-1),
812 0 : (lattSize/2)+(newConvSize/2-1), nBeamPols-1, nBeamChans-1,ndishpair-1);
813 0 : IPosition shp(5, newConvSize, newConvSize, nBeamPols, nBeamChans, ndishpair);
814 :
815 0 : convFunctions_p[actualConvIndex_p]= new Array<Complex>(IPosition(5, newRealConvSize, newRealConvSize, nBeamPols, nBeamChans, ndishpair ));
816 0 : convWeights_p[actualConvIndex_p]= new Array<Complex>(IPosition(5, newRealConvSize, newRealConvSize, nBeamPols, nBeamChans, ndishpair ));
817 0 : (*convFunctions_p[actualConvIndex_p])=resample(convFuncTemp.getSlice(blc,shp),Double(convSamp)/Double(convSampling));
818 0 : convSize_p=newRealConvSize;
819 0 : (*convWeights_p[actualConvIndex_p])=resample(weightConvFuncTemp.getSlice(blc, shp),Double(convSamp)/Double(convSampling));
820 : //cerr << "nchan " << nchan_p << " getconj " << getConjConvFunc << endl;
821 :
822 0 : }
823 : else {
824 0 : newRealConvSize=lattSize* Int(Double(convSamp)/Double(convSampling));
825 0 : convFunctions_p[actualConvIndex_p]= new Array<Complex>(IPosition(5, newRealConvSize, newRealConvSize, nBeamPols, nBeamChans, ndishpair ));
826 0 : convWeights_p[actualConvIndex_p]= new Array<Complex>(IPosition(5, newRealConvSize, newRealConvSize, nBeamPols, nBeamChans, ndishpair ));
827 :
828 0 : (*convFunctions_p[actualConvIndex_p])=resample(convFuncTemp.get(), Double(convSamp)/Double(convSampling));
829 0 : (*convWeights_p[actualConvIndex_p])=resample(weightConvFuncTemp.get(), Double(convSamp)/Double(convSampling));
830 0 : convSize_p=newRealConvSize;
831 : }
832 :
833 :
834 0 : if((nchan_p == 1) && getConjConvFunc) {
835 0 : fillConjConvFunc(beamFreqs);
836 : /////////////////////////TESTOOO
837 : /*PagedImage<Complex> SCREEN2(TiledShape(convFunctions_p[actualConvIndex_p]->shape()), TMP, "CONJU"+String::toString(actualConvIndex_p));
838 : SCREEN2.put(*convFunctionsConjFreq_p[actualConvIndex_p] );
839 : */
840 : ////////////////////////
841 : }
842 :
843 0 : convFunc_p.resize();
844 0 : weightConvFunc_p.resize();
845 :
846 0 : }
847 : else {
848 0 : convSize_p=max(*convSizes_p[actualConvIndex_p]);
849 0 : convSupport_p.resize();
850 0 : convSupport_p=*convSupportBlock_p[actualConvIndex_p];
851 : }
852 : /*
853 : rowMap.resize(vb.nRow());
854 : for (Int k=0; k < vb.nRow(); ++k){
855 : //plane of convfunc that match this pair of antennas
856 : rowMap(k)=antIndexToDiamIndex_p(vb.antenna1()(k))*ndish+
857 : antIndexToDiamIndex_p(vb.antenna2()(k));
858 :
859 : }
860 : */
861 : ////////////////TESTOOO
862 : // CoordinateSystem TMP = coords;
863 : // CoordinateUtil::addLinearAxes(TMP, Vector<String>(1,"gulu"), IPosition(1,nBeamChans));
864 : // PagedImage<Complex> SCREEN(TiledShape(convFunctions_p[actualConvIndex_p]->shape()), TMP, "NONCONJUVI2"+String::toString(actualConvIndex_p));
865 : // SCREEN.put(*convFunctions_p[actualConvIndex_p] );
866 : // PagedImage<Complex> SCREEN3(TiledShape(convWeights_p[actualConvIndex_p]->shape()), TMP, "FTWEIGHTVI2"+String::toString(actualConvIndex_p));
867 : // SCREEN3.put(*convWeights_p[actualConvIndex_p] );
868 :
869 : /////////////////
870 :
871 0 : makerowmap(vb, convFuncRowMap);
872 : ///need to deal with only the maximum of different baselines available in this
873 : ///vb
874 0 : ndishpair=max(convFuncRowMap)+1;
875 :
876 : //convSupportBlock_p.resize(actualConvIndex_p+1);
877 0 : convSizes_p.resize(actualConvIndex_p+1);
878 : //convSupportBlock_p[actualConvIndex_p]=new Vector<Int>(ndishpair);
879 : //(*convSupportBlock_p[actualConvIndex_p])=convSupport_p;
880 0 : convSizes_p[actualConvIndex_p]=new Vector<Int> (ndishpair);
881 :
882 : /* convFunctions_p[actualConvIndex_p]->resize(convSize_p, convSize_p, ndishpair);
883 : *(convFunctions_p[actualConvIndex_p])=convSave_p;
884 : convWeights_p[actualConvIndex_p]->resize(convSize_p, convSize_p, ndishpair);
885 : *(convWeights_p[actualConvIndex_p])=weightSave_p;
886 : */
887 :
888 0 : convFunc_p.resize();
889 0 : if((nchan_p == 1) && getConjConvFunc) {
890 : // cerr << this << " recovering " << actualConvIndex_p << " " <<convFunctionsConjFreq_p.size() << endl;
891 0 : if(Int(convFunctionsConjFreq_p.size()) <= actualConvIndex_p){
892 0 : fillConjConvFunc(beamFreqs);
893 :
894 : }
895 0 : convFunc_p=(*convFunctionsConjFreq_p[actualConvIndex_p]);
896 : }
897 : else{
898 :
899 0 : convFunc_p=(*convFunctions_p[actualConvIndex_p]);
900 : }
901 :
902 0 : weightConvFunc_p.resize();
903 0 : weightConvFunc_p=(*convWeights_p[actualConvIndex_p]);
904 :
905 :
906 : // cerr << "convfunc shapes " << convFunc_p.shape() << " " << weightConvFunc_p.shape() << " " << convSize_p << " pol " << nBeamPols << " chan " << nBeamChans << " ndishpair " << ndishpair << endl;
907 : /////Due to a bug in buildCoordSysCore...sometimes an image bigger
908 : ///than the spw selection chosen is made
909 0 : if(nBeamChans > convFunc_p.shape()[3])
910 0 : nBeamChans = convFunc_p.shape()[3];
911 : //convSupport_p.resize();
912 : //convSupport_p=(*convSupportBlock_p[actualConvIndex_p]);
913 : Bool delc;
914 : Bool delw;
915 0 : Double dirX=pixFieldDir(0);
916 0 : Double dirY=pixFieldDir(1);
917 0 : Complex *convstor=convFunc_p.getStorage(delc);
918 0 : Complex *weightstor=weightConvFunc_p.getStorage(delw);
919 0 : Int elconvsize=convSize_p;
920 :
921 0 : #pragma omp parallel default(none) firstprivate(convstor, weightstor, dirX, dirY, elconvsize, ndishpair, nBeamChans, nBeamPols)
922 : {
923 :
924 : #pragma omp for
925 : for(Int iy=0; iy<elconvsize; ++iy) {
926 : applyGradientToYLine(iy, convstor, weightstor, dirX, dirY, elconvsize, ndishpair, nBeamChans, nBeamPols);
927 :
928 : }
929 : }///End of pragma
930 :
931 0 : convFunc_p.putStorage(convstor, delc);
932 0 : weightConvFunc_p.putStorage(weightstor, delw);
933 :
934 :
935 :
936 : //For now all have the same size convsize;
937 0 : convSizes_p[actualConvIndex_p]->set(convSize_p);
938 :
939 : //We have to get the references right now
940 : // convFunc_p.resize();
941 : //convFunc_p.reference(*convFunctions_p[actualConvIndex_p]);
942 : //weightConvFunc_p.resize();
943 : //weightConvFunc_p.reference(*convWeights_p[actualConvIndex_p]);
944 :
945 0 : convFunc.reference(convFunc_p);
946 0 : weightConvFunc.reference(weightConvFunc_p);
947 0 : convsize=*convSizes_p[actualConvIndex_p];
948 0 : convSupport=convSupport_p;
949 :
950 :
951 0 : }
952 0 : void HetArrayConvFunc::rephaseConvFunc(const ImageInterface<Complex>& iimage,
953 : const vi::VisBuffer2& vb,const Int& convSampling,Array<Complex>& convFunc,
954 : Array<Complex>& weightConvFunc, const vector<Int>& polmap, const vector<Int>& chanmap, const vector<Int>& rowmap, const MVDirection& extraShift, const Bool useExtraShift){
955 0 : storeImageParams(iimage,vb);
956 0 : toPix(vb, extraShift, useExtraShift);
957 0 : Vector<Double> pixFieldDir(2);
958 0 : pixFieldDir=thePix_p;
959 0 : pixFieldDir(0)=pixFieldDir(0)- Double(nx_p / 2);
960 0 : pixFieldDir(1)=pixFieldDir(1)- Double(ny_p / 2);
961 0 : pixFieldDir(0)=-pixFieldDir(0)*2.0*C::pi/Double(nx_p)/Double(convSampling);
962 0 : pixFieldDir(1)=-pixFieldDir(1)*2.0*C::pi/Double(ny_p)/Double(convSampling);
963 0 : Int nconvrow=convFunc.shape()(4);
964 0 : Int nconvchan=convFunc.shape()(3);
965 0 : Int nconvpol=convFunc.shape()(2);
966 0 : Int convsize=convFunc.shape()(0);
967 : Bool delc;
968 : Bool delw;
969 0 : Double dirX=pixFieldDir(0);
970 0 : Double dirY=pixFieldDir(1);
971 0 : Complex *convstor=convFunc.getStorage(delc);
972 0 : Complex *weightstor=weightConvFunc.getStorage(delw);
973 : //Vector<Int> pmap(polmap);
974 : //Vector<Int> cmap(chanmap);
975 : //Vector<Int> rmap(rowmap);
976 0 : #pragma omp parallel default(none) firstprivate(convstor, weightstor, dirX, dirY, convsize, nconvrow, nconvchan, nconvpol) shared(polmap, chanmap, rowmap)
977 : {
978 :
979 : #pragma omp for
980 : for(Int iy=0; iy<convsize; ++iy) {
981 : applyGradientToYLine(iy, convstor, weightstor, dirX, dirY, convsize, nconvrow, nconvchan, nconvpol, polmap, chanmap, rowmap);
982 :
983 : }
984 : }///End of pragma
985 0 : convFunc.putStorage(convstor, delc);
986 0 : weightConvFunc.putStorage(weightstor, delw);
987 :
988 0 : }
989 :
990 : typedef unsigned long long ooLong;
991 :
992 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) {
993 : Double cy, sy;
994 :
995 0 : SINCOS(Double(iy-convSize/2)*pixYdir, sy, cy);
996 0 : Complex phy(cy,sy) ;
997 0 : for (Int ix=0; ix<convSize; ix++) {
998 : Double cx, sx;
999 0 : SINCOS(Double(ix-convSize/2)*pixXdir, sx, cx);
1000 0 : Complex phx(cx,sx) ;
1001 0 : for (Int ipol=0; ipol< nPol; ++ipol) {
1002 : //Int poloffset=ipol*nChan*ndishpair*convSize*convSize;
1003 0 : for (Int ichan=0; ichan < nChan; ++ichan) {
1004 : //Int chanoffset=ichan*ndishpair*convSize*convSize;
1005 0 : for(Int iz=0; iz <ndishpair; ++iz) {
1006 0 : ooLong index=((ooLong(iz*nChan+ichan)*nPol+ipol)*ooLong(convSize)+ooLong(iy))*ooLong(convSize)+ooLong(ix);
1007 0 : convFunctions[index]= convFunctions[index]*phx*phy;
1008 0 : convWeights[index]= convWeights[index]*phx*phy;
1009 : }
1010 : }
1011 : }
1012 :
1013 : }
1014 0 : }
1015 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 ) {
1016 : Double cy, sy;
1017 :
1018 0 : SINCOS(Double(iy-convSize/2)*pixYdir, sy, cy);
1019 0 : Complex phy(cy,sy) ;
1020 0 : for (Int ix=0; ix<convSize; ix++) {
1021 : Double cx, sx;
1022 0 : SINCOS(Double(ix-convSize/2)*pixXdir, sx, cx);
1023 0 : Complex phx(cx,sx) ;
1024 0 : for (uint p=0; p< polmap.size(); ++p) {
1025 : //for (uint p=0; p < nPol; ++p) {
1026 0 : Int ipol=polmap[p];
1027 : //Int ipol=p;
1028 0 : for (uint c=0; c < chanmap.size(); ++c) {
1029 : //for (uint c=0; c < nChan; ++c) {
1030 0 : Int ichan=chanmap[c];
1031 : //Int ichan=c;
1032 0 : for (uint z=0; z < rowmap.size(); ++z) {
1033 : //for (uint z=0; z < ndishpair; ++z) {
1034 0 : Int iz=rowmap[z];
1035 : //Int iz=z;
1036 0 : ooLong index=((ooLong(iz*nChan+ichan)*nPol+ipol)*ooLong(convSize)+ooLong(iy))*ooLong(convSize)+ooLong(ix);
1037 0 : convFunctions[index]= convFunctions[index]*phx*phy;
1038 0 : convWeights[index]= convWeights[index]*phx*phy;
1039 : }
1040 : }
1041 : }
1042 :
1043 : }
1044 0 : }
1045 :
1046 0 : Int HetArrayConvFunc::conjSupport(const casacore::Vector<casacore::Double>& freqs){
1047 0 : Double centerFreq=SpectralImageUtil::worldFreq(csys_p, 0.0);
1048 0 : Double maxRatio=-1.0;
1049 0 : for (Int k=0; k < freqs.shape()[0]; ++k) {
1050 : //Double conjFreq=(centerFreq-freqs[k])+centerFreq;
1051 0 : Double conjFreq=SynthesisUtils::conjFreq(freqs[k], centerFreq);
1052 0 : if(maxRatio < conjFreq/freqs[k] )
1053 0 : maxRatio=conjFreq/freqs[k];
1054 : }
1055 0 : return Int(max(convSupport_p)*sqrt(maxRatio)/2.0)*2;
1056 : }
1057 0 : void HetArrayConvFunc::fillConjConvFunc(const Vector<Double>& freqs) {
1058 : //cerr << "Actualconv index " << actualConvIndex_p << endl;
1059 0 : convFunctionsConjFreq_p.resize(actualConvIndex_p+1);
1060 0 : Double centerFreq=SpectralImageUtil::worldFreq(csys_p, 0.0);
1061 0 : IPosition shp=convFunctions_p[actualConvIndex_p]->shape();
1062 0 : convFunctionsConjFreq_p[actualConvIndex_p]=new Array<Complex>(shp, Complex(0.0));
1063 : //cerr << "convsize " << convSize_p << " convsupport " << convSupport_p << endl;
1064 : /*
1065 : Double maxRatio=-1.0;
1066 : for (Int k=0; k < freqs.shape()[0]; ++k) {
1067 : Double conjFreq=(centerFreq-freqs[k])+centerFreq;
1068 : if(maxRatio < conjFreq/freqs[k] )
1069 : maxRatio=conjFreq/freqs[k];
1070 : }
1071 : */
1072 0 : IPosition blc(5,0,0,0,0,0);
1073 0 : IPosition trc=shp-1;
1074 : /*
1075 : IPosition trcOut=trc;
1076 : IPosition trcOut(0)= Int(shp(0)*maxRatio/2.0)*2-1;
1077 : IPosition trcOut(1)= Int(shp(1)*maxRatio/2.0)*2-1;
1078 : */
1079 0 : for (Int k=0; k < freqs.shape()[0]; ++k) {
1080 : //Double conjFreq=(centerFreq-freqs[k])+centerFreq;
1081 0 : Double conjFreq=SynthesisUtils::conjFreq(freqs[k], centerFreq);
1082 0 : blc[3]=k;
1083 0 : trc[3]=k;
1084 : //cerr << "blc " << blc << " trc "<< trc << " ratio " << conjFreq/freqs[k] << endl;
1085 : //Matrix<Complex> convSlice((*convFunctions_p[actualConvIndex_p])(blc, trc).reform(IPosition(2, shp[0], shp[1])));
1086 0 : Array<Complex> convSlice((*convFunctions_p[actualConvIndex_p])(blc, trc));
1087 : //Array<Complex> weightSlice((*convWeights_p[actualConvIndex_p])(blc,trc));
1088 0 : Array<Complex> conjFreqSlice(resample(convSlice, conjFreq/freqs[k]));
1089 0 : Double ratio1= Double(Int(Double(convSlice.shape()(0))*conjFreq/freqs[k]/2.0)*2)/Double(convSlice.shape()(0));
1090 0 : Double ratio2= Double(Int(Double(convSlice.shape()(1))*conjFreq/freqs[k]/2.0)*2)/Double(convSlice.shape()(1));
1091 : //cerr << "resample shape " << conjFreqSlice.shape() << " ratio " << ratio1*ratio2 << " trc " << trc << endl;
1092 0 : Array<Complex> conjSlice=(*convFunctionsConjFreq_p[actualConvIndex_p])(blc, trc);
1093 0 : if(conjFreq > freqs[k]) {
1094 0 : IPosition end=shp-1;
1095 0 : IPosition beg(5,0,0,0,0,0);
1096 0 : beg(0)=(conjFreqSlice.shape()(0)-shp(0))/2;
1097 0 : beg(1)=(conjFreqSlice.shape()(1)-shp(1))/2;
1098 0 : end(0)=beg(0)+shp(0)-1;
1099 0 : end(1)=beg(1)+shp(1)-1;
1100 0 : end[3]=0;
1101 0 : conjSlice=conjFreqSlice(beg, end);
1102 0 : }
1103 : else {
1104 0 : IPosition end=conjFreqSlice.shape()-1;
1105 0 : end[3]=0;
1106 0 : IPosition beg(5,0,0,0,0,0);
1107 0 : beg(0)=(shp(0)-conjFreqSlice.shape()(0))/2;
1108 0 : beg(1)=(shp(1)-conjFreqSlice.shape()(1))/2;
1109 0 : end(0)+=beg(0);
1110 0 : end(1)+=beg(1);
1111 0 : conjSlice(beg, end)=conjFreqSlice;
1112 0 : }
1113 : //cerr << "SUMS " << sum(real(convSlice)) << " new " << sum(real(conjSlice))/ratio1/ratio2 << endl; //" weight " << sum(real(weightSlice))/ratio1/ratio2<< endl;
1114 0 : Complex renorm( 1.0/(ratio1*ratio2),0.0);
1115 0 : conjSlice=conjSlice*renorm;
1116 : //weightSlice=weightSlice*Complex(1.0/(ratio1*ratio2),0.0);
1117 :
1118 0 : }
1119 :
1120 :
1121 0 : }
1122 0 : Bool HetArrayConvFunc::toRecord(RecordInterface& rec) {
1123 :
1124 : try {
1125 0 : rec.define("name", "HetArrayConvFunc");
1126 0 : Int numConv=convFunctions_p.nelements();
1127 0 : rec.define("ndefined", numConv);
1128 : //rec.define("convfunctionmap", convFunctionMap_p);
1129 0 : std::map<String, Int>::iterator it=vbConvIndex_p.begin();
1130 0 : for (Int64 k=0; k < numConv; ++k) {
1131 0 : rec.define("convfunctions"+String::toString(k), *(convFunctions_p[k]));
1132 0 : rec.define("convweights"+String::toString(k), *(convWeights_p[k]));
1133 0 : rec.define("convsizes"+String::toString(k), *(convSizes_p[k]));
1134 0 : rec.define("convsupportblock"+String::toString(k), *(convSupportBlock_p[k]));
1135 0 : rec.define(String("key")+String::toString(k), it->first);
1136 0 : rec.define(String("val")+String::toString(k), it->second);
1137 0 : it++;
1138 : }
1139 0 : rec.define("actualconvindex", actualConvIndex_p);
1140 0 : rec.define("donemainconv", doneMainConv_p);
1141 0 : rec.define("vptable", vpTable_p);
1142 0 : rec.define("pbclass", Int(pbClass_p));
1143 :
1144 : }
1145 0 : catch(AipsError& x) {
1146 0 : return false;
1147 0 : }
1148 0 : return true;
1149 :
1150 : }
1151 :
1152 0 : Bool HetArrayConvFunc::fromRecord(String& err, const RecordInterface& rec, Bool calcfluxscale) {
1153 : //Force pbmath stuff and saved image stuff
1154 0 : nchan_p=0;
1155 0 : msId_p=-1;
1156 : try {
1157 0 : if(!rec.isDefined("name") || rec.asString("name") != "HetArrayConvFunc") {
1158 0 : throw(AipsError("Wrong record to recover HetArray from"));
1159 : }
1160 0 : nDefined_p=rec.asInt("ndefined");
1161 : //rec.get("convfunctionmap", convFunctionMap_p);
1162 0 : convFunctions_p.resize(nDefined_p, true, false);
1163 0 : convSupportBlock_p.resize(nDefined_p, true, false);
1164 0 : convWeights_p.resize(nDefined_p, true, false);
1165 0 : convSizes_p.resize(nDefined_p, true, false);
1166 0 : vbConvIndex_p.erase(vbConvIndex_p.begin(), vbConvIndex_p.end());
1167 0 : for (Int64 k=0; k < nDefined_p; ++k) {
1168 0 : convFunctions_p[k]=new Array<Complex>();
1169 0 : convWeights_p[k]=new Array<Complex>();
1170 0 : convSizes_p[k]=new Vector<Int>();
1171 0 : convSupportBlock_p[k]=new Vector<Int>();
1172 0 : rec.get("convfunctions"+String::toString(k), *(convFunctions_p[k]));
1173 0 : rec.get("convweights"+String::toString(k), *(convWeights_p[k]));
1174 0 : rec.get("convsizes"+String::toString(k), *(convSizes_p[k]));
1175 0 : rec.get("convsupportblock"+String::toString(k), *(convSupportBlock_p[k]));
1176 0 : String key;
1177 : Int val;
1178 0 : rec.get(String("key")+String::toString(k), key);
1179 0 : rec.get(String("val")+String::toString(k), val);
1180 0 : vbConvIndex_p[key]=val;
1181 0 : }
1182 : //Now that we are calculating all phase gradients on the fly ..
1183 : ///we should clean up some and get rid of the cached variables
1184 :
1185 0 : convSize_p= nDefined_p > 0 ? (*(convSizes_p[0]))[0] : 0;
1186 : //convSave_p.resize();
1187 : //rec.get("convsave", convSave_p);
1188 : //weightSave_p.resize();
1189 : //rec.get("weightsave", weightSave_p);
1190 0 : rec.get("vptable", vpTable_p);
1191 0 : rec.get("donemainconv", doneMainConv_p);
1192 : //convSupport_p.resize();
1193 : //rec.get("convsupport", convSupport_p);
1194 0 : pbClass_p=static_cast<PBMathInterface::PBClass>(rec.asInt("pbclass"));
1195 0 : calcFluxScale_p=calcfluxscale;
1196 : }
1197 0 : catch(AipsError& x) {
1198 0 : err=x.getMesg();
1199 0 : return false;
1200 0 : }
1201 :
1202 0 : return true;
1203 : }
1204 :
1205 :
1206 0 : void HetArrayConvFunc::supportAndNormalize(Int plane, Int convSampling) {
1207 :
1208 0 : LogIO os;
1209 0 : os << LogOrigin("HetArrConvFunc", "suppAndNorm") << LogIO::NORMAL;
1210 : // Locate support
1211 0 : Int convSupport=-1;
1212 0 : IPosition begin(5, 0, 0, 0, 0, plane);
1213 0 : IPosition end(5, convFunc_p.shape()[0]-1, convFunc_p.shape()[1]-1, 0, 0, plane);
1214 0 : Matrix<Complex> convPlane(convFunc_p(begin, end).reform(IPosition(2,convFunc_p.shape()[0], convFunc_p.shape()[1]))) ;
1215 0 : Float maxAbsConvFunc=max(amplitude(convPlane));
1216 0 : Float minAbsConvFunc=min(amplitude(convPlane));
1217 0 : Bool found=false;
1218 0 : Int trial=0;
1219 0 : for (trial=convSize_p/2-2; trial>0; trial--) {
1220 : //Searching down a diagonal
1221 0 : if(abs(convPlane(convSize_p/2-trial,convSize_p/2-trial)) > (1.0e-2*maxAbsConvFunc) ) {
1222 0 : found=true;
1223 0 : trial=Int(sqrt(2.0*Float(trial*trial)));
1224 :
1225 0 : break;
1226 : }
1227 : }
1228 0 : if(!found) {
1229 0 : if((maxAbsConvFunc-minAbsConvFunc) > (1.0e-2*maxAbsConvFunc))
1230 0 : found=true;
1231 : // if it drops by more than 2 magnitudes per pixel
1232 0 : trial=( (10*convSampling) < convSize_p) ? 5*convSampling : (convSize_p/2 - 4*convSampling);
1233 : }
1234 :
1235 :
1236 0 : if(found) {
1237 0 : if(trial < 5*convSampling)
1238 0 : trial= ( (10*convSampling) < convSize_p) ? 5*convSampling : (convSize_p/2 - 4*convSampling);
1239 0 : convSupport=Int(0.5+Float(trial)/Float(convSampling))+1;
1240 : //support is really over the edge
1241 0 : if( (convSupport*convSampling) >= convSize_p/2) {
1242 0 : convSupport=convSize_p/2/convSampling-1;
1243 : }
1244 : }
1245 : else {
1246 : /*
1247 : os << "Convolution function is misbehaved - support seems to be zero\n"
1248 : << "Reasons can be: \nThe image definition not covering one or more of the pointings selected \n"
1249 : << "Or no unflagged data in a given pointing"
1250 :
1251 : << LogIO::EXCEPTION;
1252 : */
1253 : //OTF may have flagged stuff ...
1254 0 : convSupport=0;
1255 : }
1256 : //cerr << "trial " << trial << " convSupport " << convSupport << " convSize " << convSize_p << endl;
1257 0 : convSupport_p(plane)=convSupport;
1258 0 : Double pbSum=0.0;
1259 : /*
1260 : Double pbSum1=0.0;
1261 :
1262 : for (Int iy=-convSupport;iy<=convSupport;iy++) {
1263 : for (Int ix=-convSupport;ix<=convSupport;ix++) {
1264 : Complex val=convFunc_p.xyPlane(plane)(ix*convSampling+convSize_p/2,
1265 : iy*convSampling+convSize_p/2);
1266 :
1267 : pbSum1+=sqrt(real(val)*real(val)+ imag(val)*imag(val));
1268 : }
1269 : }
1270 :
1271 : */
1272 0 : if(convSupport >0) {
1273 0 : IPosition blc(2, -convSupport*convSampling+convSize_p/2, -convSupport*convSampling+convSize_p/2);
1274 0 : IPosition trc(2, convSupport*convSampling+convSize_p/2, convSupport*convSampling+convSize_p/2);
1275 0 : for (Int chan=0; chan < convFunc_p.shape()[3]; ++chan) {
1276 0 : begin[3]=chan;
1277 0 : end[3]=chan;
1278 0 : convPlane.resize();
1279 0 : convPlane.reference(convFunc_p(begin, end).reform(IPosition(2,convFunc_p.shape()[0], convFunc_p.shape()[1])));
1280 0 : pbSum=real(sum(convPlane(blc,trc)))/Double(convSampling)/Double(convSampling);
1281 0 : if(pbSum>0.0) {
1282 0 : (convPlane)=convPlane*Complex(1.0/pbSum,0.0);
1283 0 : convPlane.resize();
1284 0 : convPlane.reference(weightConvFunc_p(begin, end).reform(IPosition(2,convFunc_p.shape()[0], convFunc_p.shape()[1])));
1285 :
1286 0 : (convPlane) =(convPlane)*Complex(1.0/pbSum,0.0);
1287 : }
1288 : else {
1289 : os << "Convolution function integral is not positive"
1290 0 : << LogIO::EXCEPTION;
1291 : }
1292 : }
1293 0 : }
1294 : else {
1295 : //no valid convolution for this pointing
1296 0 : for (Int chan=0; chan < convFunc_p.shape()[3]; ++chan) {
1297 0 : begin[3]=chan;
1298 0 : end[3]=chan;
1299 0 : convFunc_p(begin, end).set(Complex(0.0));
1300 0 : weightConvFunc_p(begin, end).set(Complex(0.0));
1301 : //convFunc_p.xyPlane(plane).set(0.0);
1302 : //weightConvFunc_p.xyPlane(plane).set(0.0);
1303 : }
1304 : }
1305 :
1306 0 : }
1307 :
1308 0 : void HetArrayConvFunc::supportAndNormalizeLatt(Int plane, Int convSampling, TempLattice<Complex>& convFuncLat,
1309 : TempLattice<Complex>& weightConvFuncLat) {
1310 :
1311 0 : LogIO os;
1312 0 : os << LogOrigin("HetArrConvFunc", "suppAndNorm") << LogIO::NORMAL;
1313 : // Locate support
1314 0 : Int convSupport=-1;
1315 : ///Use largest channel as highest freq thus largest conv func
1316 0 : IPosition begin(5, 0, 0, 0, convFuncLat.shape()(3)-1, plane);
1317 0 : IPosition shape(5, convFuncLat.shape()[0], convFuncLat.shape()[1], 1, 1, 1);
1318 : //Int convSize=convSize_p;
1319 0 : Int convSize=shape(0);
1320 : ///use FT weightconvlat as it is wider
1321 0 : Matrix<Complex> convPlane=weightConvFuncLat.getSlice(begin, shape, true);
1322 : Float maxAbsConvFunc, minAbsConvFunc;
1323 0 : IPosition minpos, maxpos;
1324 0 : minMax(minAbsConvFunc, maxAbsConvFunc, minpos, maxpos, amplitude(convPlane));
1325 0 : Bool found=false;
1326 0 : Int trial=0;
1327 0 : Float cutlevel=2.5e-2;
1328 : //numeric needs a larger ft
1329 0 : for (uInt k=0; k < antMath_p.nelements() ; ++k){
1330 0 : if((antMath_p[k]->whichPBClass()) == PBMathInterface::NUMERIC)
1331 0 : cutlevel=5e-3;
1332 : }
1333 :
1334 0 : for (trial=0; trial< (convSize-max(maxpos.asVector())-2); ++trial) {
1335 : ///largest along either axis
1336 : //cerr << "rat1 " << abs(convPlane(maxpos[0]-trial,maxpos[1]))/maxAbsConvFunc << " rat2 " << abs(convPlane(maxpos[0],maxpos[1]-trial))/maxAbsConvFunc << endl;
1337 0 : if((abs(convPlane(maxpos[0]-trial, maxpos[1])) < (cutlevel*maxAbsConvFunc)) &&(abs(convPlane(maxpos[0],maxpos[1]-trial)) < (cutlevel*maxAbsConvFunc)) )
1338 : {
1339 :
1340 0 : found=true;
1341 : //trial=Int(sqrt(2.0*Float(trial*trial)));
1342 :
1343 0 : break;
1344 : }
1345 : }
1346 0 : if(!found) {
1347 0 : if((maxAbsConvFunc-minAbsConvFunc) > (cutlevel*maxAbsConvFunc))
1348 0 : found=true;
1349 : // if it drops by more than 2 magnitudes per pixel
1350 : //trial=( (10*convSampling) < convSize) ? 5*convSampling : (convSize/2 - 4*convSampling);
1351 0 : trial=convSize/2 - 4*convSampling;
1352 : }
1353 :
1354 0 : if(found) {
1355 0 : if(trial < 5*convSampling)
1356 0 : trial= ( (10*convSampling) < convSize) ? 5*convSampling : (convSize/2 - 4*convSampling);
1357 0 : convSupport=(Int(0.5+Float(trial)/Float(convSampling)))+1 ;
1358 : //cerr << "convsupp " << convSupport << endl;
1359 : //support is really over the edge
1360 0 : if( (convSupport*convSampling) >= convSize/2) {
1361 0 : convSupport=convSize/2/convSampling-1;
1362 : }
1363 : }
1364 : else {
1365 : /*
1366 : os << "Convolution function is misbehaved - support seems to be zero\n"
1367 : << "Reasons can be: \nThe image definition not covering one or more of the pointings selected \n"
1368 : << "Or no unflagged data in a given pointing"
1369 :
1370 : << LogIO::EXCEPTION;
1371 : */
1372 : //OTF may have flagged stuff ...
1373 0 : convSupport=0;
1374 : }
1375 0 : convSupport_p(plane)=convSupport;
1376 0 : Double pbSum=0.0;
1377 : /*
1378 : Double pbSum1=0.0;
1379 :
1380 : for (Int iy=-convSupport;iy<=convSupport;iy++) {
1381 : for (Int ix=-convSupport;ix<=convSupport;ix++) {
1382 : Complex val=convFunc_p.xyPlane(plane)(ix*convSampling+convSize_p/2,
1383 : iy*convSampling+convSize_p/2);
1384 :
1385 : pbSum1+=sqrt(real(val)*real(val)+ imag(val)*imag(val));
1386 : }
1387 : }
1388 :
1389 : */
1390 : //cerr << "convSize_p " << convSize_p << " convSize " << convSize << endl;
1391 0 : if(convSupport >0) {
1392 0 : IPosition blc(2, -convSupport*convSampling+convSize/2, -convSupport*convSampling+convSize/2);
1393 0 : IPosition trc(2, convSupport*convSampling+convSize/2, convSupport*convSampling+convSize/2);
1394 0 : for (Int chan=0; chan < convFuncLat.shape()[3]; ++chan) {
1395 0 : begin[3]=chan;
1396 : //end[3]=chan;
1397 0 : convPlane.resize();
1398 0 : convPlane=convFuncLat.getSlice(begin, shape, true);
1399 0 : pbSum=real(sum(convPlane(blc,trc)))/Double(convSampling)/Double(convSampling);
1400 0 : if(pbSum>0.0) {
1401 0 : (convPlane)=convPlane*Complex(1.0/pbSum,0.0);
1402 0 : convFuncLat.putSlice(convPlane, begin);
1403 0 : convPlane.resize();
1404 0 : convPlane=weightConvFuncLat.getSlice(begin, shape, true);
1405 0 : Double pbSum1=0.0;
1406 0 : pbSum1=real(sum(convPlane(blc,trc)))/Double(convSampling)/Double(convSampling);
1407 0 : (convPlane) =(convPlane)*Complex(1.0/pbSum1,0.0);
1408 0 : weightConvFuncLat.putSlice(convPlane, begin);
1409 : }
1410 : else {
1411 : os << "Convolution function integral is not positive"
1412 0 : << LogIO::EXCEPTION;
1413 : }
1414 : }
1415 0 : }
1416 : else {
1417 : //no valid convolution for this pointing
1418 0 : for (Int chan=0; chan < convFuncLat.shape()[3]; ++chan) {
1419 0 : begin[3]=chan;
1420 : //end[3]=chan;
1421 0 : convPlane.resize(shape[0], shape[1]);
1422 0 : convPlane.set(Complex(0.0));
1423 0 : convFuncLat.putSlice(convPlane, begin);
1424 0 : weightConvFuncLat.putSlice(convPlane, begin);
1425 : //convFunc_p.xyPlane(plane).set(0.0);
1426 : //weightConvFunc_p.xyPlane(plane).set(0.0);
1427 : }
1428 : }
1429 :
1430 0 : }
1431 :
1432 0 : Int HetArrayConvFunc::factorial(Int n) {
1433 0 : Int fact=1;
1434 0 : for (Int k=1; k<=n; ++k)
1435 0 : fact *=k;
1436 0 : return fact;
1437 : }
1438 :
1439 :
1440 0 : Int HetArrayConvFunc::checkPBOfField(const vi::VisBuffer2& vb,
1441 : Vector<Int>& /*rowMap*/, const MVDirection& extraShift, const Bool useExtraShift) {
1442 :
1443 0 : toPix(vb, extraShift, useExtraShift);
1444 0 : Vector<Int> pixdepoint(2);
1445 0 : convertArray(pixdepoint, thePix_p);
1446 0 : if((pixdepoint(0) < 0) || pixdepoint(0) >= nx_p || pixdepoint(1) < 0 ||
1447 0 : pixdepoint(1) >=ny_p) {
1448 : //cout << "in pix de point off " << pixdepoint << endl;
1449 0 : return 2;
1450 : }
1451 0 : String pointingid=String::toString(pixdepoint(0))+"_"+String::toString(pixdepoint(1));
1452 0 : String msid=vb.msName(true);
1453 :
1454 :
1455 0 : if(convFunctionMap_p.nelements() == 0) {
1456 0 : convFunctionMap_p.resize(nx_p*ny_p);
1457 0 : convFunctionMap_p.set(-1);
1458 0 : convFunctionMap_p[pixdepoint[1]*nx_p+pixdepoint[0]]=0;
1459 0 : nDefined_p=1;
1460 0 : actualConvIndex_p=0;
1461 0 : return -1;
1462 : }
1463 :
1464 0 : if(convFunctionMap_p[pixdepoint[1]*nx_p+pixdepoint[0]] <0) {
1465 0 : actualConvIndex_p=nDefined_p;
1466 0 : convFunctionMap_p[pixdepoint[1]*nx_p+pixdepoint[0]]=nDefined_p;
1467 : // ++nDefined_p;
1468 0 : nDefined_p=1;
1469 0 : return -1;
1470 : }
1471 : else {
1472 0 : actualConvIndex_p=0;
1473 0 : return -1;
1474 : }
1475 :
1476 : return 1;
1477 :
1478 :
1479 0 : }
1480 :
1481 0 : void HetArrayConvFunc::makerowmap(const vi::VisBuffer2& vb,
1482 : Vector<Int>& rowMap) {
1483 :
1484 0 : uInt ndish=antMath_p.nelements();
1485 0 : rowMap.resize(vb.nRows());
1486 0 : for (rownr_t k=0; k < vb.nRows(); ++k) {
1487 0 : Int index1=antIndexToDiamIndex_p(vb.antenna1()(k));
1488 0 : Int index2=antIndexToDiamIndex_p(vb.antenna2()(k));
1489 0 : if(index2 < index1) {
1490 0 : index1=index2;
1491 0 : index2=antIndexToDiamIndex_p(vb.antenna1()(k));
1492 : }
1493 0 : Int plane=0;
1494 0 : for (Int jj=0; jj < index1; ++jj)
1495 0 : plane=plane+ndish-jj-1;
1496 0 : plane=plane+index2;
1497 : //plane of convfunc that match this pair of antennas
1498 0 : rowMap(k)=plane;
1499 :
1500 : }
1501 :
1502 0 : }
1503 :
1504 0 : Array<Complex> HetArrayConvFunc::resample(const Array<Complex>& inarray, const Double factor) {
1505 :
1506 0 : Double nx=Double(inarray.shape()(0));
1507 0 : Double ny=Double(inarray.shape()(1));
1508 0 : IPosition shp=inarray.shape();
1509 0 : shp(0)=Int(nx*factor/2.0)*2;
1510 0 : shp(1)=Int(ny*factor/2.0)*2;
1511 0 : Int newNx=shp(0);
1512 0 : Int newNy=shp(1);
1513 :
1514 0 : Array<Complex> out(shp, Complex(0.0));
1515 : // cerr << "SHP " << shp << endl;
1516 :
1517 0 : IPosition incursor=IPosition(inarray.shape().nelements(),1);
1518 0 : incursor[0]=nx;
1519 0 : incursor[1]=ny;
1520 0 : IPosition outcursor=IPosition(inarray.shape().nelements(),1);
1521 0 : outcursor[0]=shp[0];
1522 0 : outcursor[1]=shp[1];
1523 0 : ArrayIterator<Complex> inIt(inarray, IPosition(2,0,1), True);
1524 0 : ArrayIterator<Complex> outIt(out, IPosition(2,0,1),True);
1525 0 : inIt.origin();
1526 0 : outIt.origin();
1527 : //for (zzz=0; zzz< shp.(4); ++zzz){
1528 : // for(yyy=0; yyy< shp.(3); ++yyy){
1529 : // for(xxx=0; xxx< shp.(2); ++xxx){
1530 0 : while(!inIt.pastEnd()) {
1531 : // cerr << "Iter shape " << inIt.array().shape() << endl;
1532 0 : Matrix<Complex> inmat;
1533 0 : inmat=inIt.array();
1534 : //Matrix<Float> leReal=real(Matrix<Complex>(inIt.array()));
1535 : //Matrix<Float> leImag=imag(Matrix<Complex>(inIt.array()));
1536 0 : Matrix<Float> leReal=real(inmat);
1537 0 : Matrix<Float> leImag=imag(inmat);
1538 : Bool leRealCopy, leImagCopy;
1539 0 : Float *realptr=leReal.getStorage(leRealCopy);
1540 0 : Float *imagptr=leImag.getStorage(leImagCopy);
1541 : Bool isCopy;
1542 0 : Matrix<Complex> outMat(outIt.array());
1543 0 : Complex *intPtr=outMat.getStorage(isCopy);
1544 : Float realval, imagval;
1545 : #ifdef _OPENMP
1546 0 : omp_set_nested(0);
1547 : #endif
1548 0 : #pragma omp parallel for default(none) private(realval, imagval) firstprivate(intPtr, realptr, imagptr, nx, ny, newNx, newNy) shared(leReal, leImag)
1549 :
1550 : for (Int k =0; k < newNy; ++k) {
1551 : Double y =Double(k)/Double(newNy)*Double(ny);
1552 :
1553 : for (Int j=0; j < newNx; ++j) {
1554 : // Interpolate2D interp(Interpolate2D::LANCZOS);
1555 : Double x=Double(j)/Double(newNx)*Double(nx);
1556 : //interp.interp(realval, where, leReal);
1557 : realval=interpLanczos(x , y, nx, ny,
1558 : realptr);
1559 : imagval=interpLanczos(x , y, nx, ny,
1560 : imagptr);
1561 : //interp.interp(imagval, where, leImag);
1562 : intPtr[k*Int(newNx)+j]=Complex(realval, imagval);
1563 : }
1564 :
1565 : }
1566 0 : outMat.putStorage(intPtr, isCopy);
1567 0 : leReal.putStorage(realptr, leRealCopy);
1568 0 : leImag.putStorage(imagptr, leImagCopy);
1569 0 : inIt.next();
1570 0 : outIt.next();
1571 0 : }
1572 0 : return out;
1573 0 : }
1574 0 : Matrix <Complex> HetArrayConvFunc::resample2(const Matrix<Complex>& inarray, const Double factor) {
1575 :
1576 0 : Double nx=Double(inarray.shape()(0));
1577 0 : Double ny=Double(inarray.shape()(1));
1578 0 : IPosition shp=inarray.shape();
1579 0 : shp(0)=Int(nx*factor/2.0)*2;
1580 0 : shp(1)=Int(ny*factor/2.0)*2;
1581 :
1582 :
1583 0 : Matrix<Complex> outMat(shp, Complex(0.0));
1584 :
1585 :
1586 :
1587 : {
1588 : //cerr << "Iter shape " << inarray.shape() << endl;
1589 :
1590 0 : Matrix<Float> leReal=real(inarray);
1591 0 : Matrix<Float> leImag=imag(inarray);
1592 : Bool leRealCopy, leImagCopy;
1593 0 : Float *realptr=leReal.getStorage(leRealCopy);
1594 0 : Float *imagptr=leImag.getStorage(leImagCopy);
1595 : Bool isCopy;
1596 0 : Complex *intPtr=outMat.getStorage(isCopy);
1597 : Float realval, imagval;
1598 : #ifdef _OPENMP
1599 : // omp_set_nested(0);
1600 : #endif
1601 : // #pragma omp parallel for default(none) private(realval, imagval) firstprivate(intPtr, realptr, imagptr, nx, ny) shared(leReal, leImag)
1602 :
1603 0 : for (Int k =0; k < shp(1); ++k) {
1604 0 : Double y =Double(k)/Double(shp(1))*Double(ny);
1605 :
1606 0 : for (Int j=0; j < Int(nx*factor); ++j) {
1607 : // Interpolate2D interp(Interpolate2D::LANCZOS);
1608 0 : Double x=Double(j)/Double(factor);
1609 : //interp.interp(realval, where, leReal);
1610 0 : realval=interpLanczos(x , y, nx, ny,
1611 : realptr);
1612 0 : imagval=interpLanczos(x , y, nx, ny,
1613 : imagptr);
1614 : //interp.interp(imagval, where, leImag);
1615 0 : intPtr[k*Int(nx*factor)+j]=Complex(realval, imagval);
1616 : }
1617 :
1618 : }
1619 0 : outMat.putStorage(intPtr, isCopy);
1620 0 : leReal.putStorage(realptr, leRealCopy);
1621 0 : leImag.putStorage(imagptr, leImagCopy);
1622 :
1623 :
1624 0 : }
1625 0 : return outMat;
1626 0 : }
1627 0 : Float HetArrayConvFunc::sinc(const Float x) {
1628 0 : if (x == 0) {
1629 0 : return 1;
1630 : }
1631 0 : return sin(C::pi * x) / (C::pi * x);
1632 : }
1633 0 : Float HetArrayConvFunc::interpLanczos( const Double& x , const Double& y, const Double& nx, const Double& ny, const Float* data, const Float a) {
1634 0 : Double floorx = floor(x);
1635 0 : Double floory = floor(y);
1636 0 : Float result=0.0;
1637 0 : if (floorx < a || floorx >= nx - a || floory < a || floory >= ny - a) {
1638 0 : result = 0;
1639 0 : return result;
1640 : }
1641 0 : for (Float i = floorx - a + 1; i <= floorx + a; ++i) {
1642 0 : for (Float j = floory - a + 1; j <= floory + a; ++j) {
1643 0 : result += Float(Double(data[Int(j*nx+i)]) * sinc(x - i)*sinc((x-i)/ a) * sinc(y - j)*sinc((y-j)/ a));
1644 : }
1645 : }
1646 0 : return result;
1647 : }
1648 :
1649 0 : ImageInterface<Float>& HetArrayConvFunc::getFluxScaleImage() {
1650 0 : if(!calcFluxScale_p)
1651 0 : throw(AipsError("Programmer Error: flux image cannot be retrieved"));
1652 0 : if(!filledFluxScale_p) {
1653 : //The best flux image for a heterogenous array is the weighted coverage
1654 0 : fluxScale_p=TempImage<Float>(IPosition(4, nx_p, ny_p, npol_p, nchan_p), csys_p);
1655 0 : fluxScale_p.copyData(*(convWeightImage_p));
1656 0 : IPosition blc(4,nx_p, ny_p, npol_p, nchan_p);
1657 0 : IPosition trc(4, ny_p, ny_p, npol_p, nchan_p);
1658 0 : blc(0)=0;
1659 0 : blc(1)=0;
1660 0 : trc(0)=nx_p-1;
1661 0 : trc(1)=ny_p-1;
1662 :
1663 0 : for (Int j=0; j < npol_p; ++j) {
1664 0 : for (Int k=0; k < nchan_p ; ++k) {
1665 :
1666 0 : blc(2)=j;
1667 0 : trc(2)=j;
1668 0 : blc(3)=k;
1669 0 : trc(3)=k;
1670 0 : Slicer sl(blc, trc, Slicer::endIsLast);
1671 0 : SubImage<Float> fscalesub(fluxScale_p, sl, true);
1672 : Float planeMax;
1673 0 : LatticeExprNode LEN = max( fscalesub );
1674 0 : planeMax = LEN.getFloat();
1675 0 : if(planeMax !=0) {
1676 0 : fscalesub.copyData( (LatticeExpr<Float>) (fscalesub/planeMax));
1677 :
1678 : }
1679 0 : }
1680 : }
1681 0 : filledFluxScale_p=true;
1682 0 : }
1683 :
1684 :
1685 0 : return fluxScale_p;
1686 :
1687 : }
1688 :
1689 0 : void HetArrayConvFunc::sliceFluxScale(Int npol) {
1690 0 : IPosition fshp=fluxScale_p.shape();
1691 0 : if (fshp(2)>npol) {
1692 0 : npol_p=npol;
1693 : // use first npol planes...
1694 0 : IPosition blc(4,0,0,0,0);
1695 0 : IPosition trc(4,fluxScale_p.shape()(0)-1, fluxScale_p.shape()(1)-1,npol-1,fluxScale_p.shape()(3)-1);
1696 0 : Slicer sl=Slicer(blc, trc, Slicer::endIsLast);
1697 : //writeable if possible
1698 0 : SubImage<Float> fluxScaleSub = SubImage<Float> (fluxScale_p, sl, true);
1699 0 : SubImage<Float> convWeightImageSub = SubImage<Float> (*convWeightImage_p, sl, true);
1700 0 : fluxScale_p = TempImage<Float>(fluxScaleSub.shape(),fluxScaleSub.coordinates());
1701 0 : convWeightImage_p = new TempImage<Float> (convWeightImageSub.shape(),convWeightImageSub.coordinates());
1702 0 : LatticeExpr<Float> le(fluxScaleSub);
1703 0 : fluxScale_p.copyData(le);
1704 0 : LatticeExpr<Float> le2(convWeightImageSub);
1705 0 : convWeightImage_p->copyData(le2);
1706 0 : }
1707 0 : }
1708 : } // namespace refim end
1709 : } //# NAMESPACE CASA - END
1710 :
1711 :
1712 :
1713 :
|
