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 Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be adressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //#
27 : //# $Id$
28 :
29 : #include <casacore/casa/Arrays/ArrayMath.h>
30 : #include <casacore/casa/Arrays/ArrayLogical.h>
31 : #include <casacore/casa/Arrays/Array.h>
32 : #include <casacore/casa/Arrays/MaskedArray.h>
33 : #include <casacore/casa/Arrays/Vector.h>
34 : #include <casacore/casa/Arrays/Slice.h>
35 : #include <casacore/casa/Arrays/Matrix.h>
36 : #include <casacore/casa/Arrays/Cube.h>
37 : #include <casacore/scimath/Mathematics/FFTServer.h>
38 : #include <casacore/measures/Measures/MeasTable.h>
39 : #include <casacore/scimath/Mathematics/MathFunc.h>
40 : #include <casacore/scimath/Mathematics/ConvolveGridder.h>
41 : #include <casacore/casa/Utilities/Assert.h>
42 : #include <casacore/casa/Utilities/CompositeNumber.h>
43 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
44 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
45 :
46 : #include <casacore/images/Images/ImageInterface.h>
47 : #include <casacore/images/Images/PagedImage.h>
48 : #include <casacore/images/Images/SubImage.h>
49 : #include <casacore/images/Images/TempImage.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/VisBuffer.h>
66 : #include <msvis/MSVis/VisibilityIterator.h>
67 :
68 : #include <synthesis/TransformMachines/Utils.h>
69 : #include <synthesis/TransformMachines/PBMath1DAiry.h>
70 : #include <synthesis/TransformMachines/PBMath1DNumeric.h>
71 : #include <synthesis/TransformMachines/PBMath2DImage.h>
72 : #include <synthesis/TransformMachines/PBMath.h>
73 : #include <synthesis/TransformMachines/PBMathInterface.h>
74 : #include <synthesis/TransformMachines/HetArrayConvFunc.h>
75 : #include <synthesis/MeasurementEquations/VPManager.h>
76 :
77 : #include <casacore/casa/OS/Timer.h>
78 : #ifdef _OPENMP
79 : #include <omp.h>
80 : #endif
81 :
82 :
83 : using namespace casacore;
84 : namespace casa { //# NAMESPACE CASA - BEGIN
85 :
86 0 : HetArrayConvFunc::HetArrayConvFunc() : SimplePBConvFunc(), convFunctionMap_p(0), nDefined_p(0), msId_p(-1), actualConvIndex_p(-1)
87 : {
88 0 : calcFluxScale_p=true;
89 0 : init(PBMathInterface::AIRY);
90 0 : }
91 :
92 0 : HetArrayConvFunc::HetArrayConvFunc(const PBMathInterface::PBClass typeToUse, const String vpTable): SimplePBConvFunc(),
93 0 : convFunctionMap_p(0), nDefined_p(0), msId_p(-1), actualConvIndex_p(-1), vpTable_p(vpTable)
94 : {
95 0 : calcFluxScale_p=true;
96 0 : init(typeToUse);
97 :
98 0 : }
99 :
100 0 : HetArrayConvFunc::HetArrayConvFunc(const RecordInterface& rec, Bool calcFluxneeded): SimplePBConvFunc(), convFunctionMap_p(0), nDefined_p(0), msId_p(-1), actualConvIndex_p(-1) {
101 0 : String err;
102 0 : fromRecord(err, rec, calcFluxneeded);
103 0 : }
104 :
105 0 : HetArrayConvFunc::~HetArrayConvFunc(){
106 : //
107 0 : }
108 :
109 0 : void HetArrayConvFunc::init(const PBMathInterface::PBClass typeTouse){
110 0 : doneMainConv_p=false;
111 0 : filledFluxScale_p=false;
112 0 : pbClass_p=typeTouse;
113 0 : }
114 :
115 :
116 :
117 0 : void HetArrayConvFunc::findAntennaSizes(const VisBuffer& vb){
118 0 : if(msId_p != vb.msId()){
119 0 : msId_p=vb.msId();
120 :
121 0 : const MSAntennaColumns& ac=vb.msColumns().antenna();
122 0 : antIndexToDiamIndex_p.resize(ac.nrow());
123 0 : antIndexToDiamIndex_p.set(-1);
124 0 : Int diamIndex=antDiam2IndexMap_p.size( );
125 0 : Vector<Double> dishDiam=ac.dishDiameter().getColumn();
126 0 : Vector<String>dishName=ac.name().getColumn();
127 0 : String telescop=vb.msColumns().observation().telescopeName()(0);
128 : PBMath::CommonPB whichPB;
129 0 : if(pbClass_p==PBMathInterface::COMMONPB){
130 0 : String band;
131 0 : String commonPBName;
132 : // This frequency is ONLY required to determine which PB model to use:
133 : // The VLA, the ATNF, and WSRT have frequency - dependent PB models
134 0 : Quantity freq( vb.msColumns().spectralWindow().refFrequency()(0), "Hz");
135 :
136 :
137 0 : PBMath::whichCommonPBtoUse( telescop, freq, band, whichPB, commonPBName );
138 :
139 : ////////////////////
140 0 : String commonpbstr;
141 0 : PBMath::nameCommonPB(whichPB, commonpbstr);
142 : //cerr << "ms " << vb.msName() << " pb " << commonpbstr << endl;
143 :
144 :
145 : ////////////////
146 : //Revert to using AIRY for unknown common telescope
147 0 : if(whichPB==PBMath::UNKNOWN)
148 0 : pbClass_p=PBMathInterface::AIRY;
149 :
150 0 : }
151 0 : if(pbClass_p== PBMathInterface::AIRY){
152 0 : LogIO os;
153 0 : os << LogOrigin("HetArrConvFunc", "findAntennaSizes") << LogIO::NORMAL;
154 : ////////We'll be using dish diameter as key
155 0 : for (uInt k=0; k < dishDiam.nelements(); ++k){
156 0 : if((diamIndex !=0) && antDiam2IndexMap_p.find(String::toString(dishDiam(k))) != antDiam2IndexMap_p.end( )){
157 0 : antIndexToDiamIndex_p(k)=antDiam2IndexMap_p[String::toString(dishDiam(k))];
158 : }
159 : else{
160 0 : if(dishDiam[k] > 0.0){ //there may be stations with no dish on
161 0 : antDiam2IndexMap_p.insert(std::pair<casacore::String, casacore::Int>(String::toString(dishDiam(k)), diamIndex));
162 0 : antIndexToDiamIndex_p(k)=diamIndex;
163 0 : antMath_p.resize(diamIndex+1);
164 : //cerr << "diamindex " << diamIndex << " antMath_p,size " << antMath_p.nelements() << endl;
165 0 : if(pbClass_p== PBMathInterface::AIRY){
166 0 : Quantity qdiam= Quantity (dishDiam(k),"m");
167 : //ALMA ratio of blockage to dish
168 0 : Quantity blockDiam= Quantity(dishDiam(k)/12.0*.75, "m");
169 :
170 0 : if((vb.msColumns().observation().telescopeName()(0) =="ALMA") || (vb.msColumns().observation().telescopeName()(0) =="ACA")){
171 0 : if (abs(dishDiam[k] - 12.0) < 0.5) {
172 0 : qdiam= Quantity(10.7, "m");
173 0 : blockDiam= Quantity(0.75, "m");
174 : } else {
175 0 : qdiam= Quantity(6.25,"m");
176 0 : blockDiam = Quantity(0.75,"m");
177 : }
178 : }
179 0 : os << "Overriding PB with Airy of diam,blockage="<<qdiam<<","<<blockDiam<<" starting with antenna "<<k<<LogIO::POST;
180 0 : antMath_p[diamIndex]=new PBMath1DAiry(qdiam, blockDiam,
181 0 : Quantity(150,"arcsec"),
182 0 : Quantity(100.0,"GHz"));
183 :
184 :
185 0 : }
186 :
187 : //////Will no longer support this
188 : /*else if(pbClass_p== PBMathInterface::IMAGE){
189 : //Get the image name by calling code for the antenna name and array name
190 : //For now hard wired to ALMA as this part of the code will not be accessed for non-ALMA
191 : //see Imager::setMosaicFTMachine
192 : // When ready to generalize then code that calls with telescope name, antenna name
193 : //(via vb.msColumns) and/or diameter and frequency via vb.frequency (indexing will need to
194 : //be upgraded to account for frequency too) should be done to return the
195 : //right voltage pattern image.
196 : String vpImageName="";
197 : if (abs(dishDiam[k]-7.0) < 1.0)
198 : Aipsrc::find(vpImageName, "alma.vp.7m", "");
199 : else
200 : Aipsrc::find(vpImageName, "alma.vp.12m", "") ;
201 : //cerr << "first vpImagename " << vpImageName << endl;
202 : if(vpImageName==""){
203 : String beamPath;
204 : if(!MeasTable::AntennaResponsesPath(beamPath, "ALMA")){
205 : throw(AipsError("Alma beam images requested cannot be found "));
206 : }
207 : else{
208 : beamPath=beamPath.before(String("AntennaResponses"));
209 : vpImageName= (abs(dishDiam[k]-7.0) < 1.0) ? beamPath
210 : +String("/ALMA_AIRY_7M.VP") :
211 : beamPath+String("/ALMA_AIRY_12M.VP");
212 : }
213 :
214 :
215 : }
216 : //cerr << "Using the image VPs " << vpImageName << endl;
217 : if(Table::isReadable(vpImageName))
218 : antMath_p[diamIndex]=new PBMath2DImage(PagedImage<Complex>(vpImageName));
219 : else
220 : throw(AipsError(String("Cannot find voltage pattern image ") + vpImageName));
221 : }
222 : else{
223 :
224 : throw(AipsError("Do not deal with non airy dishes or images of VP yet "));
225 : }
226 : */
227 0 : ++diamIndex;
228 : }
229 : }
230 :
231 : }
232 :
233 :
234 0 : }
235 0 : else if(pbClass_p== PBMathInterface::IMAGE) {
236 :
237 0 : VPManager *vpman=VPManager::Instance();
238 0 : if(vpTable_p != String(""))
239 0 : vpman->loadfromtable(vpTable_p);
240 : ///else it is already loaded in the static object
241 0 : Vector<Record> recs;
242 0 : Vector<Vector<String> > antnames;
243 :
244 0 : if(vpman->imagepbinfo(antnames, recs)){
245 0 : Vector<Bool> dishDefined(dishName.nelements(), false);
246 0 : Int nbeams=antnames.nelements();
247 : ///will be keying on file image file name here
248 0 : for (uInt k=0; k < dishDiam.nelements(); ++k){
249 0 : String key;
250 0 : Bool beamDone=false;
251 0 : Int recordToUse=0;
252 0 : for (Int j =0; j < nbeams; ++j){
253 0 : key=recs[j].isDefined("realimage") ? recs[j].asString("realimage") : recs[j].asString("compleximage");
254 0 : if(antnames[j][0]=="*" || anyEQ(dishName[k], antnames[j])){
255 0 : dishDefined[k]=true;
256 0 : recordToUse=j;
257 :
258 0 : if((diamIndex !=0) && antDiam2IndexMap_p.find(key) != antDiam2IndexMap_p.end( )){
259 0 : antIndexToDiamIndex_p(k)=antDiam2IndexMap_p[key];
260 0 : beamDone=true;
261 : }
262 : }
263 : }
264 0 : if(!beamDone && dishDefined[k]){
265 0 : key=recs[recordToUse].isDefined("realimage") ? recs[recordToUse].asString("realimage") : recs[recordToUse].asString("compleximage");
266 0 : antDiam2IndexMap_p.insert(std::pair<casacore::String, casacore::Int>(key, diamIndex));
267 0 : antIndexToDiamIndex_p(k)=diamIndex;
268 0 : antMath_p.resize(diamIndex+1);
269 0 : if(recs[recordToUse].isDefined("realimage") && recs[recordToUse].isDefined("imagimage")){
270 0 : PagedImage<Float> realim(recs[recordToUse].asString("realimage"));
271 0 : PagedImage<Float> imagim(recs[recordToUse].asString("imagim"));
272 0 : antMath_p[diamIndex]=new PBMath2DImage(realim, imagim);
273 0 : }
274 : else {
275 0 : antMath_p[diamIndex]=new PBMath2DImage(PagedImage<Complex>(recs[recordToUse].asString("compleximage")));
276 : }
277 0 : ++diamIndex;
278 : }
279 0 : }
280 0 : if(!allTrue(dishDefined)){
281 : //cerr << "dishDefined " << dishDefined << endl;
282 0 : throw(AipsError("Some Antennas in the MS did not have a VP defined"));
283 : }
284 0 : }
285 : else{
286 0 : throw(AipsError("Mosaic does not support non-image voltage patterns yet"));
287 : }
288 :
289 : //Get rid of the static class
290 0 : vpman->reset();
291 0 : }
292 :
293 0 : else if(pbClass_p==PBMathInterface::COMMONPB){
294 : ///Have to use telescop part as string as in multims case different
295 : //telescopes may have same dish size but different commonpb
296 : // VLA and EVLA for e.g.
297 0 : if((diamIndex !=0) && antDiam2IndexMap_p.find(telescop+String("_")+String::toString(dishDiam(0))) != antDiam2IndexMap_p.end( )){
298 0 : antIndexToDiamIndex_p.set(antDiam2IndexMap_p[telescop+String("_")+String::toString(dishDiam(0))]);
299 : //cerr << " 0 telescop " << telescop << " index " << antDiam2IndexMap_p(telescop+String("_")+String::toString(dishDiam(0))) << endl;
300 : }
301 : else{
302 0 : antDiam2IndexMap_p.insert(std::pair<casacore::String, casacore::Int>(telescop+"_"+String::toString(dishDiam(0)), diamIndex));
303 0 : antIndexToDiamIndex_p.set(diamIndex);
304 0 : antMath_p.resize(diamIndex+1);
305 0 : antMath_p[diamIndex]=PBMath::pbMathInterfaceForCommonPB(whichPB, True);
306 : //cerr << " 1 telescop " << telescop << " index " << antDiam2IndexMap_p(telescop+String("_")+String::toString(dishDiam(0))) << endl;
307 :
308 : }
309 :
310 : }
311 : else{
312 :
313 0 : throw(AipsError("Mosaic supports image based or Airy voltage patterns only for now"));
314 :
315 : }
316 :
317 : //cerr << "antIndexTodiamIndex " << antIndexToDiamIndex_p << endl;
318 0 : }
319 :
320 :
321 :
322 :
323 :
324 0 : }
325 :
326 0 : void HetArrayConvFunc::reset(){
327 0 : doneMainConv_p=false;
328 0 : convFunctions_p.resize(0, true);
329 0 : convWeights_p.resize(0, true);
330 0 : convSizes_p.resize(0, true);
331 0 : convSupportBlock_p.resize(0, true);
332 0 : convFunctionMap_p.resize(0);
333 0 : vbConvIndex_p.clear();
334 :
335 0 : }
336 :
337 :
338 :
339 0 : void HetArrayConvFunc::findConvFunction(const ImageInterface<Complex>& iimage,
340 : const VisBuffer& vb,
341 : const Int& convSamp, const Vector<Double>& visFreq,
342 : Array<Complex>& convFunc,
343 : Array<Complex>& weightConvFunc,
344 : Vector<Int>& convsize,
345 : Vector<Int>& convSupport,
346 : Vector<Int>& convFuncPolMap,
347 : Vector<Int>& convFuncChanMap,
348 : Vector<Int>& convFuncRowMap,
349 : const Bool /*conjugateFreqFuncs*/)
350 : {
351 :
352 0 : storeImageParams(iimage,vb);
353 0 : convFuncChanMap.resize(vb.nChannel());
354 0 : Vector<Double> beamFreqs;
355 0 : findUsefulChannels(convFuncChanMap, beamFreqs, vb, visFreq);
356 0 : Int nBeamChans=beamFreqs.nelements();
357 : /////For now not doing beam rotation or squints but to be enabled easily
358 0 : convFuncPolMap.resize(vb.nCorr());
359 0 : Int nBeamPols=1;
360 0 : convFuncPolMap.set(0);
361 0 : Bool msChanged= msId_p != vb.msId();
362 0 : findAntennaSizes(vb);
363 0 : uInt ndish=antMath_p.nelements();
364 0 : if(ndish==0)
365 0 : throw(AipsError("Don't have dishsize"));
366 : Int ndishpair;
367 0 : if(ndish==1)
368 0 : ndishpair=1;
369 : else
370 0 : ndishpair=factorial(ndish)/factorial(ndish-2)/2 + ndish;
371 :
372 0 : convFunc.resize();
373 0 : weightConvFunc.resize();
374 0 : convFuncRowMap.resize();
375 0 : convsize.resize();
376 0 : convSupport.resize();
377 :
378 0 : Int isCached=checkPBOfField(vb, convFuncRowMap);
379 : //cout << "isCached " << isCached << endl;
380 0 : if(isCached==1 && (convFuncRowMap.shape()[0]==vb.nRow())){
381 : /*convFunc.reference(convFunc_p);
382 : weightConvFunc.reference(weightConvFunc_p);
383 : convsize=*convSizes_p[actualConvIndex_p];
384 : convSupport=convSupport_p;
385 : return;
386 : */
387 : }
388 0 : else if(isCached ==2){
389 0 : convFunc.resize();
390 0 : weightConvFunc.resize();
391 0 : convsize.resize();
392 0 : convSupport.resize();
393 0 : convFuncRowMap.resize();
394 0 : return;
395 :
396 : }
397 0 : actualConvIndex_p=convIndex(vb);
398 0 : if(doneMainConv_p.shape()[0] < (actualConvIndex_p+1)){
399 : //cerr << "resizing DONEMAIN " << doneMainConv_p.shape()[0] << endl;
400 0 : doneMainConv_p.resize(actualConvIndex_p+1, true);
401 0 : doneMainConv_p[actualConvIndex_p]=false;
402 0 : convFunctions_p.resize(actualConvIndex_p+1);
403 0 : convFunctions_p[actualConvIndex_p]=nullptr;
404 : }
405 : ///// In multi ms mode ndishpair may change when meeting a new ms
406 : //// redo the calculation then
407 0 : if(msChanged){
408 : //cerr << "ms " << vb.msName() << " spw " << vb.spectralWindow() << endl;
409 0 : doneMainConv_p[actualConvIndex_p]=false;
410 : //cerr << "invalidating doneMainConv " << endl; //<< convFunctions_p[actualConvIndex_p]->shape()[3] << " =? " << nBeamChans << " convsupp " << convSupport_p.nelements() << endl;
411 : }
412 :
413 : ////Trap for cases when the selection seem to have changed
414 0 : if(doneMainConv_p[actualConvIndex_p]){
415 0 : if(nBeamChans != (*convFunctions_p[actualConvIndex_p]).shape()[3])
416 0 : doneMainConv_p[actualConvIndex_p]=False;
417 :
418 : }
419 :
420 : // Get the coordinate system
421 :
422 0 : CoordinateSystem coords;
423 0 : coords=iimage.coordinates();
424 0 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
425 0 : AlwaysAssert(directionIndex>=0, AipsError);
426 : // Set up the convolution function.
427 0 : Int nx=nx_p;
428 0 : Int ny=ny_p;
429 0 : Int support=0;
430 0 : Int convSampling=1;
431 0 : if(!doneMainConv_p[actualConvIndex_p]){
432 0 : for (uInt ii=0; ii < ndish; ++ii){
433 0 : support=max((antMath_p[ii])->support(coords), support);
434 : }
435 0 : support=Int(min(Float(support), max(Float(nx_p), Float(ny_p)))*2.0)/2;
436 0 : convSize_p=support*convSampling;
437 : // Make this a nice composite number, to speed up FFTs
438 0 : CompositeNumber cn(uInt(convSize_p*2.0));
439 0 : convSize_p = cn.nextLargerEven(Int(convSize_p));
440 0 : convSize_p=(convSize_p/16)*16; // need it to be divisible by 8 in places
441 0 : }
442 :
443 :
444 0 : DirectionCoordinate dc=dc_p;
445 : //where in the image in pixels is this pointing
446 0 : Vector<Double> pixFieldDir(2);
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 : }
452 : //no need to call toPix here as its been done already above in checkOFPB
453 : //thus the values are still current.
454 0 : pixFieldDir=thePix_p;
455 : //toPix(pixFieldDir, vb);
456 0 : MDirection fieldDir=direction1_p;
457 : //shift from center
458 0 : pixFieldDir(0)=pixFieldDir(0)- Double(nx / 2);
459 0 : pixFieldDir(1)=pixFieldDir(1)- Double(ny / 2);
460 : //phase gradient per pixel to apply
461 0 : pixFieldDir(0)=-pixFieldDir(0)*2.0*C::pi/Double(nx)/Double(convSamp);
462 0 : pixFieldDir(1)=-pixFieldDir(1)*2.0*C::pi/Double(ny)/Double(convSamp);
463 :
464 :
465 :
466 0 : if(!doneMainConv_p[actualConvIndex_p]){
467 0 : Vector<Double> sampling;
468 0 : sampling = dc.increment();
469 0 : sampling*=Double(convSampling);
470 0 : sampling(0)*=Double(nx)/Double(convSize_p);
471 0 : sampling(1)*=Double(ny)/Double(convSize_p);
472 0 : dc.setIncrement(sampling);
473 :
474 0 : Vector<Double> unitVec(2);
475 0 : unitVec=convSize_p/2;
476 0 : dc.setReferencePixel(unitVec);
477 : //make sure we are using the same units
478 0 : fieldDir.set(dc.worldAxisUnits()(0));
479 0 : dc.setReferenceValue(fieldDir.getAngle().getValue());
480 0 : coords.replaceCoordinate(dc, directionIndex);
481 0 : Int spind=coords.findCoordinate(Coordinate::SPECTRAL);
482 0 : SpectralCoordinate spCoord=coords.spectralCoordinate(spind);
483 0 : spCoord.setReferencePixel(Vector<Double>(1,0.0));
484 0 : spCoord.setReferenceValue(Vector<Double>(1, beamFreqs(0)));
485 0 : if(beamFreqs.nelements() >1)
486 0 : spCoord.setIncrement(Vector<Double>(1, beamFreqs(1)-beamFreqs(0)));
487 :
488 0 : coords.replaceCoordinate(spCoord, spind);
489 :
490 0 : IPosition pbShape(4, convSize_p, convSize_p, 1, nBeamChans);
491 : //TempImage<Complex> twoDPB(pbShape, coords);
492 :
493 :
494 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);
495 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);
496 : //convFunc_p.resize(IPosition(5, convSize_p, convSize_p, nBeamPols, nBeamChans, ndishpair));
497 :
498 : // convFunc_p=0.0;
499 : //weightConvFunc_p.resize(IPosition(5, convSize_p, convSize_p, nBeamPols, nBeamChans, ndishpair));
500 : //weightConvFunc_p=0.0;
501 0 : IPosition begin(5, 0, 0, 0, 0, 0);
502 0 : IPosition end(5, convFuncTemp.shape()[0]-1, convFuncTemp.shape()[1]-1, nBeamPols-1, nBeamChans-1, 0);
503 : //FFTServer<Float, Complex> fft(IPosition(2, convSize_p, convSize_p));
504 : //TempImage<Complex> pBScreen(TiledShape(pbShape, IPosition(4, convSize_p, convSize_p, 1, 1)), coords, 0);
505 : //TempImage<Complex> pB2Screen(TiledShape(pbShape, IPosition(4, convSize_p, convSize_p, 1, 1)), coords, 0);
506 0 : Long memtot=HostInfo::memoryFree();
507 0 : Double memtobeused= Double(memtot)*1024.0;
508 0 : if(memtot <= 2000000)
509 0 : memtobeused=0.0;
510 : //cerr << "Shape " << pbShape << " memtobe used " << memtobeused << endl;
511 0 : TempImage<Complex> pBScreen(TiledShape(pbShape), coords, memtobeused/2.2);
512 0 : TempImage<Complex> pB2Screen(TiledShape(pbShape), coords, memtobeused/2.2);
513 0 : IPosition start(4, 0, 0, 0, 0);
514 0 : convSupport_p.resize(ndishpair);
515 :
516 0 : for (uInt k=0; k < ndish; ++k){
517 :
518 0 : for (uInt j =k ; j < ndish; ++j){
519 : //Timer tim;
520 : //Matrix<Complex> screen(convSize_p, convSize_p);
521 : //screen=1.0;
522 : //pBScreen.putSlice(screen, start);
523 : //cerr << "k " << k << " shape " << pBScreen.shape() << " direction1 " << direction1_p << " direction2 " << direction2_p << endl;
524 :
525 : //pBScreen.set(Complex(1.0, 0.0));
526 : //one antenna
527 0 : IPosition blcin(4, 0, 0, 0, 0);
528 0 : IPosition trcin(4, convSize_p-1, convSize_p-1, 0, 0);
529 0 : for (Int kk=0; kk < nBeamChans; ++kk){
530 0 : blcin[3]=kk;
531 0 : trcin[3]=kk;
532 : // tim.mark();
533 : // Double wtime0=omp_get_wtime();
534 : //cerr << "Doing channel " << kk << endl;
535 0 : Slicer slin(blcin, trcin, Slicer::endIsLast);
536 0 : SubImage<Complex> subim(pBScreen, slin, true);
537 0 : subim.set(Complex(1.0, 0.0));
538 0 : (antMath_p[k])->applyVP(subim, subim, direction1_p);
539 :
540 : //Then the other
541 0 : (antMath_p[j])->applyVP(subim, subim, direction2_p);
542 : //tim.show("After Apply ");
543 : //tim.mark();
544 : //pB2Screen.set(Complex(1.0,0.0));
545 0 : SubImage<Complex> subim2(pB2Screen, slin, true);
546 0 : subim2.set(Complex(1.0,0.0));
547 : //one antenna
548 0 : (antMath_p[k])->applyPB(subim2, subim2, direction1_p);
549 : //Then the other
550 0 : (antMath_p[j])->applyPB(subim2, subim2, direction2_p);
551 :
552 : //Double wtime1=omp_get_wtime();
553 : //tim.show("After Apply2 ");
554 : //tim.mark();
555 : //subim.copyData((LatticeExpr<Complex>) (iif(abs(subim)> 5e-2, subim, 0)));
556 : //subim2.copyData((LatticeExpr<Complex>) (iif(abs(subim2)> 25e-4, subim2, 0)));
557 :
558 :
559 : //ft_p.c2cFFT(subim);
560 :
561 : //ft_p.c2cFFT(subim2);
562 0 : ft_p.c2cFFTInDouble(subim);
563 0 : ft_p.c2cFFTInDouble(subim2);
564 : //LatticeFFT::cfft2d(subim2);
565 : // tim.show("after ffts ");
566 : // cerr << kk << " applies " << wtime1-wtime0 << " ffts " << omp_get_wtime()-wtime1 << endl;
567 :
568 :
569 0 : }
570 : /*
571 : if(nBeamChans >1){
572 : Slicer slin(blcin, trcin, Slicer::endIsLast);
573 : SubImage<Complex> origPB(pBScreen, slin, false);
574 : IPosition elshape= origPB.shape();
575 : Matrix<Complex> i1=origPB.get(true);
576 : SubImage<Complex> origPB2(pB2Screen, slin, false);
577 : Matrix<Complex> i2=origPB2.get(true);
578 : Int cenX=i1.shape()(0)/2;
579 : Int cenY=i1.shape()(1)/2;
580 :
581 : for (Int kk=0; kk < (nBeamChans-1); ++kk){
582 : Double fratio=beamFreqs(kk)/beamFreqs(nBeamChans-1);
583 : cerr << "fratio " << fratio << endl;
584 : blcin[3]=kk;
585 : trcin[3]=kk;
586 : //Slicer slout(blcin, trcin, Slicer::endIsLast);
587 : Matrix<Complex> o1(i1.shape(), Complex(0.0));
588 : Matrix<Complex> o2(i2.shape(), Complex(0.0));
589 : for (Int yy=0; yy < i1.shape()(1); ++yy){
590 : //Int nyy= (Double(yy-cenY)*fratio) + cenY;
591 : Double nyy= (Double(yy-cenY)/fratio) + cenY;
592 : Double cyy=ceil(nyy);
593 : Double fyy= floor(nyy);
594 : Int iy=nyy > fyy+0.5 ? cyy : fyy;
595 : if(cyy <2*cenY && fyy >=0.0)
596 : for(Int xx=0; xx < i1.shape()(0); ++ xx){
597 : //Int nxx= Int(Double(xx-cenX)*fratio) + cenX;
598 : Double nxx= Int(Double(xx-cenX)/fratio) + cenX;
599 : Double cxx=ceil(nxx);
600 : Double fxx= floor(nxx);
601 : Int ix=nxx > fxx+0.5 ? cxx : fxx ;
602 : if(cxx < 2*cenX && fxx >=0.0 ){
603 : //Double dist=sqrt((nxx-cxx)*(nxx-cxx)+(nyy-cyy)*(nyy-cyy))/sqrt(2.0);
604 : //o1(xx, yy)=float(1-dist)*i1(fxx, fyy)+ dist*i1(cxx,cyy);
605 : o1(xx, yy)=i1( ix, iy);
606 : //o2(xx, yy)=i2(nxx, nyy);
607 : //o2(xx, yy)=float(1-dist)*i2(fxx, fyy)+ dist*i2(cxx,cyy);
608 : o2(xx, yy)=i2(ix, iy);
609 : }
610 : }
611 : }
612 : pBScreen.putSlice(o1.reform(elshape), blcin);
613 : pB2Screen.putSlice(o2.reform(elshape), blcin);
614 : }
615 :
616 : }
617 : */
618 :
619 : //tim.show("after apply+apply2+masking+fft ");
620 : //tim.mark();
621 : //LatticeFFT::cfft2d(pBScreen);
622 : //LatticeFFT::cfft2d(pB2Screen);
623 :
624 : //Matrix<Complex> lala=pBScreen.get(true);
625 : //fft.fft0(lala, true);
626 : //fft.flip(lala, false, false);
627 : // pBScreen.put(lala.reform(IPosition(4, convSize_p, convSize_p, 1, 1)));
628 : //lala=pB2Screen.get(true);
629 : //fft.fft0(lala, true);
630 : //fft.flip(lala, false, false);
631 : //pB2Screen.put(lala.reform(IPosition(4, convSize_p, convSize_p, 1, 1)));
632 :
633 : //////////*****************
634 : /*if(0){
635 : ostringstream os1;
636 : os1 << "PB_field_" << Int(thePix_p[0]) << "_" << Int(thePix_p[1]) << "_antpair_" << k <<"_"<<j ;
637 : PagedImage<Float> thisScreen(pbShape, coords, String(os1));
638 : LatticeExpr<Float> le(abs(pBScreen));
639 : thisScreen.copyData(le);
640 :
641 : }*/
642 : ////////*****************/
643 :
644 : //tim.show("after FFT ");
645 : //tim.mark();
646 0 : Int plane=0;
647 0 : for (uInt jj=0; jj < k; ++jj)
648 0 : plane=plane+ndish-jj-1;
649 0 : plane=plane+j;
650 0 : begin[4]=plane;
651 0 : end[4]=plane;
652 : //cerr << "ms" << vb.msName() << " spw " << vb.spectralWindow() << " k " << k << " j " << j << " plane " << plane << endl;
653 0 : Slicer slplane(begin, end, Slicer::endIsLast);
654 : //cerr << "SHAPES " << convFunc_p(begin, end).shape() << " " << pBScreen.get(false).shape() << " begin and end " << begin << " " << end << endl;
655 : //convFunc_p(begin, end).copyMatchingPart(pBScreen.get(false));
656 : //weightConvFunc_p(begin, end).copyMatchingPart(pB2Screen.get(false));
657 0 : IPosition blcQ(4, pbShape(0)/8*3, pbShape(1)/8*3, 0, 0);
658 0 : IPosition trcQ(4, blcQ[0]+ pbShape(0)/4-1, blcQ[1]+pbShape(1)/4-1 , nBeamPols-1, nBeamChans-1);
659 :
660 : //cerr << "blcQ " << blcQ << " trcQ " << trcQ << " pbShape " << pbShape << endl;
661 0 : Slicer slQ(blcQ, trcQ, Slicer::endIsLast);
662 : {
663 0 : SubImage<Complex> pBSSub(pBScreen, slQ, false);
664 0 : SubLattice<Complex> cFTempSub(convFuncTemp, slplane, true);
665 0 : LatticeConcat<Complex> lc(4);
666 0 : lc.setLattice(pBSSub);
667 : //cerr << "SHAPES " << cFTempSub.shape() << " " << lc.shape() << endl;
668 0 : cFTempSub.copyData(lc);
669 : //cFTempSub.copyData(pBScreen);
670 0 : }
671 : {
672 0 : SubImage<Complex> pB2SSub(pB2Screen, slQ, false);
673 0 : SubLattice<Complex> cFTempSub2(weightConvFuncTemp, slplane, true);
674 0 : LatticeConcat<Complex> lc(4);
675 0 : lc.setLattice(pB2SSub);
676 0 : cFTempSub2.copyData(lc);
677 : // cFTempSub2.copyData(pB2Screen);
678 : //weightConvFuncTemp.putSlice(pB2Screen.get(false), begin);
679 :
680 0 : }
681 : // supportAndNormalize(plane, convSampling);
682 0 : supportAndNormalizeLatt( plane, convSampling, convFuncTemp, weightConvFuncTemp);
683 :
684 :
685 :
686 : // tim.show("After search of support ");
687 0 : }
688 :
689 : }
690 :
691 :
692 0 : doneMainConv_p[actualConvIndex_p]=true;
693 0 : convFunctions_p.resize(actualConvIndex_p+1);
694 0 : convWeights_p.resize(actualConvIndex_p+1);
695 0 : convSupportBlock_p.resize(actualConvIndex_p+1);
696 :
697 0 : convFunctions_p[actualConvIndex_p]= new Array<Complex>();
698 0 : convWeights_p[actualConvIndex_p]= new Array<Complex>();
699 0 : convSupportBlock_p[actualConvIndex_p]=new Vector<Int>();
700 0 : Int newConvSize=2*(max(convSupport_p)+2)*convSampling;
701 0 : Int newRealConvSize=newConvSize* Int(Double(convSamp)/Double(convSampling));
702 0 : Int lattSize=convFuncTemp.shape()(0);
703 0 : (*convSupportBlock_p[actualConvIndex_p])=convSupport_p;
704 :
705 : //cerr << "convSize " << newConvSize << " " << newRealConvSize<< " lattSize " << lattSize << endl;
706 0 : LogIO os(LogOrigin("HetArrConvFunc", "findConvFunction", WHERE));
707 0 : os << "convolution function support: " << convSupport_p << LogIO::POST;
708 :
709 0 : if(newConvSize < lattSize){
710 0 : IPosition blc(5, (lattSize/2)-(newConvSize/2),
711 0 : (lattSize/2)-(newConvSize/2),0,0,0);
712 0 : IPosition trc(5, (lattSize/2)+(newConvSize/2-1),
713 0 : (lattSize/2)+(newConvSize/2-1), nBeamPols-1, nBeamChans-1,ndishpair-1);
714 0 : IPosition shp(5, newConvSize, newConvSize, nBeamPols, nBeamChans, ndishpair);
715 :
716 0 : convFunctions_p[actualConvIndex_p]= new Array<Complex>(IPosition(5, newRealConvSize, newRealConvSize, nBeamPols, nBeamChans, ndishpair ));
717 0 : convWeights_p[actualConvIndex_p]= new Array<Complex>(IPosition(5, newRealConvSize, newRealConvSize, nBeamPols, nBeamChans, ndishpair ));
718 0 : (*convFunctions_p[actualConvIndex_p])=resample(convFuncTemp.getSlice(blc,shp), Double(convSamp)/Double(convSampling));
719 0 : convSize_p=newRealConvSize;
720 0 : (*convWeights_p[actualConvIndex_p])=resample(weightConvFuncTemp.getSlice(blc, shp),Double(convSamp)/Double(convSampling)) ;
721 0 : }
722 : else{
723 0 : newRealConvSize=lattSize* Int(Double(convSamp)/Double(convSampling));
724 0 : convFunctions_p[actualConvIndex_p]= new Array<Complex>(IPosition(5, newRealConvSize, newRealConvSize, nBeamPols, nBeamChans, ndishpair ));
725 0 : convWeights_p[actualConvIndex_p]= new Array<Complex>(IPosition(5, newRealConvSize, newRealConvSize, nBeamPols, nBeamChans, ndishpair ));
726 :
727 0 : (*convFunctions_p[actualConvIndex_p])=resample(convFuncTemp.get(), Double(convSamp)/Double(convSampling));
728 0 : (*convWeights_p[actualConvIndex_p])=resample(weightConvFuncTemp.get(), Double(convSamp)/Double(convSampling));
729 0 : convSize_p=newRealConvSize;
730 : }
731 0 : convFunc_p.resize();
732 0 : weightConvFunc_p.resize();
733 :
734 0 : }
735 : else{
736 0 : convSize_p=max(*convSizes_p[actualConvIndex_p]);
737 0 : convSupport_p.resize();
738 0 : convSupport_p=*convSupportBlock_p[actualConvIndex_p];
739 : }
740 : /*
741 : rowMap.resize(vb.nRow());
742 : for (Int k=0; k < vb.nRow(); ++k){
743 : //plane of convfunc that match this pair of antennas
744 : rowMap(k)=antIndexToDiamIndex_p(vb.antenna1()(k))*ndish+
745 : antIndexToDiamIndex_p(vb.antenna2()(k));
746 :
747 : }
748 : */
749 :
750 :
751 0 : makerowmap(vb, convFuncRowMap);
752 : ///need to deal with only the maximum of different baselines available in this
753 : ///vb
754 : //cerr << "ms " << vb.msName() << " convFuncRowMap " << convFuncRowMap[0] << endl;
755 0 : ndishpair=max(convFuncRowMap)+1;
756 :
757 : // convSupportBlock_p.resize(actualConvIndex_p+1);
758 0 : convSizes_p.resize(actualConvIndex_p+1);
759 : //convSupportBlock_p[actualConvIndex_p]=new Vector<Int>(ndishpair);
760 : //(*convSupportBlock_p[actualConvIndex_p])=convSupport_p;
761 0 : convSizes_p[actualConvIndex_p]=new Vector<Int> (ndishpair);
762 :
763 : /* convFunctions_p[actualConvIndex_p]->resize(convSize_p, convSize_p, ndishpair);
764 : *(convFunctions_p[actualConvIndex_p])=convSave_p;
765 : convWeights_p[actualConvIndex_p]->resize(convSize_p, convSize_p, ndishpair);
766 : *(convWeights_p[actualConvIndex_p])=weightSave_p;
767 : */
768 :
769 0 : convFunc_p.resize();
770 0 : convFunc_p=(*convFunctions_p[actualConvIndex_p]);
771 0 : weightConvFunc_p.resize();
772 0 : weightConvFunc_p=(*convWeights_p[actualConvIndex_p]);
773 :
774 : //cerr << "convfunc shapes " << convFunc_p.shape() << " " << weightConvFunc_p.shape() << " " << convSize_p << " pol " << nBeamPols << " chan " << nBeamChans << " ndishpair " << ndishpair << endl;
775 : /////Due to a bug in buildCoordSysCore...sometimes an image bigger
776 : ///than the spw selection chosen is made
777 0 : if(nBeamChans > convFunc_p.shape()[3])
778 0 : nBeamChans = convFunc_p.shape()[3];
779 : //convSupport_p.resize();
780 : //convSupport_p=(*convSupportBlock_p[actualConvIndex_p]);
781 : Bool delc; Bool delw;
782 0 : Double dirX=pixFieldDir(0);
783 0 : Double dirY=pixFieldDir(1);
784 0 : Complex *convstor=convFunc_p.getStorage(delc);
785 0 : Complex *weightstor=weightConvFunc_p.getStorage(delw);
786 0 : Int elconvsize=convSize_p;
787 :
788 0 : #pragma omp parallel default(none) firstprivate(convstor, weightstor, dirX, dirY, elconvsize, ndishpair, nBeamChans, nBeamPols)
789 : {
790 : #pragma omp for
791 : for(Int iy=0; iy<elconvsize; ++iy){
792 : applyGradientToYLine(iy, convstor, weightstor, dirX, dirY, elconvsize, ndishpair, nBeamChans, nBeamPols);
793 :
794 : }
795 : }///End of pragma
796 :
797 0 : convFunc_p.putStorage(convstor, delc);
798 0 : weightConvFunc_p.putStorage(weightstor, delw);
799 :
800 :
801 :
802 : //For now all have the same size convsize;
803 0 : convSizes_p[actualConvIndex_p]->set(convSize_p);
804 :
805 : //We have to get the references right now
806 : // convFunc_p.resize();
807 : //convFunc_p.reference(*convFunctions_p[actualConvIndex_p]);
808 : //weightConvFunc_p.resize();
809 : //weightConvFunc_p.reference(*convWeights_p[actualConvIndex_p]);
810 :
811 0 : convFunc.reference(convFunc_p);
812 0 : weightConvFunc.reference(weightConvFunc_p);
813 0 : convsize=*convSizes_p[actualConvIndex_p];
814 0 : convSupport=convSupport_p;
815 :
816 :
817 0 : }
818 :
819 : typedef unsigned long long ooLong;
820 :
821 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){
822 : Double cy, sy;
823 :
824 0 : SINCOS(Double(iy-convSize/2)*pixYdir, sy, cy);
825 0 : Complex phy(cy,sy) ;
826 0 : for (Int ix=0;ix<convSize;ix++) {
827 : Double cx, sx;
828 0 : SINCOS(Double(ix-convSize/2)*pixXdir, sx, cx);
829 0 : Complex phx(cx,sx) ;
830 0 : for (Int ipol=0; ipol< nPol; ++ipol){
831 : //Int poloffset=ipol*nChan*ndishpair*convSize*convSize;
832 0 : for (Int ichan=0; ichan < nChan; ++ichan){
833 : //Int chanoffset=ichan*ndishpair*convSize*convSize;
834 0 : for(Int iz=0; iz <ndishpair; ++iz){
835 0 : ooLong index=((ooLong(iz*nChan+ichan)*nPol+ipol)*ooLong(convSize)+ooLong(iy))*ooLong(convSize)+ooLong(ix);
836 0 : convFunctions[index]= convFunctions[index]*phx*phy;
837 0 : convWeights[index]= convWeights[index]*phx*phy;
838 : }
839 : }
840 : }
841 :
842 : }
843 0 : }
844 0 : Bool HetArrayConvFunc::toRecord(RecordInterface& rec){
845 :
846 : try{
847 0 : rec.define("name", "HetArrayConvFunc");
848 0 : Int numConv=convFunctions_p.nelements();
849 0 : rec.define("ndefined", numConv);
850 : //rec.define("convfunctionmap", convFunctionMap_p);
851 0 : std::map<String, Int>::iterator it=vbConvIndex_p.begin();
852 0 : for (Int64 k=0; k < numConv; ++k){
853 0 : rec.define("convfunctions"+String::toString(k), *(convFunctions_p[k]));
854 0 : rec.define("convweights"+String::toString(k), *(convWeights_p[k]));
855 0 : rec.define("convsizes"+String::toString(k), *(convSizes_p[k]));
856 0 : rec.define("convsupportblock"+String::toString(k), *(convSupportBlock_p[k]));
857 0 : rec.define(String("key")+String::toString(k), it->first);
858 0 : rec.define(String("val")+String::toString(k), it->second);
859 0 : it++;
860 : }
861 0 : rec.define("actualconvindex", actualConvIndex_p);
862 0 : rec.define("donemainconv", doneMainConv_p);
863 0 : rec.define("vptable", vpTable_p);
864 0 : rec.define("pbclass", Int(pbClass_p));
865 :
866 : }
867 0 : catch(AipsError& x) {
868 0 : return false;
869 0 : }
870 0 : return true;
871 :
872 : }
873 :
874 0 : Bool HetArrayConvFunc::fromRecord(String& err, const RecordInterface& rec, Bool calcfluxscale){
875 : //Force pbmath stuff and saved image stuff
876 0 : nchan_p=0;
877 0 : msId_p=-1;
878 : try{
879 0 : if(!rec.isDefined("name") || rec.asString("name") != "HetArrayConvFunc"){
880 0 : throw(AipsError("Wrong record to recover HetArray from"));
881 : }
882 0 : nDefined_p=rec.asInt("ndefined");
883 : //rec.get("convfunctionmap", convFunctionMap_p);
884 0 : convFunctions_p.resize(nDefined_p, true, false);
885 0 : convSupportBlock_p.resize(nDefined_p, true, false);
886 0 : convWeights_p.resize(nDefined_p, true, false);
887 0 : convSizes_p.resize(nDefined_p, true, false);
888 0 : vbConvIndex_p.erase(vbConvIndex_p.begin(), vbConvIndex_p.end());
889 0 : for (Int64 k=0; k < nDefined_p; ++k){
890 0 : convFunctions_p[k]=new Array<Complex>();
891 0 : convWeights_p[k]=new Array<Complex>();
892 0 : convSizes_p[k]=new Vector<Int>();
893 0 : convSupportBlock_p[k]=new Vector<Int>();
894 0 : rec.get("convfunctions"+String::toString(k), *(convFunctions_p[k]));
895 0 : rec.get("convweights"+String::toString(k), *(convWeights_p[k]));
896 0 : rec.get("convsizes"+String::toString(k), *(convSizes_p[k]));
897 0 : rec.get("convsupportblock"+String::toString(k), *(convSupportBlock_p[k]));
898 0 : String key;
899 : Int val;
900 0 : rec.get(String("key")+String::toString(k), key);
901 0 : rec.get(String("val")+String::toString(k), val);
902 0 : vbConvIndex_p[key]=val;
903 0 : }
904 : //Now that we are calculating all phase gradients on the fly ..
905 : ///we should clean up some and get rid of the cached variables
906 :
907 0 : convSize_p= nDefined_p > 0 ? (*(convSizes_p[0]))[0] : 0;
908 : //convSave_p.resize();
909 : //rec.get("convsave", convSave_p);
910 : //weightSave_p.resize();
911 : //rec.get("weightsave", weightSave_p);
912 0 : rec.get("vptable", vpTable_p);
913 0 : rec.get("donemainconv", doneMainConv_p);
914 : //convSupport_p.resize();
915 : //rec.get("convsupport", convSupport_p);
916 0 : pbClass_p=static_cast<PBMathInterface::PBClass>(rec.asInt("pbclass"));
917 0 : calcFluxScale_p=calcfluxscale;
918 : }
919 0 : catch(AipsError& x) {
920 0 : err=x.getMesg();
921 0 : return false;
922 0 : }
923 :
924 0 : return true;
925 : }
926 :
927 :
928 0 : void HetArrayConvFunc::supportAndNormalize(Int plane, Int convSampling){
929 :
930 0 : LogIO os;
931 0 : os << LogOrigin("HetArrConvFunc", "suppAndNorm") << LogIO::NORMAL;
932 : // Locate support
933 0 : Int convSupport=-1;
934 0 : IPosition begin(5, 0, 0, 0, 0, plane);
935 0 : IPosition end(5, convFunc_p.shape()[0]-1, convFunc_p.shape()[1]-1, 0, 0, plane);
936 0 : Matrix<Complex> convPlane(convFunc_p(begin, end).reform(IPosition(2,convFunc_p.shape()[0], convFunc_p.shape()[1]))) ;
937 0 : Float maxAbsConvFunc=max(amplitude(convPlane));
938 0 : Float minAbsConvFunc=min(amplitude(convPlane));
939 0 : Bool found=false;
940 0 : Int trial=0;
941 0 : for (trial=convSize_p/2-2;trial>0;trial--) {
942 : //Searching down a diagonal
943 0 : if(abs(convPlane(convSize_p/2-trial,convSize_p/2-trial)) > (1.0e-2*maxAbsConvFunc)) {
944 0 : found=true;
945 0 : trial=Int(sqrt(2.0*Float(trial*trial)));
946 0 : break;
947 : }
948 : }
949 0 : if(!found){
950 0 : if((maxAbsConvFunc-minAbsConvFunc) > (1.0e-2*maxAbsConvFunc))
951 0 : found=true;
952 : // if it drops by more than 2 magnitudes per pixel
953 0 : trial=( (10*convSampling) < convSize_p) ? 5*convSampling : (convSize_p/2 - 4*convSampling);
954 : }
955 :
956 :
957 0 : if(found) {
958 0 : if(trial < 5*convSampling)
959 0 : trial= ( (10*convSampling) < convSize_p) ? 5*convSampling : (convSize_p/2 - 4*convSampling);
960 0 : convSupport=Int(0.5+Float(trial)/Float(convSampling))+1;
961 : //support is really over the edge
962 0 : if( (convSupport*convSampling) >= convSize_p/2){
963 0 : convSupport=convSize_p/2/convSampling-1;
964 : }
965 : }
966 : else {
967 : /*
968 : os << "Convolution function is misbehaved - support seems to be zero\n"
969 : << "Reasons can be: \nThe image definition not covering one or more of the pointings selected \n"
970 : << "Or no unflagged data in a given pointing"
971 :
972 : << LogIO::EXCEPTION;
973 : */
974 : //OTF may have flagged stuff ...
975 0 : convSupport=0;
976 : }
977 : //cerr << "trial " << trial << " convSupport " << convSupport << " convSize " << convSize_p << endl;
978 0 : convSupport_p(plane)=convSupport;
979 0 : Double pbSum=0.0;
980 : /*
981 : Double pbSum1=0.0;
982 :
983 : for (Int iy=-convSupport;iy<=convSupport;iy++) {
984 : for (Int ix=-convSupport;ix<=convSupport;ix++) {
985 : Complex val=convFunc_p.xyPlane(plane)(ix*convSampling+convSize_p/2,
986 : iy*convSampling+convSize_p/2);
987 :
988 : pbSum1+=sqrt(real(val)*real(val)+ imag(val)*imag(val));
989 : }
990 : }
991 :
992 : */
993 0 : if(convSupport >0){
994 0 : IPosition blc(2, -convSupport*convSampling+convSize_p/2, -convSupport*convSampling+convSize_p/2);
995 0 : IPosition trc(2, convSupport*convSampling+convSize_p/2, convSupport*convSampling+convSize_p/2);
996 0 : for (Int chan=0; chan < convFunc_p.shape()[3]; ++chan){
997 0 : begin[3]=chan;
998 0 : end[3]=chan;
999 0 : convPlane.resize();
1000 0 : convPlane.reference(convFunc_p(begin, end).reform(IPosition(2,convFunc_p.shape()[0], convFunc_p.shape()[1])));
1001 0 : pbSum=real(sum(convPlane(blc,trc)))/Double(convSampling)/Double(convSampling);
1002 0 : if(pbSum>0.0) {
1003 0 : (convPlane)=convPlane*Complex(1.0/pbSum,0.0);
1004 0 : convPlane.resize();
1005 0 : convPlane.reference(weightConvFunc_p(begin, end).reform(IPosition(2,convFunc_p.shape()[0], convFunc_p.shape()[1])));
1006 0 : (convPlane) =(convPlane)*Complex(1.0/pbSum,0.0);
1007 : }
1008 : else {
1009 : os << "Convolution function integral is not positive"
1010 0 : << LogIO::EXCEPTION;
1011 : }
1012 : }
1013 0 : }
1014 : else{
1015 : //no valid convolution for this pointing
1016 0 : for (Int chan=0; chan < convFunc_p.shape()[3]; ++chan){
1017 0 : begin[3]=chan;
1018 0 : end[3]=chan;
1019 0 : convFunc_p(begin, end).set(Complex(0.0));
1020 0 : weightConvFunc_p(begin, end).set(Complex(0.0));
1021 : //convFunc_p.xyPlane(plane).set(0.0);
1022 : //weightConvFunc_p.xyPlane(plane).set(0.0);
1023 : }
1024 : }
1025 :
1026 0 : }
1027 :
1028 0 : void HetArrayConvFunc::supportAndNormalizeLatt(Int plane, Int convSampling, TempLattice<Complex>& convFuncLat,
1029 : TempLattice<Complex>& weightConvFuncLat){
1030 :
1031 0 : LogIO os;
1032 0 : os << LogOrigin("HetArrConvFunc", "suppAndNorm") << LogIO::NORMAL;
1033 : // Locate support
1034 0 : Int convSupport=-1;
1035 : ///Use largest channel as highest freq thus largest conv func
1036 0 : IPosition begin(5, 0, 0, 0, convFuncLat.shape()(3)-1, plane);
1037 0 : IPosition shape(5, convFuncLat.shape()[0], convFuncLat.shape()[1], 1, 1, 1);
1038 : //Int convSize=convSize_p;
1039 0 : Int convSize=shape(0);
1040 : ///use FT weightconvlat as it is wider
1041 0 : Matrix<Complex> convPlane=weightConvFuncLat.getSlice(begin, shape, true);
1042 0 : Float maxAbsConvFunc=max(amplitude(convPlane));
1043 0 : Float minAbsConvFunc=min(amplitude(convPlane));
1044 0 : Bool found=false;
1045 0 : Int trial=0;
1046 0 : Float cutlevel=2.5e-2;
1047 : //numeric needs a larger ft
1048 0 : for (uInt k=0; k < antMath_p.nelements() ; ++k){
1049 0 : if((antMath_p[k]->whichPBClass()) == PBMathInterface::NUMERIC)
1050 0 : cutlevel=1e-3;
1051 : }
1052 0 : for (trial=convSize/2-2;trial>0;trial--) {
1053 : //largest of either
1054 0 : if((abs(convPlane(convSize/2-trial-1,convSize/2-1)) > (cutlevel*maxAbsConvFunc)) || (abs(convPlane(convSize/2-1,convSize/2-trial-1)) > (cutlevel*maxAbsConvFunc))) {
1055 0 : found=true;
1056 0 : trial=Int(sqrt(2.0*Float(trial*trial)));
1057 0 : break;
1058 : }
1059 : }
1060 0 : if(!found){
1061 0 : if((maxAbsConvFunc-minAbsConvFunc) > (cutlevel*maxAbsConvFunc))
1062 0 : found=true;
1063 : // if it drops by more than 2 magnitudes per pixel
1064 0 : trial=( (10*convSampling) < convSize) ? 5*convSampling : (convSize/2 - 4*convSampling);
1065 : }
1066 :
1067 :
1068 0 : if(found) {
1069 0 : if(trial < 5*convSampling)
1070 0 : trial= ( (10*convSampling) < convSize) ? 5*convSampling : (convSize/2 - 4*convSampling);
1071 0 : convSupport=Int(0.5+Float(trial)/Float(convSampling))+1;
1072 : //support is really over the edge
1073 0 : if( (convSupport*convSampling) >= convSize/2){
1074 0 : convSupport=convSize/2/convSampling-1;
1075 : }
1076 : }
1077 : else {
1078 : /*
1079 : os << "Convolution function is misbehaved - support seems to be zero\n"
1080 : << "Reasons can be: \nThe image definition not covering one or more of the pointings selected \n"
1081 : << "Or no unflagged data in a given pointing"
1082 :
1083 : << LogIO::EXCEPTION;
1084 : */
1085 : //OTF may have flagged stuff ...
1086 0 : convSupport=0;
1087 : }
1088 : //cerr << "trial " << trial << " convSupport " << convSupport << " convSize " << convSize_p << endl;
1089 0 : convSupport_p(plane)=convSupport;
1090 0 : Double pbSum=0.0;
1091 : /*
1092 : Double pbSum1=0.0;
1093 :
1094 : for (Int iy=-convSupport;iy<=convSupport;iy++) {
1095 : for (Int ix=-convSupport;ix<=convSupport;ix++) {
1096 : Complex val=convFunc_p.xyPlane(plane)(ix*convSampling+convSize_p/2,
1097 : iy*convSampling+convSize_p/2);
1098 :
1099 : pbSum1+=sqrt(real(val)*real(val)+ imag(val)*imag(val));
1100 : }
1101 : }
1102 :
1103 : */
1104 : //cerr << "convSize_p " << convSize_p << " convSize " << convSize << endl;
1105 0 : if(convSupport >0){
1106 0 : IPosition blc(2, -convSupport*convSampling+convSize/2, -convSupport*convSampling+convSize/2);
1107 0 : IPosition trc(2, convSupport*convSampling+convSize/2, convSupport*convSampling+convSize/2);
1108 0 : for (Int chan=0; chan < convFuncLat.shape()[3]; ++chan){
1109 0 : begin[3]=chan;
1110 : //end[3]=chan;
1111 0 : convPlane.resize();
1112 0 : convPlane=convFuncLat.getSlice(begin, shape, true);
1113 : ///////
1114 : //Matrix<Float> ampi=amplitude(convPlane) ;
1115 : //LogicalArray mask(ampi < Float(1.e-3* maxAbsConvFunc));
1116 :
1117 : //convPlane(mask)=0.0;
1118 : /////
1119 0 : pbSum=real(sum(convPlane(blc,trc)))/Double(convSampling)/Double(convSampling);
1120 0 : if(pbSum>0.0) {
1121 0 : (convPlane)=convPlane*Complex(1.0/pbSum,0.0);
1122 0 : convFuncLat.putSlice(convPlane, begin);
1123 : //convPlane.resize();
1124 0 : Matrix<Complex> convPlane2;
1125 0 : convPlane2=weightConvFuncLat.getSlice(begin, shape, true);
1126 : //convPlane2(mask)=0.0;
1127 0 : (convPlane2) =(convPlane2)*Complex(1.0/pbSum,0.0);
1128 0 : weightConvFuncLat.putSlice(convPlane2, begin);
1129 0 : }
1130 : else {
1131 : os << "Convolution function integral is not positive"
1132 0 : << LogIO::EXCEPTION;
1133 : }
1134 : }
1135 0 : }
1136 : else{
1137 : //no valid convolution for this pointing
1138 0 : for (Int chan=0; chan < convFuncLat.shape()[3]; ++chan){
1139 0 : begin[3]=chan;
1140 : //end[3]=chan;
1141 0 : convPlane.resize(shape[0], shape[1]);
1142 0 : convPlane.set(Complex(0.0));
1143 0 : convFuncLat.putSlice(convPlane, begin);
1144 0 : weightConvFuncLat.putSlice(convPlane, begin);
1145 : //convFunc_p.xyPlane(plane).set(0.0);
1146 : //weightConvFunc_p.xyPlane(plane).set(0.0);
1147 : }
1148 : }
1149 :
1150 0 : }
1151 :
1152 0 : Int HetArrayConvFunc::factorial(Int n){
1153 0 : Int fact=1;
1154 0 : for (Int k=1; k<=n; ++k)
1155 0 : fact *=k;
1156 0 : return fact;
1157 : }
1158 :
1159 :
1160 0 : Int HetArrayConvFunc::checkPBOfField(const VisBuffer& vb,
1161 : Vector<Int>& /*rowMap*/){
1162 :
1163 0 : toPix(vb);
1164 0 : Vector<Int> pixdepoint(2);
1165 0 : convertArray(pixdepoint, thePix_p);
1166 0 : if((pixdepoint(0) < 0) || pixdepoint(0) >= nx_p || pixdepoint(1) < 0 ||
1167 0 : pixdepoint(1) >=ny_p){
1168 : //cout << "in pix de point off " << pixdepoint << endl;
1169 0 : return 2;
1170 : }
1171 0 : String pointingid=String::toString(pixdepoint(0))+"_"+String::toString(pixdepoint(1));
1172 : //Int fieldid=vb.fieldId();
1173 0 : String msid=vb.msName(true);
1174 : //If channel or pol length has changed underneath...then its time to
1175 : //restart the map
1176 : /*
1177 : if(convFunctionMap_p.ndefined() > 0){
1178 : if ((fluxScale_p.shape()[3] != nchan_p) || (fluxScale_p.shape()[2] != npol_p)){
1179 : convFunctionMap_p.clear();
1180 : }
1181 : }
1182 :
1183 : */
1184 0 : if(convFunctionMap_p.nelements() > 0){
1185 0 : if (calcFluxScale_p && ((fluxScale_p.shape()[3] != nchan_p) || (fluxScale_p.shape()[2] != npol_p))){
1186 0 : convFunctionMap_p.resize();
1187 0 : nDefined_p=0;
1188 : }
1189 : }
1190 : //String mapid=msid+String("_")+pointingid;
1191 : /*
1192 : if(convFunctionMap_p.ndefined() == 0){
1193 : convFunctionMap_p.define(mapid, 0);
1194 : actualConvIndex_p=0;
1195 : fluxScale_p=TempImage<Float>(IPosition(4,nx_p,ny_p,npol_p,nchan_p), csys_p);
1196 : filledFluxScale_p=false;
1197 : fluxScale_p.set(0.0);
1198 : return -1;
1199 : }
1200 : */
1201 0 : if(convFunctionMap_p.nelements() == 0){
1202 0 : convFunctionMap_p.resize(nx_p*ny_p);
1203 0 : convFunctionMap_p.set(-1);
1204 0 : convFunctionMap_p[pixdepoint[1]*nx_p+pixdepoint[0]]=0;
1205 0 : nDefined_p=1;
1206 0 : actualConvIndex_p=0;
1207 0 : if(calcFluxScale_p){
1208 0 : fluxScale_p=TempImage<Float>(IPosition(4,nx_p,ny_p,npol_p,nchan_p), csys_p);
1209 0 : filledFluxScale_p=false;
1210 0 : fluxScale_p.set(0.0);
1211 : }
1212 0 : return -1;
1213 : }
1214 :
1215 : // if(!convFunctionMap_p.isDefined(mapid)){
1216 : // actualConvIndex_p=convFunctionMap_p.ndefined();
1217 : // convFunctionMap_p.define(mapid, actualConvIndex_p);
1218 0 : if(convFunctionMap_p[pixdepoint[1]*nx_p+pixdepoint[0]] <0){
1219 0 : actualConvIndex_p=nDefined_p;
1220 0 : convFunctionMap_p[pixdepoint[1]*nx_p+pixdepoint[0]]=nDefined_p;
1221 : // ++nDefined_p;
1222 0 : nDefined_p=1;
1223 0 : return -1;
1224 : }
1225 : else{
1226 : /*
1227 : actualConvIndex_p=convFunctionMap_p[pixdepoint[1]*nx_p+pixdepoint[0]];
1228 : convFunc_p.resize(); // break any reference
1229 : weightConvFunc_p.resize();
1230 : convSupport_p.resize();
1231 : //Here we will need to use the right xyPlane for different PA range
1232 : //and frequency may be
1233 : convFunc_p.reference(*convFunctions_p[actualConvIndex_p]);
1234 : weightConvFunc_p.reference(*convWeights_p[actualConvIndex_p]);
1235 : //Again this for one time of antenna only later should be fixed for all
1236 : // antennas independently
1237 : //these are not really needed right now
1238 : convSupport_p=(*convSupportBlock_p[actualConvIndex_p]);
1239 : convSize_p=(*convSizes_p[actualConvIndex_p])[0];
1240 : makerowmap(vb, rowMap);
1241 : */
1242 0 : actualConvIndex_p=0;
1243 0 : return -1;
1244 : }
1245 :
1246 : return 1;
1247 :
1248 :
1249 0 : }
1250 :
1251 0 : void HetArrayConvFunc::makerowmap(const VisBuffer& vb,
1252 : Vector<Int>& rowMap){
1253 :
1254 0 : uInt ndish=antMath_p.nelements();
1255 0 : rowMap.resize(vb.nRow());
1256 0 : for (Int k=0; k < vb.nRow(); ++k){
1257 0 : Int index1=antIndexToDiamIndex_p(vb.antenna1()(k));
1258 0 : Int index2=antIndexToDiamIndex_p(vb.antenna2()(k));
1259 0 : if(index2 < index1){
1260 0 : index1=index2;
1261 0 : index2=antIndexToDiamIndex_p(vb.antenna1()(k));
1262 : }
1263 0 : Int plane=0;
1264 0 : for (Int jj=0; jj < index1; ++jj)
1265 0 : plane=plane+ndish-jj-1;
1266 0 : plane=plane+index2;
1267 : //plane of convfunc that match this pair of antennas
1268 0 : rowMap(k)=plane;
1269 :
1270 : }
1271 :
1272 0 : }
1273 :
1274 0 : Array<Complex> HetArrayConvFunc::resample(const Array<Complex>& inarray, const Double factor){
1275 :
1276 0 : Double nx=Double(inarray.shape()(0));
1277 0 : Double ny=Double(inarray.shape()(1));
1278 0 : IPosition shp=inarray.shape();
1279 0 : shp(0)=Int(nx*factor);
1280 0 : shp(1)=Int(ny*factor);
1281 0 : cerr << "nx " << nx << " ny " << ny << " shape " << shp << endl;
1282 :
1283 0 : Array<Complex> out(shp);
1284 0 : ArrayIterator<Complex> inIt(inarray, 2);
1285 0 : ArrayIterator<Complex> outIt(out, 2);
1286 : //for (zzz=0; zzz< shp.(4); ++zzz){
1287 : // for(yyy=0; yyy< shp.(3); ++yyy){
1288 : // for(xxx=0; xxx< shp.(2); ++xxx){
1289 0 : while(!inIt.pastEnd()){
1290 0 : Matrix<Float> leReal=real(Matrix<Complex>(inIt.array()));
1291 0 : Matrix<Float> leImag=imag(Matrix<Complex>(inIt.array()));
1292 : Bool leRealCopy, leImagCopy;
1293 0 : Float *realptr=leReal.getStorage(leRealCopy);
1294 0 : Float *imagptr=leImag.getStorage(leImagCopy);
1295 : Bool isCopy;
1296 0 : Matrix<Complex> outMat(outIt.array());
1297 0 : Complex *intPtr=outMat.getStorage(isCopy);
1298 : Float realval, imagval;
1299 : #ifdef _OPENMP
1300 0 : omp_set_nested(0);
1301 : #endif
1302 0 : #pragma omp parallel for default(none) private(realval, imagval) firstprivate(intPtr, realptr, imagptr, nx, ny, factor) shared(leReal, leImag)
1303 :
1304 : for (Int k =0; k < Int(ny*factor); ++k){
1305 : Double y =Double(k)/factor;
1306 :
1307 : for (Int j=0; j < Int(nx*factor); ++j){
1308 : // Interpolate2D interp(Interpolate2D::LANCZOS);
1309 : Double x=Double(j)/Double(factor);
1310 : //interp.interp(realval, where, leReal);
1311 : realval=interpLanczos(x , y, nx, ny,
1312 : realptr);
1313 : imagval=interpLanczos(x , y, nx, ny,
1314 : imagptr);
1315 : //interp.interp(imagval, where, leImag);
1316 : intPtr[k*Int(nx*factor)+j]=Complex(realval, imagval);
1317 : }
1318 :
1319 : }
1320 0 : outMat.putStorage(intPtr, isCopy);
1321 0 : leReal.putStorage(realptr, leRealCopy);
1322 0 : leImag.putStorage(imagptr, leImagCopy);
1323 0 : inIt.next();
1324 0 : outIt.next();
1325 0 : }
1326 0 : return out;
1327 0 : }
1328 :
1329 0 : Float HetArrayConvFunc::sinc(const Float x) {
1330 0 : if (x == 0) {
1331 0 : return 1;
1332 : }
1333 0 : return sin(C::pi * x) / (C::pi * x);
1334 : }
1335 0 : Float HetArrayConvFunc::interpLanczos( const Double& x , const Double& y, const Double& nx, const Double& ny, const Float* data, const Float a){
1336 0 : Double floorx = floor(x);
1337 0 : Double floory = floor(y);
1338 0 : Float result=0.0;
1339 0 : if (floorx < a || floorx >= nx - a || floory < a || floory >= ny - a) {
1340 0 : result = 0;
1341 0 : return result;
1342 : }
1343 0 : for (Float i = floorx - a + 1; i <= floorx + a; ++i) {
1344 0 : for (Float j = floory - a + 1; j <= floory + a; ++j) {
1345 0 : result += Float(Double(data[Int(j*nx+i)]) * sinc(x - i)*sinc((x-i)/ a) * sinc(y - j)*sinc((y-j)/ a));
1346 : }
1347 : }
1348 0 : return result;
1349 : }
1350 :
1351 0 : ImageInterface<Float>& HetArrayConvFunc::getFluxScaleImage(){
1352 0 : if(!calcFluxScale_p)
1353 0 : throw(AipsError("Programmer Error: flux image cannot be retrieved"));
1354 0 : if(!filledFluxScale_p){
1355 : //The best flux image for a heterogenous array is the weighted coverage
1356 0 : fluxScale_p.copyData(*(convWeightImage_p));
1357 0 : IPosition blc(4,nx_p, ny_p, npol_p, nchan_p);
1358 0 : IPosition trc(4, ny_p, ny_p, npol_p, nchan_p);
1359 0 : blc(0)=0; blc(1)=0; trc(0)=nx_p-1; trc(1)=ny_p-1;
1360 :
1361 0 : for (Int j=0; j < npol_p; ++j){
1362 0 : for (Int k=0; k < nchan_p ; ++k){
1363 :
1364 0 : blc(2)=j; trc(2)=j;
1365 0 : blc(3)=k; trc(3)=k;
1366 0 : Slicer sl(blc, trc, Slicer::endIsLast);
1367 0 : SubImage<Float> fscalesub(fluxScale_p, sl, true);
1368 : Float planeMax;
1369 0 : LatticeExprNode LEN = max( fscalesub );
1370 0 : planeMax = LEN.getFloat();
1371 0 : if(planeMax !=0){
1372 0 : fscalesub.copyData( (LatticeExpr<Float>) (fscalesub/planeMax));
1373 :
1374 : }
1375 0 : }
1376 : }
1377 0 : filledFluxScale_p=true;
1378 0 : }
1379 :
1380 :
1381 0 : return fluxScale_p;
1382 :
1383 : }
1384 :
1385 0 : void HetArrayConvFunc::sliceFluxScale(Int npol) {
1386 0 : IPosition fshp=fluxScale_p.shape();
1387 0 : if (fshp(2)>npol){
1388 0 : npol_p=npol;
1389 : // use first npol planes...
1390 0 : IPosition blc(4,0,0,0,0);
1391 0 : IPosition trc(4,fluxScale_p.shape()(0)-1, fluxScale_p.shape()(1)-1,npol-1,fluxScale_p.shape()(3)-1);
1392 0 : Slicer sl=Slicer(blc, trc, Slicer::endIsLast);
1393 : //writeable if possible
1394 0 : SubImage<Float> fluxScaleSub = SubImage<Float> (fluxScale_p, sl, true);
1395 0 : SubImage<Float> convWeightImageSub = SubImage<Float> (*convWeightImage_p, sl, true);
1396 0 : fluxScale_p = TempImage<Float>(fluxScaleSub.shape(),fluxScaleSub.coordinates());
1397 0 : convWeightImage_p = new TempImage<Float> (convWeightImageSub.shape(),convWeightImageSub.coordinates());
1398 0 : LatticeExpr<Float> le(fluxScaleSub);
1399 0 : fluxScale_p.copyData(le);
1400 0 : LatticeExpr<Float> le2(convWeightImageSub);
1401 0 : convWeightImage_p->copyData(le2);
1402 0 : }
1403 0 : }
1404 :
1405 : } //# NAMESPACE CASA - END
1406 :
1407 :
1408 :
1409 :
|
