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