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