Line data Source code
1 : // -*- C++ -*-
2 : //# VLAIlluminationConvFunc.cc: Implementation for VLAIlluminationConvFunc
3 : //# Copyright (C) 1996,1997,1998,1999,2000,2002
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 adressed 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 : //#
28 : //# $Id$
29 :
30 : #define USETABLES 1
31 : #include <synthesis/TransformMachines2/VLAIlluminationConvFunc.h>
32 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
33 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
34 : #include <casacore/coordinates/Coordinates/LinearCoordinate.h>
35 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
36 : #include <synthesis/TransformMachines2/Utils.h>
37 : #include <casacore/lattices/LEL/LatticeExpr.h>
38 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
39 : #include <casacore/images/Images/ImageRegrid.h>
40 : #include <casacore/images/Images/PagedImage.h>
41 : #include <casacore/casa/Arrays/ArrayMath.h>
42 : #include <casacore/casa/OS/File.h>
43 : #include <fstream>
44 :
45 : using namespace casacore;
46 : namespace casa{
47 : namespace refim{
48 : using namespace SynthesisUtils;
49 : //
50 : //------------------------------------------------------------------------
51 : //
52 0 : void VLAIlluminationConvFunc::getIdealConvFunc(Array<Complex>& buf)
53 : {
54 0 : convFunc_p.get(buf);
55 0 : }
56 : //
57 : //------------------------------------------------------------------------
58 : //
59 0 : VLAIlluminationConvFunc::VLAIlluminationConvFunc(String fileName):
60 0 : IlluminationConvFunc(),convFunc_p(),resolution()
61 : {
62 0 : pbRead_p=false;
63 0 : String parts="Re"+fileName;
64 0 : PagedImage<Float> reApp(parts);
65 0 : parts="Im"+fileName;
66 0 : PagedImage<Float> imApp(parts);
67 0 : }
68 : //
69 : //------------------------------------------------------------------------
70 : //
71 0 : void VLAIlluminationConvFunc::load(String& fileName,
72 : Vector<Int>& whichStokes,
73 : Float overSampling,
74 : Bool putCoords)
75 : {
76 : Int Nx,Ny,NPol, NCol;
77 0 : Vector<Int> origin(2);
78 0 : ifstream inp(fileName.c_str());
79 0 : Vector<Float> line;
80 : Float Freq, Dia;
81 :
82 0 : Int nStokes=whichStokes.nelements();
83 0 : Vector<Int> polnMap(nStokes);
84 0 : for(Int i=0;i<nStokes;i++)
85 : {
86 0 : switch(whichStokes(i))
87 : {
88 0 : case Stokes::RR:
89 : {
90 0 : polnMap(i)=3;
91 0 : break;
92 : };
93 0 : case Stokes::RL:
94 : {
95 0 : polnMap(i)=5;
96 0 : break;
97 : }
98 0 : case Stokes::LR:
99 : {
100 0 : polnMap(i)=7;
101 0 : break;
102 : }
103 0 : case Stokes::LL:
104 : {
105 0 : polnMap(i)=9;
106 0 : break;
107 : }
108 : }
109 : }
110 0 : inp >> Nx >> Ny >> NPol >> NCol >> freq_p >> Dia;
111 0 : Freq = freq_p*1E9;
112 0 : line.resize(NCol);
113 0 : IPosition shape(4);
114 0 : shape(0)=(Int)(Nx*overSampling);
115 0 : shape(1)=(Int)(Ny*overSampling);
116 0 : shape(0)=shape(1)=512;
117 0 : shape(2)=nStokes;
118 0 : shape(3)=1;
119 0 : origin(0) = shape(0)/2+1;
120 0 : origin(1) = shape(1)/2+1;
121 : // origin(0) = Int((shape(0)+1)/2);
122 : // origin(1) = Int((shape(1)+1)/2);
123 :
124 0 : convFunc_p.resize(shape);
125 : // reAperture_p.resize(shape);
126 : // imAperture_p.resize(shape);
127 :
128 0 : convFunc_p.set(Complex(0,0));
129 : // reAperture_p.set(0.0);
130 : // imAperture_p.set(0.0);
131 :
132 : Double rx,ry,blockage;
133 : Double re,im;
134 : Double dx0,dy0,dx1,dy1,dx,dy;
135 0 : IPosition ndx(4);
136 0 : dx0 = dy0 = dx1 = dy1 = dx = dy = 0;
137 0 : ndx(3)=0;
138 :
139 0 : for (ndx(0)=0;ndx(0)<Nx;ndx(0)++)
140 0 : for(ndx(1)=0;ndx(1)<Ny;ndx(1)++)
141 : {
142 0 : IPosition ndx1(ndx);
143 : Int i;
144 0 : for(i=0;i<NCol;i++) inp >> line(i);
145 0 : rx = line(0);ry=line(1); blockage=line(2);
146 :
147 0 : if (ndx(0)==0) dx0=rx; else {dx0=dx1;dx1=rx;dx += (dx1-dx0);}
148 0 : if (ndx(1)==0) dy0=ry; else {dy0=dy1;dy1=ry;dy += (dy1-dy0);}
149 :
150 0 : i=3;
151 :
152 0 : if ((ndx(0) < Nx-1) && (ndx(1) < Ny-1))
153 0 : for(ndx(2)=0;ndx(2)<shape(2);ndx(2)++)
154 : {
155 : Float reVal, imVal;
156 0 : i = polnMap(ndx(2));
157 0 : re = line(i);
158 0 : im = line(i+1);
159 :
160 0 : ndx1=ndx;
161 0 : ndx1(0)=origin(0)-(Nx/2)+ndx1(0);
162 0 : ndx1(1)=origin(1)-(Ny/2)+ndx1(1);
163 : // convFunc_p.putAt(Complex(re,im)*blockage,ndx1);
164 : // Float amp=sqrt(re*re+im*im);
165 : // Float phs=-atan2(im,re);
166 0 : if (blockage != 0.0)
167 : {
168 : // reVal = amp*cos(phs);
169 : // imVal = amp*sin(phs);
170 0 : reVal = re; imVal = -im;
171 0 : convFunc_p.putAt(Complex(reVal,imVal),ndx1);
172 : // reAperture_p.putAt(reVal,ndx1);
173 : // imAperture_p.putAt(imVal,ndx1);
174 : }
175 : else
176 : {
177 0 : convFunc_p.putAt(Complex(0,0),ndx1);
178 : }
179 : // if (amp > 1e-4) amp=1.0; else amp=0.0;
180 : // reVal = amp*cos(phs);
181 : // imVal = amp*sin(phs);
182 : // convFunc_p.putAt(Complex(reVal,imVal),ndx1);
183 : // reAperture_p.putAt(reVal,ndx1);
184 : // imAperture_p.putAt(imVal,ndx1);
185 : }
186 : }
187 :
188 0 : dx /= Nx/2; dy /= Ny/2;
189 0 : dy = dx;
190 : //
191 : // Make coordinate systems: 2 Linear, 1 Full Stokes, and 1 spectral
192 : // axis.
193 : //
194 0 : Int nAxis = 2;
195 0 : Vector<String> axisNames(nAxis), axisUnits(nAxis);
196 0 : Vector<Double> refPixel(nAxis), increment(nAxis);
197 0 : Vector<Double> refValue(nAxis);
198 0 : axisNames(0)="UU"; axisNames(1)="VV";
199 0 : axisUnits(0)=axisUnits(1)="lambda";
200 0 : Double Lambda = C::c/(Freq), R=Dia/2;
201 0 : increment(0)=dx*R;
202 0 : increment(1)=dy*R;
203 0 : resolution.resize(2);
204 0 : resolution(0)=-increment(0)/Lambda;resolution(1)=increment(1)/Lambda;
205 :
206 0 : refPixel(0)=origin(0);refPixel(1)=origin(1);
207 0 : refValue(0)=refValue(1)=0.0;
208 :
209 :
210 : //
211 : // Probably a silly way of making a co-ordinate system for UV-plane.
212 : //
213 0 : Matrix<Double> xform(2,2);
214 0 : xform = 0.0; xform.diagonal() = 1.0;
215 : DirectionCoordinate dirCoords(MDirection::J2000,
216 0 : Projection(Projection::SIN),
217 0 : 135*C::pi/180.0, 60*C::pi/180.0,
218 0 : -1*C::pi/180.0, 1*C::pi/180,
219 : xform,
220 0 : 128, 128);
221 0 : Vector<Bool> diraxes(2); diraxes=true;
222 0 : Vector<Int> dirShape(2);
223 0 : dirShape(0)=convFunc_p.shape()(0);
224 0 : dirShape(1)=convFunc_p.shape()(1);
225 0 : Coordinate* FTdc=dirCoords.makeFourierCoordinate(diraxes,dirShape);
226 0 : FTdc->setReferencePixel(refPixel);
227 0 : FTdc->setIncrement(resolution);
228 : /*
229 : LinearCoordinate linearCoords(nAxis);
230 : linearCoords.setWorldAxisNames(axisNames);
231 : linearCoords.setWorldAxisUnits(axisUnits);
232 : linearCoords.setReferencePixel(refPixel);
233 : linearCoords.setReferenceValue(refValue);
234 : linearCoords.setIncrement(resolution);
235 : */
236 :
237 :
238 0 : Vector<Int> stokes(4);
239 0 : stokes(0)=Stokes::RR;
240 0 : stokes(1)=Stokes::RL;
241 0 : stokes(2)=Stokes::LR;
242 0 : stokes(3)=Stokes::LL;
243 0 : StokesCoordinate stokesCoords(whichStokes);
244 :
245 0 : SpectralCoordinate spectralCoords(MFrequency::TOPO,Freq,1.0,0.0);
246 : //
247 : // Make the full coordinate system for the image
248 : //
249 0 : CoordinateSystem cs;
250 : // cs.addCoordinate(linearCoords);
251 0 : cs.addCoordinate(*FTdc);
252 0 : cs.addCoordinate(stokesCoords);
253 0 : cs.addCoordinate(spectralCoords);
254 0 : if (putCoords)
255 0 : convFunc_p.setCoordinateInfo(cs);
256 0 : delete FTdc;
257 : // String fn="vlaAperture.im";
258 : // storeImg(fn,convFunc_p);
259 : // exit(0);
260 0 : };
261 :
262 :
263 0 : CoordinateSystem VLAIlluminationConvFunc::makeUVCoords(CoordinateSystem& imageCoordSys,
264 : IPosition& shape)
265 : {
266 0 : CoordinateSystem FTCoords = imageCoordSys;
267 :
268 0 : Int dirIndex=FTCoords.findCoordinate(Coordinate::DIRECTION);
269 0 : DirectionCoordinate dc=imageCoordSys.directionCoordinate(dirIndex);
270 0 : Vector<Bool> axes(2); axes=true;
271 0 : Vector<Int> dirShape(2); dirShape(0)=shape(0);dirShape(1)=shape(1);
272 0 : Coordinate* FTdc=dc.makeFourierCoordinate(axes,dirShape);
273 0 : FTCoords.replaceCoordinate(*FTdc,dirIndex);
274 0 : delete FTdc;
275 :
276 0 : return FTCoords;
277 : }
278 :
279 0 : void VLAIlluminationConvFunc::applyPB(ImageInterface<Float>& pbImage,
280 : const VisBuffer2& vb)
281 : {
282 0 : CoordinateSystem skyCS(pbImage.coordinates());
283 0 : IPosition skyShape(pbImage.shape());
284 0 : TempImage<Complex> uvGrid(skyShape, skyCS);
285 0 : regridApeture(skyCS, skyShape, uvGrid, vb,false);
286 0 : fillPB(uvGrid,pbImage);
287 0 : }
288 :
289 0 : void VLAIlluminationConvFunc::applyPB(ImageInterface<Complex>& pbImage,
290 : const VisBuffer2& vb)
291 : {
292 0 : CoordinateSystem skyCS(pbImage.coordinates());
293 0 : IPosition skyShape(pbImage.shape());
294 0 : TempImage<Complex> uvGrid(skyShape, skyCS);
295 0 : regridApeture(skyCS, skyShape, uvGrid, vb);
296 0 : fillPB(uvGrid,pbImage);
297 0 : }
298 :
299 0 : void VLAIlluminationConvFunc::regridApeture(CoordinateSystem& skyCS,
300 : IPosition& skyShape,
301 : TempImage<Complex>& uvGrid,
302 : const VisBuffer2& vb,
303 : Bool doSquint)
304 : {
305 0 : CoordinateSystem skyCoords(skyCS);
306 :
307 : Int index;
308 : //Double timeValue = getCurrentTimeStamp(vb);
309 0 : cerr << "###########" << endl << "vb2.feedPa1()(0) used in VLAIlluminationConvFunc::regridAperture()" << endl;
310 0 : Float pa = vb.feedPa1()(0);
311 :
312 0 : IPosition imsize(skyShape);
313 0 : CoordinateSystem uvCoords = SynthesisUtils::makeUVCoords(skyCoords,imsize);
314 0 : CoordinateSystem appCoords(convFunc_p.coordinates());
315 0 : index=uvCoords.findCoordinate(Coordinate::LINEAR);
316 0 : LinearCoordinate uvLC=uvCoords.linearCoordinate(index);
317 0 : index=appCoords.findCoordinate(Coordinate::LINEAR);
318 0 : LinearCoordinate appLC=appCoords.linearCoordinate(index);
319 :
320 0 : Vector<Double> incrUVGrid(2), incrApp(2), ratio(2);
321 0 : incrApp = appLC.increment();
322 0 : incrUVGrid = uvLC.increment();
323 0 : ratio = incrUVGrid/incrApp;
324 :
325 0 : IPosition ndx(imsize),uvndx(2,0,0);
326 :
327 0 : IPosition appShape(convFunc_p.shape());
328 0 : IPosition uvSize(2,imsize(0),imsize(1));
329 0 : IPosition appSize(2,appShape(0),appShape(1));
330 : //
331 : // Extract the linear axes from the UV-coordinate system. Make
332 : // co-ordinate system with only two Linear axes with +PA rotation.
333 : //
334 0 : Matrix<Double> paRot(2,2);
335 0 : paRot(0,0) = cos(pa); paRot(1,0) = +sin(pa);
336 0 : paRot(0,1) = -sin(pa); paRot(1,1) = cos(pa);
337 :
338 0 : Vector<Double> refVal(2);refVal = 0.0;
339 :
340 0 : uvLC.setReferenceValue(refVal);
341 0 : CoordinateSystem onlyUVLinCoords;
342 0 : onlyUVLinCoords.addCoordinate(uvLC);
343 :
344 : //
345 : // Make a co-ordinate system with only 2 Linear axes with the
346 : // resolution of the finer sampled aperture function.
347 0 : index=appCoords.findCoordinate(Coordinate::LINEAR);
348 0 : LinearCoordinate dc=appCoords.linearCoordinate(index);
349 0 : dc.setLinearTransform(paRot);
350 0 : CoordinateSystem onlyAppLinCoords;
351 0 : onlyAppLinCoords.addCoordinate(dc);
352 : //
353 : // Make images with PA rotated co-ordinate system. This is the
354 : // UV-grid consistent with the SkyImage, but rotated by PA and has
355 : // only 2 Linear axes (holds only one poln.).
356 : //
357 : // Put this in a scope so that when the code gets to the FFT, the
358 : // big temp. mem. (regriddedUVGrid) is released.
359 : //
360 : {
361 0 : TempImage<Float> regriddedUVGrid(uvSize, onlyUVLinCoords);
362 : //
363 : // Make a TempImage to hold the real or imag parts of the
364 : // aperture function for this polarization product.
365 : //
366 0 : TempImage<Float> theApp(appSize, onlyAppLinCoords);
367 : //
368 : // Re-grid convFunc_p on uvGrid one polarization axis at a time.
369 : //
370 0 : regriddedUVGrid.set(0.0);
371 0 : uvGrid.set(Complex(0,0));
372 0 : ndx = convFunc_p.shape();
373 0 : ndx(3)=0;
374 :
375 0 : index = uvCoords.findCoordinate(Coordinate::STOKES);
376 0 : StokesCoordinate skyStokesCo=uvCoords.stokesCoordinate(index);
377 0 : Vector<Int> skyStokes = skyStokesCo.stokes();
378 : // cout << "Sky stokes = " << skyStokes << endl;
379 :
380 0 : index = appCoords.findCoordinate(Coordinate::STOKES);
381 0 : StokesCoordinate appStokesCo=appCoords.stokesCoordinate(index);
382 0 : Vector<Int> appStokes = appStokesCo.stokes();
383 : // cout << "Aperture stokes = " << appStokes << endl;
384 :
385 :
386 0 : std::complex<double> aperture;
387 0 : ImageRegrid<Float> ir;
388 0 : IPosition ndx2d(2,0,0);
389 : //char Roter[6] = {'-','|','/','-','\\','|'};
390 : //int RotNdx=0;
391 0 : for(ndx(2)=0;ndx(2)<imsize(2);ndx(2)++) // The Poln. axes
392 : {
393 : //
394 : // Extract the real and imag parts of the aperture function
395 : // for this Poln.
396 : //
397 0 : for(Int F=0;F<2;F++)
398 : {
399 : // cerr << Roter[RotNdx%6] << "\b"; RotNdx++;
400 0 : theApp.set(0);
401 0 : for(ndx2d(0)=0,ndx(0)=0;ndx(0)<appSize(0);ndx(0)++,ndx2d(0)++) // The spatial axes
402 0 : if(F==0) {
403 0 : for(ndx2d(1)=0,ndx(1)=0;ndx(1)<appSize(1);ndx(1)++,ndx2d(1)++)
404 : {
405 0 : aperture = convFunc_p.getAt(ndx);
406 0 : if (!doSquint) aperture=Complex(abs(aperture),0.0);
407 0 : theApp.putAt((std::real)(aperture),ndx2d);
408 : }
409 : }
410 : else {
411 0 : for(ndx2d(1)=0,ndx(1)=0;ndx(1)<appSize(1);ndx(1)++,ndx2d(1)++)
412 : {
413 0 : aperture = convFunc_p.getAt(ndx);
414 0 : if (!doSquint) aperture=Complex(abs(aperture),0.0);
415 0 : theApp.putAt((std::imag)(aperture),ndx2d);
416 : }
417 : }
418 : //
419 : // Re-grid the real and imag parts of the aperture function
420 : // onto the real imag parts of the uvGrid.
421 : //
422 : // cerr << Roter[RotNdx%6] << "\b"; RotNdx++;
423 0 : IPosition whichAxes(2, 0, 1);
424 0 : ir.regrid(regriddedUVGrid, Interpolate2D::LINEAR, whichAxes, theApp);
425 : // ir.regrid(imUVGrid, Interpolate2D::LINEAR, whichAxes, imApp);
426 :
427 : //
428 : // Copy the re-gridded real and imag parts to a complex uvGrid.
429 : //
430 : // cerr << Roter[RotNdx%6] << "\b"; RotNdx++;
431 0 : for(uvndx(0)=0,ndx(0)=0;ndx(0)<imsize(0);ndx(0)++,uvndx(0)++) // The spatial axes
432 0 : for(uvndx(1)=0,ndx(1)=0;ndx(1)<imsize(1);ndx(1)++,uvndx(1)++)
433 : {
434 : // Float re,im;
435 : // re = reUVGrid.getAt(uvndx);
436 : // im = imUVGrid.getAt(uvndx);
437 0 : Complex tmp;
438 0 : tmp = uvGrid.getAt(ndx);
439 0 : if (F==0) tmp = Complex(regriddedUVGrid.getAt(uvndx),imag(tmp));
440 0 : else tmp = Complex(real(tmp),regriddedUVGrid.getAt(uvndx));
441 :
442 0 : uvGrid.putAt(tmp,ndx);
443 : }
444 : // cerr << Roter[RotNdx%6] << "\b"; RotNdx++;
445 : }
446 : }
447 : }
448 : //
449 : // Now FT the re-gridded Fourier plane to get the primary beam.
450 : //
451 0 : ftAperture(uvGrid);
452 0 : }
453 :
454 :
455 0 : void VLAIlluminationConvFunc::fillPB(ImageInterface<Complex>& inImg,
456 : ImageInterface<Complex>& outImg)
457 : {
458 0 : IPosition imsize(outImg.shape());
459 0 : IPosition ndx(outImg.shape());
460 0 : ndx(3)=0;
461 0 : for(ndx(2)=0;ndx(2)<imsize(2);ndx(2)++) // The poln axes
462 0 : for(ndx(0)=0;ndx(0)<imsize(0);ndx(0)++)
463 0 : for(ndx(1)=0;ndx(1)<imsize(1);ndx(1)++)
464 : {
465 0 : Complex cval;
466 0 : cval = inImg.getAt(ndx);
467 0 : outImg.putAt(cval*outImg.getAt(ndx),ndx);
468 : }
469 0 : }
470 :
471 0 : void VLAIlluminationConvFunc::fillPB(ImageInterface<Complex>& inImg,
472 : ImageInterface<Float>& outImg)
473 : {
474 0 : IPosition imsize(outImg.shape());
475 0 : IPosition ndx(outImg.shape());
476 0 : ndx(3)=0;
477 0 : for(ndx(0)=0;ndx(0)<imsize(0);ndx(0)++)
478 0 : for(ndx(1)=0;ndx(1)<imsize(1);ndx(1)++)
479 : {
480 : //
481 : // Average along the polarization axes and fillin the
482 : // amp. of the average in the output image.
483 : //
484 0 : Complex cval=0.0;
485 : // for(ndx(2)=0;ndx(2)<imsize(2);ndx(2)++)
486 : // cval += inImg.getAt(ndx);
487 : // cval /= imsize(2);
488 0 : ndx(2)=0;cval = inImg.getAt(ndx);
489 0 : for(ndx(2)=0;ndx(2)<imsize(2);ndx(2)++) // The poln axes
490 0 : outImg.putAt(abs(cval*outImg.getAt(ndx)),ndx);
491 : }
492 0 : }
493 :
494 0 : void VLAIlluminationConvFunc::ftAperture(String& fileName,
495 : Vector<Int>& whichStokes,
496 : Float& overSampling,
497 : const CoordinateSystem& coordSys)
498 : {
499 0 : load(fileName,whichStokes,overSampling,false);
500 0 : CoordinateSystem pbCoords(coordSys);
501 0 : Int dirIndex=pbCoords.findCoordinate(Coordinate::DIRECTION);
502 0 : DirectionCoordinate dc=coordSys.directionCoordinate(dirIndex);
503 0 : Double Lambda=C::c/(freq_p*1E9);
504 0 : IPosition shape(convFunc_p.shape());
505 :
506 0 : resolution(0) = (Lambda/(resolution(0)*shape(0)));
507 0 : resolution(1) = (Lambda/(resolution(1)*shape(1)));
508 :
509 0 : dc.setIncrement(resolution);
510 :
511 0 : Vector<Double> refPix(2),refValue(2);
512 0 : refPix(0)=shape(0)/2+1;
513 0 : refPix(1)=shape(1)/2+1;
514 0 : refValue(0)=refValue(1)=0;
515 0 : dc.setReferencePixel(refPix);
516 :
517 0 : pbCoords.replaceCoordinate(dc,dirIndex);
518 :
519 0 : convFunc_p.setCoordinateInfo(pbCoords);
520 0 : ftAperture();
521 0 : }
522 :
523 0 : void VLAIlluminationConvFunc::ftAperture(TempImage<Complex>& uvgrid)
524 : {
525 : // String fn("reUVGrid.im");
526 : // storeImg(fn,uvgrid);
527 :
528 0 : LatticeFFT::cfft2d(uvgrid);
529 :
530 0 : Array<Complex> buf=uvgrid.get();
531 0 : buf *= conj(buf);
532 :
533 : // Float peak = abs(max(buf));
534 : // buf /= Complex(peak,0.0);
535 : // cout << "Peak = " << peak << endl;
536 :
537 0 : uvgrid.put(buf);
538 :
539 : // String fName = "vlapb.im";
540 : // storeImg(fName,uvgrid);
541 :
542 0 : }
543 :
544 0 : void VLAIlluminationConvFunc::store(String& fileName){storeImg(fileName,convFunc_p);}
545 :
546 0 : void VLAIlluminationConvFunc::storeImg(String& fileName,ImageInterface<Complex>& theImg)
547 : {
548 0 : ostringstream reName,imName;
549 0 : reName << "re" << fileName;
550 0 : imName << "im" << fileName;
551 : {
552 0 : PagedImage<Float> tmp(theImg.shape(), theImg.coordinates(), reName);
553 0 : LatticeExpr<Float> le(abs(theImg));
554 0 : tmp.copyData(le);
555 : }
556 : {
557 0 : PagedImage<Float> tmp(theImg.shape(), theImg.coordinates(), imName);
558 0 : LatticeExpr<Float> le(arg(theImg));
559 0 : tmp.copyData(le);
560 : }
561 0 : }
562 :
563 0 : void VLAIlluminationConvFunc::storeImg(String& fileName,ImageInterface<Float>& theImg)
564 : {
565 0 : PagedImage<Float> tmp(theImg.shape(), theImg.coordinates(), fileName);
566 0 : LatticeExpr<Float> le(theImg);
567 0 : tmp.copyData(le);
568 0 : }
569 :
570 0 : void VLAIlluminationConvFunc::storePB(String& fileName)
571 : {
572 : {
573 0 : ostringstream Name;
574 0 : Name << "re" << fileName;
575 0 : IPosition newShape(convFunc_p.shape());
576 0 : newShape(0)=newShape(1)=200;
577 0 : PagedImage<Float> tmp(newShape, convFunc_p.coordinates(), Name);
578 : // PagedImage<Float> tmp(convFunc_p.shape(), FTCoords, Name);
579 0 : LatticeExpr<Float> le(real(convFunc_p));
580 0 : tmp.copyData(le);
581 : }
582 : {
583 0 : ostringstream Name;
584 0 : Name << "im" << fileName;
585 :
586 0 : IPosition newShape(convFunc_p.shape());
587 0 : newShape(0)=newShape(1)=200;
588 0 : PagedImage<Float> tmp(newShape, convFunc_p.coordinates(), Name);
589 : // PagedImage<Float> tmp(convFunc_p.shape(), FTCoords, Name);
590 0 : LatticeExpr<Float> le(imag(convFunc_p));
591 0 : tmp.copyData(le);
592 : }
593 0 : }
594 0 : void VLAIlluminationConvFunc::loadFromImage(String& /*fileName*/) {};
595 : };
596 : };
|
