Line data Source code
1 : // -*- C++ -*-
2 : //# EVLAConvFunc.cc: Implementation of the EVLAConvFunc class
3 : //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
4 : //# Associated Universities, Inc. Washington DC, USA.
5 : //#
6 : //# This library is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU Library General Public License as published by
8 : //# the Free Software Foundation; either version 2 of the License, or (at your
9 : //# option) any later version.
10 : //#
11 : //# This library is distributed in the hope that it will be useful, but WITHOUT
12 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
14 : //# License for more details.
15 : //#
16 : //# You should have received a copy of the GNU Library General Public License
17 : //# along with this library; if not, write to the Free Software Foundation,
18 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning AIPS++ should be addressed as follows:
21 : //# Internet email: aips2-request@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : //#
27 : //# $Id$
28 : //
29 : #include <synthesis/TransformMachines/EVLAConvFunc.h>
30 : #include <synthesis/TransformMachines/SynthesisError.h>
31 : #include <synthesis/TransformMachines/BeamCalc.h>
32 : #include <synthesis/TransformMachines/WTerm.h>
33 : //
34 : //---------------------------------------------------------------------
35 : //---------------------------------------------------------------------
36 : // TEMPS
37 : //---------------------------------------------------------------------
38 : //---------------------------------------------------------------------
39 : //
40 : #define CONVSIZE (1024*2)
41 : #define CONVWTSIZEFACTOR 1.0
42 : #define OVERSAMPLING 20
43 : #define THRESHOLD 1E-4
44 :
45 : using namespace casacore;
46 : namespace casa{
47 :
48 0 : EVLAConvFunc& EVLAConvFunc::operator=(const EVLAConvFunc& other)
49 : {
50 0 : if(this!=&other)
51 : {
52 0 : ConvolutionFunction::operator=(other);
53 0 : logIO_p = other.logIO_p;
54 : // setParams(other.polMap_p, other.feedStokes_p);
55 0 : setPolMap(other.polMap_p);
56 0 : bandID_p=other.bandID_p;
57 0 : Diameter_p=other.Diameter_p;
58 0 : Nant_p=other.Nant_p;
59 0 : HPBW=other.HPBW;
60 0 : sigma=other.sigma;
61 : }
62 0 : return *this;
63 : }
64 0 : Int EVLAConvFunc::getVLABandID(Double& freq,String&telescopeName)
65 : {
66 0 : if (telescopeName=="VLA")
67 : {
68 0 : if ((freq >=1.34E9) && (freq <=1.73E9))
69 0 : return BeamCalc_VLA_L;
70 0 : else if ((freq >=4.5E9) && (freq <=5.0E9))
71 0 : return BeamCalc_VLA_C;
72 0 : else if ((freq >=8.0E9) && (freq <=8.8E9))
73 0 : return BeamCalc_VLA_X;
74 0 : else if ((freq >=14.4E9) && (freq <=15.4E9))
75 0 : return BeamCalc_VLA_U;
76 0 : else if ((freq >=22.0E9) && (freq <=24.0E9))
77 0 : return BeamCalc_VLA_K;
78 0 : else if ((freq >=40.0E9) && (freq <=50.0E9))
79 0 : return BeamCalc_VLA_Q;
80 0 : else if ((freq >=100E6) && (freq <=300E6))
81 0 : return BeamCalc_VLA_4;
82 : }
83 : else
84 0 : if (telescopeName=="EVLA")
85 : {
86 0 : if ((freq >=0.6E9) && (freq <=2.0E9))
87 0 : return BeamCalc_EVLA_L;
88 0 : else if ((freq >=2.0E9) && (freq <=4.0E9))
89 0 : return BeamCalc_EVLA_S;
90 0 : else if ((freq >=4.0E9) && (freq <=8.0E9))
91 0 : return BeamCalc_EVLA_C;
92 0 : else if ((freq >=8.0E9) && (freq <=12.0E9))
93 0 : return BeamCalc_EVLA_X;
94 0 : else if ((freq >=12.0E9) && (freq <=18.0E9))
95 0 : return BeamCalc_EVLA_U;
96 0 : else if ((freq >=18.0E9) && (freq <=26.5E9))
97 0 : return BeamCalc_EVLA_K;
98 0 : else if ((freq >=26.5E9) && (freq <=40.8E9))
99 0 : return BeamCalc_EVLA_A;
100 0 : else if ((freq >=4.0E9) && (freq <=50.0E9))
101 0 : return BeamCalc_EVLA_Q;
102 : }
103 0 : ostringstream mesg;
104 0 : mesg << telescopeName << "/" << freq << "(Hz) combination not recognized.";
105 0 : throw(SynthesisError(mesg.str()));
106 : }
107 :
108 0 : void EVLAConvFunc::setPolMap(const Vector<Int>& polMap)
109 0 : {polMap_p.resize(0);polMap_p=polMap;};
110 0 : void EVLAConvFunc::setFeedStokes(const Vector<Int>& feedStokes)
111 0 : {feedStokes_p.resize(0);feedStokes_p=feedStokes;};
112 :
113 0 : int EVLAConvFunc::getVisParams(const VisBuffer& vb)
114 : {
115 : Double Freq;
116 :
117 0 : Vector<String> telescopeNames=vb.msColumns().observation().telescopeName().getColumn();
118 0 : for(uInt nt=0;nt<telescopeNames.nelements();nt++)
119 : {
120 0 : if ((telescopeNames(nt) != "VLA") && (telescopeNames(nt) != "EVLA"))
121 : {
122 0 : String mesg="We can handle only (E)VLA antennas for now.\n";
123 0 : mesg += "Erroneous telescope name = " + telescopeNames(nt) + ".";
124 0 : SynthesisError err(mesg);
125 0 : throw(err);
126 : }
127 0 : if (telescopeNames(nt) != telescopeNames(0))
128 : {
129 0 : String mesg="We do not (yet) handle inhomogeneous arrays for A-Projection!\n";
130 0 : mesg += "Not yet a \"priority\"!!";
131 0 : SynthesisError err(mesg);
132 0 : throw(err);
133 : }
134 : }
135 : // MSSpWindowColumns mssp(vb.msColumns().spectralWindow());
136 0 : Freq = vb.msColumns().spectralWindow().refFrequency()(0);
137 0 : Diameter_p=0;
138 0 : Nant_p = vb.msColumns().antenna().nrow();
139 0 : for (Int i=0; i < Nant_p; i++)
140 0 : if (!vb.msColumns().antenna().flagRow()(i))
141 : {
142 0 : Diameter_p = vb.msColumns().antenna().dishDiameter()(i);
143 0 : break;
144 : }
145 0 : if (Diameter_p == 0)
146 : {
147 0 : logIO() << LogOrigin("EVLAConvFunc", "getVisParams")
148 : << "No valid or finite sized antenna found in the antenna table. "
149 : << "Assuming diameter = 25m."
150 : << LogIO::WARN
151 0 : << LogIO::POST;
152 0 : Diameter_p=25.0;
153 : }
154 :
155 0 : Double Lambda=C::c/Freq;
156 0 : HPBW = Lambda/(Diameter_p*sqrt(log(2.0)));
157 0 : sigma = 1.0/(HPBW*HPBW);
158 : // awEij.setSigma(sigma);
159 0 : Int bandID = getVLABandID(Freq,telescopeNames(0));
160 0 : return bandID;
161 : }
162 :
163 0 : Int EVLAConvFunc::makePBPolnCoords(const VisBuffer&vb,
164 : const Vector<Int>& polMap,
165 : const Int& convSize,
166 : const Int& convSampling,
167 : const CoordinateSystem& skyCoord,
168 : const Int& skyNx, const Int& /*skyNy*/,
169 : CoordinateSystem& feedCoord,
170 : Vector<Int>& cfStokes)
171 : {
172 0 : feedCoord = skyCoord;
173 : //
174 : // Make a two dimensional image to calculate auto-correlation of
175 : // the ideal illumination pattern. We want this on a fine grid in
176 : // the UV plane
177 : //
178 0 : Int directionIndex=skyCoord.findCoordinate(Coordinate::DIRECTION);
179 0 : AlwaysAssert(directionIndex>=0, AipsError);
180 0 : DirectionCoordinate dc=skyCoord.directionCoordinate(directionIndex);
181 0 : Vector<Double> sampling;
182 0 : sampling = dc.increment();
183 0 : sampling*=Double(convSampling);
184 0 : sampling*=Double(skyNx)/Double(convSize);
185 0 : dc.setIncrement(sampling);
186 :
187 :
188 0 : Vector<Double> unitVec(2);
189 0 : unitVec=convSize/2;
190 0 : dc.setReferencePixel(unitVec);
191 :
192 : // Set the reference value to that of the image
193 0 : feedCoord.replaceCoordinate(dc, directionIndex);
194 :
195 : //
196 : // Make an image with circular polarization axis.
197 : //
198 0 : Int NPol=0,M,N=0;
199 0 : M=polMap.nelements();
200 0 : for(Int i=0;i<M;i++) if (polMap(i) > -1) NPol++;
201 0 : Vector<Int> poln(NPol);
202 :
203 : Int index;
204 0 : Vector<Int> inStokes;
205 0 : index = feedCoord.findCoordinate(Coordinate::STOKES);
206 0 : inStokes = feedCoord.stokesCoordinate(index).stokes();
207 0 : N = 0;
208 : try
209 : {
210 0 : for(Int i=0;i<M;i++) if (polMap(i) > -1) {poln(N) = vb.corrType()(i);N++;}
211 0 : StokesCoordinate polnCoord(poln);
212 0 : Int StokesIndex = feedCoord.findCoordinate(Coordinate::STOKES);
213 0 : feedCoord.replaceCoordinate(polnCoord,StokesIndex);
214 0 : cfStokes = poln;
215 : }
216 0 : catch(AipsError& x)
217 : {
218 : throw(SynthesisFTMachineError("Likely cause: Discrepancy between the poln. "
219 0 : "axis of the data and the image specifications."));
220 : }
221 :
222 0 : return NPol;
223 : }
224 :
225 0 : Bool EVLAConvFunc::findSupport(Array<Complex>& func, Float& threshold,Int& origin, Int& R)
226 : {
227 : Double NSteps;
228 0 : Int PixInc=1;
229 0 : Vector<Complex> vals;
230 0 : IPosition ndx(4,origin,0,0,0);
231 0 : Bool found=false;
232 0 : IPosition cfShape=func.shape();
233 0 : Int convSize = cfShape(0);
234 0 : for(R=convSize/4;R>1;R--)
235 : {
236 0 : NSteps = 90*R/PixInc; //Check every PixInc pixel along a
237 : //circle of radious R
238 0 : vals.resize((Int)(NSteps+0.5));
239 0 : vals=0;
240 0 : for(Int th=0;th<NSteps;th++)
241 : {
242 0 : ndx(0)=(int)(origin + R*sin(2.0*M_PI*th*PixInc/R));
243 0 : ndx(1)=(int)(origin + R*cos(2.0*M_PI*th*PixInc/R));
244 :
245 0 : if ((ndx(0) < cfShape(0)) && (ndx(1) < cfShape(1)))
246 0 : vals(th)=func(ndx);
247 : }
248 0 : if (max(abs(vals)) > threshold)
249 0 : {found=true;break;}
250 : }
251 0 : return found;
252 : }
253 :
254 0 : void EVLAConvFunc::makeConvFunction(const ImageInterface<Complex>& image,
255 : const VisBuffer& vb,
256 : const Int wConvSize,
257 : // const Vector<Int>& polMap,
258 : Float pa,
259 : Float dpa,
260 : // Vector<Int>& cfStokes,
261 : CFStore& cfs,
262 : CFStore& cfwts,Bool /*fillCF*/)
263 : {
264 0 : LogIO log_l(LogOrigin("EVLAConvFunc", "makeConvFunction"));
265 : Int convSize, convSampling, polInUse;
266 0 : Double wScale=1;
267 0 : Array<Complex> convFunc_l, convWeights_l;
268 :
269 : (void)dpa;
270 :
271 0 : Int nx=image.shape()(0);
272 0 : if (bandID_p == -1) bandID_p=getVisParams(vb);
273 :
274 : log_l << "Making a new convolution function for PA="
275 0 : << pa*(180/C::pi) << "deg"
276 0 : << LogIO::NORMAL << LogIO::POST;
277 :
278 0 : if(wConvSize>0) {
279 0 : log_l << "Using " << wConvSize << " planes for W-projection" << LogIO::POST;
280 : Double maxUVW;
281 0 : maxUVW=0.25/abs(image.coordinates().increment()(0));
282 : log_l << "Estimating maximum possible W = " << maxUVW
283 0 : << " (wavelengths)" << LogIO::POST;
284 :
285 0 : Double invLambdaC=vb.frequency()(0)/C::c;
286 0 : Double invMinL = vb.frequency()((vb.frequency().nelements())-1)/C::c;
287 : log_l << "wavelength range = " << 1.0/invLambdaC << " (m) to "
288 0 : << 1.0/invMinL << " (m)" << LogIO::POST;
289 0 : if (wConvSize > 1)
290 : {
291 0 : wScale=Float((wConvSize-1)*(wConvSize-1))/maxUVW;
292 : log_l << "Scaling in W (at maximum W) = " << 1.0/wScale
293 0 : << " wavelengths per pixel" << LogIO::POST;
294 : }
295 : }
296 : //
297 : // Get the coordinate system
298 : //
299 0 : CoordinateSystem coords(image.coordinates());
300 :
301 : //
302 : // Set up the convolution function.
303 : //
304 0 : if(wConvSize>0)
305 : {
306 0 : if(wConvSize>256)
307 : {
308 0 : convSampling=4;
309 0 : convSize=min(nx,512);
310 : }
311 : else
312 : {
313 0 : convSampling=4;
314 0 : convSize=min(nx,2048);
315 : }
316 : }
317 : else
318 : {
319 0 : convSampling=4;
320 0 : convSize=nx;
321 : }
322 0 : convSampling=OVERSAMPLING;
323 0 : convSize=CONVSIZE;
324 : //
325 : // Make a two dimensional image to calculate auto-correlation of
326 : // the ideal illumination pattern. We want this on a fine grid in
327 : // the UV plane
328 : //
329 :
330 0 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
331 0 : AlwaysAssert(directionIndex>=0, AipsError);
332 0 : DirectionCoordinate dc=coords.directionCoordinate(directionIndex);
333 0 : Vector<Double> sampling;
334 0 : sampling = dc.increment();
335 : // sampling /= Double(2.0);
336 :
337 0 : sampling*=Double(convSampling);
338 0 : sampling*=Double(nx)/Double(convSize);
339 :
340 : // cerr << "Sampling on the sky = " << dc.increment() << " " << nx << "x" << ny << endl;
341 : // cerr << "Sampling on the PB = " << sampling << " " << convSize << "x" << convSize << endl;
342 0 : dc.setIncrement(sampling);
343 :
344 0 : Vector<Double> unitVec(2);
345 0 : unitVec=convSize/2;
346 0 : dc.setReferencePixel(unitVec);
347 :
348 : // Set the reference value to that of the image
349 0 : coords.replaceCoordinate(dc, directionIndex);
350 :
351 : //
352 : // Make an image with circular polarization axis. Return the
353 : // no. of vis. poln. planes that will be used in making the user
354 : // defined Stokes image.
355 : //
356 0 : polInUse=makePBPolnCoords(vb, polMap_p, convSize, convSampling,
357 : image.coordinates(),nx,nx,
358 0 : coords,feedStokes_p);
359 : //------------------------------------------------------------------
360 : //
361 : // Make the sky Stokes PB. This will be used in the gridding
362 : // correction
363 : //
364 : //------------------------------------------------------------------
365 0 : IPosition pbShape(4, convSize, convSize, polInUse, 1);
366 0 : TempImage<Complex> twoDPB(pbShape, coords);
367 :
368 0 : IPosition pbSqShp(pbShape);
369 : // pbSqShp[0] *= 2; pbSqShp[1] *= 2;
370 :
371 0 : unitVec=pbSqShp[0]/2;
372 0 : dc.setReferencePixel(unitVec);
373 : // sampling *= Double(2.0);
374 : // dc.setIncrement(sampling);
375 0 : coords.replaceCoordinate(dc, directionIndex);
376 :
377 0 : TempImage<Complex> twoDPBSq(pbSqShp,coords);
378 : // twoDPB.setMaximumCacheSize(cachesize);
379 0 : twoDPB.set(Complex(1.0,0.0));
380 : // twoDPBSq.setMaximumCacheSize(cachesize);
381 0 : twoDPBSq.set(Complex(1.0,0.0));
382 : //
383 : // Accumulate the various terms that constitute the gridding
384 : // convolution function.
385 : //
386 : //------------------------------------------------------------------
387 :
388 0 : Int inner=convSize/convSampling;
389 : // inner = convSize/2;
390 :
391 : // Vector<Double> cfUVScale(3,0),cfUVOffset(3,0);
392 :
393 : // cfUVScale(0)=Float(twoDPB.shape()(0))*sampling(0);
394 : // cfUVScale(1)=Float(twoDPB.shape()(1))*sampling(1);
395 : // cfUVOffset(0)=Float(twoDPB.shape()(0))/2;
396 : // cfUVOffset(1)=Float(twoDPB.shape()(1))/2;
397 : // // cerr << wScale << " " << cfUVScale << endl;
398 : // // cerr << uvOffset << " " << cfUVOffset << endl;
399 : // // ConvolveGridder<Double, Complex>
400 : // // // ggridder(IPosition(2, inner, inner), cfUVScale, cfUVOffset, "SF");
401 : // // ggridder(IPosition(2, inner, inner), uvScale, uvOffset, "SF");
402 :
403 : // // convFuncCache[PAIndex] = new Array<Complex>(IPosition(4,convSize/2,convSize/2,
404 : // // wConvSize,polInUse));
405 : // // convWeightsCache[PAIndex] = new Array<Complex>(IPosition(4,convSize/2,convSize/2,
406 : // // wConvSize,polInUse));
407 : // // convFuncCache[PAIndex] = new Array<Complex>(IPosition(4,convSize,convSize,
408 : // // wConvSize,polInUse));
409 : // // convWeightsCache[PAIndex] = new Array<Complex>(IPosition(4,convSize,convSize,
410 : // // wConvSize,polInUse));
411 0 : cfs.data = new Array<Complex>(IPosition(4,convSize,convSize,
412 0 : wConvSize,polInUse));
413 0 : cfwts.data = new Array<Complex>(IPosition(4,convSize,convSize,
414 0 : wConvSize,polInUse));
415 0 : convFunc_l.reference(*cfs.data);
416 0 : convWeights_l.reference(*cfwts.data);
417 0 : convFunc_l=0;
418 0 : convWeights_l=0.0;
419 0 : cfs.resize(wConvSize);
420 0 : cfwts.resize(wConvSize);
421 :
422 0 : IPosition start(4, 0, 0, 0, 0);
423 0 : IPosition pbSlice(4, convSize, convSize, 1, 1);
424 :
425 0 : Matrix<Complex> screen(convSize, convSize);
426 : // if (paChangeDetector.changed(vb,0)) paChangeDetector.update(vb,0);
427 0 : VLACalcIlluminationConvFunc vlaPB;
428 0 : WTerm wterm;
429 0 : Long cachesize=(HostInfo::memoryTotal(true)/8)*1024;
430 0 : vlaPB.setMaximumCacheSize(cachesize);
431 :
432 0 : for (Int iw=0;iw<wConvSize;iw++)
433 : {
434 0 : screen = 1.0;
435 :
436 : /*
437 : screen=0.0;
438 : // First the spheroidal function
439 : //
440 : // inner=convSize/2;
441 : // screen = 0.0;
442 : Vector<Complex> correction(inner);
443 : for (Int iy=-inner/2;iy<inner/2;iy++)
444 : {
445 : ggridder.correctX1D(correction, iy+inner/2);
446 : for (Int ix=-inner/2;ix<inner/2;ix++)
447 : screen(ix+convSize/2,iy+convSize/2)=correction(ix+inner/2);
448 : // if (iy==0)
449 : // for (Int ii=0;ii<inner;ii++)
450 : // cout << ii << " " << correction(ii) << endl;
451 : }
452 : */
453 : //
454 : // Now the w term
455 : //
456 0 : wterm.applySky(screen, iw, sampling, wScale, inner);
457 : //
458 : // Fill the complex image with the w-term...
459 : //
460 0 : IPosition PolnPlane(4,0,0,0,0);
461 0 : IPosition ndx(4,0,0,0,0);
462 :
463 0 : for(Int i=0;i<polInUse;i++)
464 : {
465 0 : PolnPlane(2)=i;
466 0 : twoDPB.putSlice(screen, PolnPlane);
467 0 : twoDPBSq.putSlice(screen, PolnPlane);
468 : }
469 : //
470 : // Apply the PB...
471 : //
472 0 : Bool doSquint=true;
473 0 : Double pa=getPA(vb);
474 0 : vlaPB.applyPB(twoDPB, pa, bandID_p, doSquint);
475 0 : doSquint = false;
476 : // vlaPB.applyPBSq(twoDPBSq, vb, bandID_p, doSquint);
477 0 : vlaPB.applyPB(twoDPBSq, pa, bandID_p, doSquint);
478 : /*
479 : // twoDPB.put(abs(twoDPB.get()));
480 : // twoDPBSq.put(abs(twoDPBSq.get()));
481 : */
482 :
483 : // {
484 : // String name("twoDPB.before.im");
485 : // storeImg(name,twoDPB);
486 : // }
487 :
488 0 : Complex cpeak=max(twoDPB.get());
489 0 : twoDPB.put(twoDPB.get()/cpeak);
490 0 : cpeak=max(twoDPBSq.get());
491 0 : twoDPBSq.put(twoDPBSq.get()/cpeak);
492 : // twoDPBSq.set(1.0);
493 : // {
494 : // String name("twoDPB.im");
495 : // storeImg(name,twoDPB);
496 : // }
497 :
498 0 : CoordinateSystem cs=twoDPB.coordinates();
499 0 : Int index= twoDPB.coordinates().findCoordinate(Coordinate::SPECTRAL);
500 0 : SpectralCoordinate SpCS = twoDPB.coordinates().spectralCoordinate(index);
501 :
502 0 : Double cfRefFreq=SpCS.referenceValue()(0);
503 0 : Vector<Double> refValue; refValue.resize(1); refValue(0)=cfRefFreq;
504 0 : SpCS.setReferenceValue(refValue);
505 0 : cs.replaceCoordinate(SpCS,index);
506 : //
507 : // Now FFT and get the result back
508 : //
509 : // {
510 : // String name("twoDPB.im");
511 : // storeImg(name,twoDPB);
512 : // }
513 0 : LatticeFFT::cfft2d(twoDPB);
514 0 : LatticeFFT::cfft2d(twoDPBSq);
515 : // {
516 : // String name("twoDPBFT.im");
517 : // storeImg(name,twoDPB);
518 : // }
519 : //
520 : // Fill the convolution function planes with the result.
521 : //
522 : {
523 0 : IPosition start(4, 0, 0, 0, 0),
524 0 : pbSlice(4, twoDPB.shape()[0]-1, twoDPB.shape()[1]-1, polInUse, 1);
525 0 : IPosition sliceStart(4,0,0,iw,0),
526 0 : sliceLength(4,convFunc_l.shape()[0]-1,convFunc_l.shape()[1]-1,1,polInUse);
527 :
528 0 : convFunc_l(Slicer(sliceStart,sliceLength)).nonDegenerate()
529 0 : =(twoDPB.getSlice(start, pbSlice, true));
530 :
531 0 : IPosition shp(twoDPBSq.shape());
532 :
533 0 : IPosition sqStart(4, 0, 0, 0, 0),
534 0 : pbSqSlice(4, shp[0]-1, shp[1]-1, polInUse, 1);
535 0 : IPosition sqSliceStart(4,0,0,iw,0),
536 0 : sqSliceLength(4,shp[0]-1,shp[1]-1,1,polInUse);
537 :
538 0 : convWeights_l(Slicer(sqSliceStart,sqSliceLength)).nonDegenerate()
539 0 : =(twoDPBSq.getSlice(sqStart, pbSqSlice, true));
540 :
541 : }
542 : }
543 :
544 : {
545 0 : Complex cpeak = max(convFunc_l);
546 0 : convFunc_l/=cpeak;
547 0 : cpeak=max(convWeights_l);
548 0 : convWeights_l/=cpeak;
549 : }
550 : //
551 : // Find the convolution function support size. No assumption
552 : // about the symmetry of the conv. func. can be made (except that
553 : // they are same for all poln. planes).
554 : //
555 :
556 0 : Int maxConvSupport=-1;
557 0 : Int ConvFuncOrigin=convFunc_l.shape()[0]/2; // Conv. Func. is half that size of convSize
558 0 : IPosition ndx(4,ConvFuncOrigin,0,0,0);
559 :
560 0 : Int maxConvWtSupport=0, supportBuffer;
561 0 : for (Int iw=0;iw<wConvSize;iw++)
562 : {
563 0 : Bool found=false;
564 : Float threshold;
565 : Int R;
566 0 : ndx(2) = iw;
567 :
568 0 : ndx(0)=ndx(1)=ConvFuncOrigin;
569 0 : ndx(2) = iw;
570 : // Complex maxVal = max(convFunc);
571 0 : threshold = abs(convFunc_l(ndx))*THRESHOLD;
572 : //
573 : // Find the support size of the conv. function in pixels
574 : //
575 : Int wtR;
576 0 : found = findSupport(convWeights_l,threshold,ConvFuncOrigin,wtR);
577 0 : found = findSupport(convFunc_l,threshold,ConvFuncOrigin,R);
578 :
579 : // R *=2.5;
580 : //
581 : // Set the support size for each W-plane and for all
582 : // Pol-planes. Assuming that support size for all Pol-planes
583 : // is same.
584 : //
585 0 : if(found)
586 : {
587 : // Int maxR=R;//max(ndx(0),ndx(1));
588 0 : cfs.sampling(0)=cfwts.sampling(0)=convSampling;
589 0 : for(Int ipol=0;ipol<polInUse;ipol++)
590 : {
591 0 : cfs.xSupport(iw)=cfs.ySupport(iw)=Int(R/cfs.sampling(0));
592 0 : cfs.xSupport(iw)=cfs.ySupport(iw)=Int(0.5+Float(R)/cfs.sampling(0))+1;
593 :
594 0 : cfwts.xSupport(iw)=cfwts.ySupport(iw)=Int(R*CONVWTSIZEFACTOR/
595 0 : cfwts.sampling(0));
596 0 : cfwts.xSupport(iw)=cfwts.ySupport(iw)=Int(0.5+Float(R)*CONVWTSIZEFACTOR/
597 0 : cfwts.sampling(0))+1;
598 :
599 0 : if (cfs.maxXSupport == -1)
600 0 : if (cfs.xSupport(iw) > maxConvSupport)
601 0 : maxConvSupport = cfs.xSupport(iw);
602 0 : maxConvWtSupport=cfwts.xSupport(iw);//HOW CAN THIS BE RIGHT!!!!
603 : }
604 : }
605 : }
606 :
607 0 : if(cfs.xSupport(0)<1)
608 : log_l << "Convolution function is misbehaved - support seems to be zero"
609 0 : << LogIO::EXCEPTION;
610 :
611 : log_l << "Re-sizing the convolution functions"
612 0 : << LogIO::POST;
613 :
614 : {
615 0 : supportBuffer = OVERSAMPLING;
616 0 : Int bot=(Int)(ConvFuncOrigin-cfs.sampling[0]*maxConvSupport-supportBuffer),//-convSampling/2,
617 0 : top=(Int)(ConvFuncOrigin+cfs.sampling[0]*maxConvSupport+supportBuffer);//+convSampling/2;
618 0 : bot = max(0,bot);
619 0 : top = min(top, convFunc_l.shape()(0)-1);
620 : {
621 0 : Array<Complex> tmp;
622 0 : IPosition blc(4,bot,bot,0,0), trc(4,top,top,wConvSize-1,polInUse-1);
623 : //
624 : // Cut out the conv. func., copy in a temp. array, resize the
625 : // CFStore.data, and copy the cutout version to CFStore.data.
626 : //
627 0 : tmp = convFunc_l(blc,trc);
628 0 : cfs.data->resize(tmp.shape());
629 0 : *cfs.data = tmp;
630 0 : convFunc_l.reference(*cfs.data);
631 : }
632 :
633 0 : supportBuffer = (Int)(OVERSAMPLING*CONVWTSIZEFACTOR);
634 0 : bot=(Int)(ConvFuncOrigin-cfwts.sampling[0]*maxConvWtSupport-supportBuffer);
635 0 : top=(Int)(ConvFuncOrigin+cfwts.sampling[0]*maxConvWtSupport+supportBuffer);
636 0 : bot=max(0,bot);
637 0 : top=min(top,convWeights_l.shape()(0)-1);
638 : {
639 0 : Array<Complex> tmp;
640 0 : IPosition blc(4,bot,bot,0,0), trc(4,top,top,wConvSize-1,polInUse-1);
641 :
642 0 : tmp = convWeights_l(blc,trc);
643 0 : cfwts.data->resize(tmp.shape());
644 0 : *cfwts.data = tmp;
645 0 : convWeights_l.reference(*cfwts.data);
646 : }
647 : }
648 :
649 : //
650 : // Normalize such that plane 0 sums to 1 (when jumping in steps of
651 : // convSampling). This is really not necessary here since we do
652 : // the normalizing by the area more accurately in the gridder
653 : // (fpbwproj.f).
654 : //
655 0 : ndx(2)=ndx(3)=0;
656 :
657 :
658 0 : Complex pbSum=0.0;
659 0 : IPosition peakPix(ndx);
660 :
661 0 : Int Nx = convFunc_l.shape()(0), Ny=convFunc_l.shape()(1);
662 :
663 0 : for(Int nw=0;nw<wConvSize;nw++)
664 0 : for(Int np=0;np<polInUse;np++)
665 : {
666 0 : ndx(2) = nw; ndx(3)=np;
667 : {
668 : //
669 : // Locate the pixel with the peak value. That's the
670 : // origin in pixel co-ordinates.
671 : //
672 0 : Float peak=0;
673 0 : peakPix = 0;
674 0 : for(ndx(1)=0;ndx(1)<convFunc_l.shape()(1);ndx(1)++)
675 0 : for(ndx(0)=0;ndx(0)<convFunc_l.shape()(0);ndx(0)++)
676 0 : if (abs(convFunc_l(ndx)) > peak) {peakPix = ndx;peak=abs(convFunc_l(ndx));}
677 : }
678 :
679 0 : ConvFuncOrigin = peakPix(0);
680 : // ConvFuncOrigin = convFunc.shape()(0)/2+1;
681 : // Int thisConvSupport=convSampling*convSupport(nw,np,lastPASlot);
682 0 : Int thisConvSupport=cfs.xSupport(nw);
683 0 : pbSum=0.0;
684 :
685 0 : for(Int iy=-thisConvSupport;iy<thisConvSupport;iy++)
686 0 : for(Int ix=-thisConvSupport;ix<thisConvSupport;ix++)
687 : {
688 0 : ndx(0)=ix*(Int)cfs.sampling[0]+ConvFuncOrigin;
689 0 : ndx(1)=iy*(Int)cfs.sampling[0]+ConvFuncOrigin;
690 0 : pbSum += real(convFunc_l(ndx));
691 : }
692 0 : if(pbSum>0.0)
693 : {
694 : //
695 : // Normalize each Poln. plane by the area under its convfunc.
696 : //
697 0 : Nx = convFunc_l.shape()(0), Ny = convFunc_l.shape()(1);
698 0 : for (ndx(1)=0;ndx(1)<Ny;ndx(1)++)
699 0 : for (ndx(0)=0;ndx(0)<Nx;ndx(0)++)
700 : {
701 0 : convFunc_l(ndx) /= pbSum;
702 : }
703 :
704 0 : Nx = convWeights_l.shape()(0); Ny = convWeights_l.shape()(1);
705 0 : for (ndx(1)=0; ndx(1)<Ny; ndx(1)++)
706 0 : for (ndx(0)=0; ndx(0)<Nx; ndx(0)++)
707 : {
708 0 : convWeights_l(ndx) /= pbSum*pbSum;
709 : }
710 : }
711 : else
712 0 : throw(SynthesisFTMachineError("Convolution function integral is not positive"));
713 :
714 0 : Vector<Float> maxVal(convWeights_l.shape()(2));
715 0 : Vector<IPosition> posMax(convWeights_l.shape()(2));
716 0 : SynthesisUtils::findLatticeMax(convWeights_l,maxVal,posMax);
717 : }
718 :
719 0 : Int index=coords.findCoordinate(Coordinate::SPECTRAL);
720 0 : SpectralCoordinate spCS = coords.spectralCoordinate(index);
721 0 : Vector<Double> refValue; refValue.resize(1);refValue(0)=spCS.referenceValue()(0);
722 0 : spCS.setReferenceValue(refValue);
723 0 : coords.replaceCoordinate(spCS,index);
724 :
725 0 : cfs.coordSys=coords; cfwts.coordSys=coords;
726 0 : cfs.pa=Quantity(pa,"rad"); cfwts.pa=Quantity(pa,"rad");
727 0 : }
728 : //
729 : //---------------------------------------------------------------------
730 : //---------------------------------------------------------------------
731 : // TEMPS
732 : //---------------------------------------------------------------------
733 : //---------------------------------------------------------------------
734 : //
735 : };
|
