Line data Source code
1 : //# ALMAIlluminationConvFunc.cc: Implementation for ALMAIlluminationConvFunc
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2002
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //# Copyright by ESO (in the framework of the ALMA collaboration)
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/ALMACalcIlluminationConvFunc.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 <synthesis/TransformMachines/SynthesisError.h>
38 : #include <casacore/ms/MeasurementSets/MSColumns.h>
39 : #include <casacore/lattices/LEL/LatticeExpr.h>
40 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
41 : #include <casacore/images/Images/ImageRegrid.h>
42 : #include <casacore/images/Images/PagedImage.h>
43 : #include <casacore/casa/Arrays/ArrayMath.h>
44 : #include <casacore/casa/Utilities/Assert.h>
45 : #include <casacore/casa/OS/File.h>
46 : #include <fstream>
47 : #include <sstream>
48 :
49 : using namespace casacore;
50 : namespace casa{
51 :
52 : //
53 : //------------------------------------------------------------------------
54 : //
55 0 : ALMACalcIlluminationConvFunc::ALMACalcIlluminationConvFunc():
56 : IlluminationConvFunc(),
57 0 : otherAntRayPath_p("")
58 : {
59 0 : ap.oversamp = 3;
60 0 : ap.x0=-6.5; ap.y0=-6.5;
61 0 : ap.dx=0.25; ap.dy=0.25;
62 :
63 0 : ap.nx=ap.ny=52;
64 0 : ap.pa=lastPA=18000000;
65 0 : ap.freq=84.; // GHz
66 0 : ap.band = BeamCalc_ALMA_3;
67 0 : IPosition shape(4,ap.nx,ap.ny,4,1);
68 0 : ap.aperture = new TempImage<Complex>();
69 0 : if (maximumCacheSize() > 0) ap.aperture->setMaximumCacheSize(maximumCacheSize());
70 0 : ap.aperture->resize(shape);
71 0 : }
72 :
73 :
74 0 : CoordinateSystem ALMACalcIlluminationConvFunc::makeUVCoords(CoordinateSystem& imageCoordSys,
75 : IPosition& shape,
76 : Double /*refFreq*/)
77 : {
78 0 : CoordinateSystem FTCoords = imageCoordSys;
79 :
80 0 : Int dirIndex=FTCoords.findCoordinate(Coordinate::DIRECTION);
81 0 : DirectionCoordinate dc=imageCoordSys.directionCoordinate(dirIndex);
82 0 : Vector<Bool> axes(2); axes=true;
83 0 : Vector<Int> dirShape(2); dirShape(0)=shape(0);dirShape(1)=shape(1);
84 0 : Coordinate* FTdc=dc.makeFourierCoordinate(axes,dirShape);
85 :
86 0 : FTCoords.replaceCoordinate(*FTdc,dirIndex);
87 0 : delete FTdc;
88 :
89 0 : return FTCoords;
90 0 : }
91 : //----------------------------------------------------------------------
92 : // Write PB to the pbImage
93 : //
94 0 : void ALMACalcIlluminationConvFunc::applyPB(ImageInterface<Float>& pbImage,
95 : const VisBuffer& vb,
96 : Bool doSquint,
97 : Int cfKey)
98 : {
99 0 : CoordinateSystem skyCS(pbImage.coordinates());
100 0 : IPosition skyShape(pbImage.shape());
101 :
102 0 : TempImage<Complex> uvGrid;
103 0 : if (maximumCacheSize() > 0) uvGrid.setMaximumCacheSize(maximumCacheSize());
104 0 : MVFrequency freqQ(vb.msColumns().spectralWindow().refFrequencyQuant()(0));
105 0 : MEpoch obsTime(vb.msColumns().timeQuant()(0));
106 0 : String antType = ALMAAperture::antTypeStrFromType(ALMAAperture::antennaTypesFromCFKey(cfKey)[0]); // take the first antenna
107 0 : Int bandID = BeamCalc::Instance()->getBandID(freqQ.getValue(), "ALMA", ""/*bandname*/, antType, obsTime, otherAntRayPath_p);
108 :
109 0 : regridAperture(skyCS, skyShape, uvGrid, vb, doSquint, bandID);
110 0 : fillPB(*(ap.aperture),pbImage);
111 0 : }
112 0 : void ALMACalcIlluminationConvFunc::applyPB(ImageInterface<Complex>& pbImage,
113 : const VisBuffer& vb,
114 : Bool doSquint,
115 : Int cfKey)
116 : {
117 0 : CoordinateSystem skyCS(pbImage.coordinates());
118 0 : IPosition skyShape(pbImage.shape());
119 :
120 0 : TempImage<Complex> uvGrid;
121 0 : if (maximumCacheSize() > 0) uvGrid.setMaximumCacheSize(maximumCacheSize());
122 :
123 0 : MVFrequency freqQ(vb.msColumns().spectralWindow().refFrequencyQuant()(0));
124 0 : MEpoch obsTime(vb.msColumns().timeQuant()(0));
125 0 : String antType = ALMAAperture::antTypeStrFromType(ALMAAperture::antennaTypesFromCFKey(cfKey)[0]); // take the first antenna
126 0 : String antType2 = ALMAAperture::antTypeStrFromType(ALMAAperture::antennaTypesFromCFKey(cfKey)[1]); // take the first antenna
127 :
128 0 : antType=String("DA");
129 0 : antType2 = antType;
130 :
131 0 : cout << "cfkey, type1, type2 " << cfKey << " " << antType << " " << antType2 << endl;
132 0 : Int bandID = BeamCalc::Instance()->getBandID(freqQ.getValue(), "ALMA", ""/*bandname*/, antType, obsTime, otherAntRayPath_p);
133 :
134 0 : regridAperture(skyCS, skyShape, uvGrid, vb, doSquint, bandID);
135 0 : pbImage.setCoordinateInfo(skyCS);
136 0 : fillPB(*(ap.aperture),pbImage);
137 0 : }
138 :
139 0 : void ALMACalcIlluminationConvFunc::applyPB(ImageInterface<Float>& pbImage,
140 : const String& telescope, const MEpoch& obsTime,
141 : const String& antType0, const String& /*antType1*/,
142 : const MVFrequency& freqQ, Double pa,
143 : Bool doSquint)
144 : {
145 0 : CoordinateSystem skyCS(pbImage.coordinates());
146 0 : IPosition skyShape(pbImage.shape());
147 :
148 0 : TempImage<Complex> uvGrid;
149 0 : if (maximumCacheSize() > 0) uvGrid.setMaximumCacheSize(maximumCacheSize());
150 :
151 0 : Int bandID = BeamCalc::Instance()->getBandID(freqQ.getValue(), telescope, ""/*bandname*/, antType0, obsTime, otherAntRayPath_p);
152 : // antType1 ignored for the moment
153 0 : regridAperture(skyCS, skyShape, uvGrid, telescope, freqQ, pa, doSquint, bandID);
154 0 : fillPB(*(ap.aperture),pbImage);
155 0 : }
156 :
157 0 : void ALMACalcIlluminationConvFunc::applyPB(ImageInterface<Complex>& pbImage,
158 : const String& telescope, const MEpoch& obsTime,
159 : const String& antType0, const String& /*antType1*/,
160 : const MVFrequency& freqQ, Double pa,
161 : Bool doSquint)
162 : {
163 0 : CoordinateSystem skyCS(pbImage.coordinates());
164 0 : IPosition skyShape(pbImage.shape());
165 :
166 0 : TempImage<Complex> uvGrid;
167 0 : if (maximumCacheSize() > 0) uvGrid.setMaximumCacheSize(maximumCacheSize());
168 :
169 0 : Int bandID = BeamCalc::Instance()->getBandID(freqQ.getValue(), telescope, ""/*bandname*/, antType0, obsTime, otherAntRayPath_p);
170 : // antType1 ignored for the moment
171 0 : regridAperture(skyCS, skyShape, uvGrid, telescope, freqQ, pa, doSquint, bandID);
172 0 : pbImage.setCoordinateInfo(skyCS);
173 0 : fillPB(*(ap.aperture),pbImage);
174 0 : }
175 :
176 0 : void ALMACalcIlluminationConvFunc::applyVP(ImageInterface<Complex>& pbImage,
177 : const String& telescope, const MEpoch& obsTime,
178 : const String& antType0, const String& /*antType1*/,
179 : const MVFrequency& freqQ, Double pa,
180 : Bool doSquint)
181 : {
182 0 : CoordinateSystem skyCS(pbImage.coordinates());
183 0 : IPosition skyShape(pbImage.shape());
184 :
185 0 : TempImage<Complex> uvGrid;
186 0 : if (maximumCacheSize() > 0) uvGrid.setMaximumCacheSize(maximumCacheSize());
187 :
188 0 : Int bandID = BeamCalc::Instance()->getBandID(freqQ.getValue(), telescope, ""/*bandname*/, antType0, obsTime, otherAntRayPath_p);
189 : // antType1 ignored for the moment
190 0 : regridAperture(skyCS, skyShape, uvGrid, telescope, freqQ, pa, doSquint, bandID);
191 0 : pbImage.setCoordinateInfo(skyCS);
192 0 : fillVP(*(ap.aperture),pbImage);
193 0 : }
194 :
195 : //
196 : //--------------------------------------------------------------------------
197 : //
198 0 : void ALMACalcIlluminationConvFunc::regridAperture(CoordinateSystem& skyCS,
199 : IPosition& skyShape,
200 : TempImage<Complex>& uvGrid,
201 : const VisBuffer& vb,
202 : Bool doSquint, Int bandID)
203 : {
204 0 : LogIO logIO(LogOrigin("ALMACalcIlluminationConvFunc","regrid"));
205 0 : CoordinateSystem skyCoords(skyCS);
206 :
207 : //UNUSED: Int index;
208 : Float pa;
209 :
210 0 : pa = getPA(vb);
211 :
212 0 : String telescopeName=vb.msColumns().observation().telescopeName().getColumn()[0];
213 :
214 : Float freqHi;
215 0 : Vector<Double> chanFreq = vb.frequency();
216 :
217 0 : freqHi = max(chanFreq);
218 : //freqLo = min(chanFreq);
219 : //Freq = freqHi;
220 :
221 0 : regridAperture(skyCS, skyShape, uvGrid, telescopeName, MVFrequency(freqHi),
222 : pa, doSquint, bandID);
223 0 : }
224 :
225 0 : void ALMACalcIlluminationConvFunc::regridAperture(CoordinateSystem& skyCS,
226 : IPosition& skyShape,
227 : TempImage<Complex>& uvGrid,
228 : const VisBuffer &vb,
229 : const Vector<Float>& paList,
230 : Bool doSquint, Int bandID)
231 : {
232 0 : CoordinateSystem skyCoords(skyCS);
233 :
234 : Float pa, Freq;
235 0 : if (bandID != -1) ap.band = bandID;
236 0 : AlwaysAssert(ap.band>=-1, AipsError);
237 0 : Vector<Double> chanFreq = vb.frequency();
238 :
239 0 : const MSSpWindowColumns& spwCol = vb.msColumns().spectralWindow();
240 0 : ArrayColumn<Double> chanfreq = spwCol.chanFreq();
241 0 : ScalarColumn<Double> reffreq = spwCol.refFrequency();
242 : // Freq = sum(chanFreq)/chanFreq.nelements();
243 :
244 0 : Freq = max(chanfreq.getColumn());
245 0 : IPosition imsize(skyShape);
246 0 : CoordinateSystem uvCoords = makeUVCoords(skyCoords,imsize,Freq);
247 0 : uvGrid.setCoordinateInfo(uvCoords);
248 :
249 0 : Int index = uvCoords.findCoordinate(Coordinate::LINEAR);
250 0 : LinearCoordinate lc=uvCoords.linearCoordinate(index);
251 0 : Vector<Double> incr = lc.increment();
252 0 : Double Lambda = C::c/Freq;
253 0 : ap.freq = Freq/1E9;
254 :
255 0 : ap.nx = skyShape(0); ap.ny = skyShape(1);
256 0 : IPosition apertureShape(ap.aperture->shape());
257 0 : apertureShape(0) = ap.nx; apertureShape(1) = ap.ny;
258 0 : ap.aperture->resize(apertureShape);
259 :
260 0 : TempImage<Complex> tmpAperture;tmpAperture.resize(apertureShape);
261 0 : if (maximumCacheSize() > 0) tmpAperture.setMaximumCacheSize(maximumCacheSize());
262 :
263 0 : ap.dx = abs(incr(0)*Lambda); ap.dy = abs(incr(1)*Lambda);
264 : // cout << ap.dx << " " << incr(0) << endl;
265 : // ap.x0 = -(25.0/(2*ap.dx)+1)*ap.dx; ap.y0 = -(25.0/(2*ap.dy)+1)*ap.dy;
266 : // Following 3 lines go with the ANT tag in BeamCalc.cc
267 : // Double antRadius=BeamCalcGeometryDefaults[ap.band].Rant;
268 : // ap.x0 = -(antRadius/(ap.dx)+1)*ap.dx;
269 : // ap.y0 = -(antRadius/(ap.dy)+1)*ap.dy;
270 : // Following 2 lines go with the PIX tag in BeamCalc.cc
271 0 : ap.x0 = -(ap.nx/2)*ap.dx;
272 0 : ap.y0 = -(ap.ny/2)*ap.dy;
273 :
274 : //
275 : // Accumulate apertures for a list of PA
276 : //
277 0 : Vector<Int> poln(4);
278 0 : poln(0) = Stokes::XX;
279 0 : poln(1) = Stokes::XY;
280 0 : poln(2) = Stokes::YX;
281 0 : poln(3) = Stokes::YY;
282 :
283 0 : for(uInt ipa=0;ipa<paList.nelements();ipa++)
284 : {
285 0 : pa = paList[ipa];
286 :
287 0 : ap.pa=pa;
288 0 : ap.aperture->set(0.0);
289 0 : for(uInt i=0; i<4; i++){
290 0 : BeamCalc::Instance()->calculateApertureLinPol(&ap, poln(i));
291 : }
292 : //
293 : // Set the phase of the aperture function to zero if doSquint==F
294 : // Poln. axis indices
295 : // 0: XX, 1:XY, 2:YX, 3:YY
296 : // This is electic field. 0=> Poln X,
297 : // 1=> Leakage of X->Y
298 : // 2=> Leakage of Y->X
299 : // 3=> Poln Y
300 : //
301 : // The squint is removed in the following code using
302 : // honest-to-god pixel indexing. If this is not the most
303 : // efficient method of doing this in AIPS++ (i.e. instead use
304 : // slices etc.), then this cost will show up in making the
305 : // average PB. Since this goes over each pixel of a full
306 : // stokes (poln. really) complex image, look here (also) for
307 : // optimization (if required).
308 : //
309 0 : if (!doSquint)
310 : {
311 0 : IPosition PolnXIndex(4,0,0,0,0), PolnYIndex(4,0,0,3,0);
312 0 : IPosition tndx(4,0,0,0,0);
313 0 : for(tndx(3)=0;tndx(3)<apertureShape(3);tndx(3)++) // The freq. axis
314 0 : for(tndx(2)=0;tndx(2)<apertureShape(2);tndx(2)++) // The Poln. axis
315 0 : for(tndx(1)=0;tndx(1)<apertureShape(1);tndx(1)++) // The spatial
316 0 : for(tndx(0)=0;tndx(0)<apertureShape(0);tndx(0)++) // axis.
317 : {
318 0 : PolnXIndex(0)=PolnYIndex(0)=tndx(0);
319 0 : PolnXIndex(1)=PolnYIndex(1)=tndx(1);
320 0 : Complex val, Xval, Yval;
321 : Float phase;
322 0 : val = ap.aperture->getAt(tndx);
323 0 : Xval = ap.aperture->getAt(PolnXIndex);
324 0 : Yval = ap.aperture->getAt(PolnYIndex);
325 0 : phase = arg(Xval); Xval=Complex(cos(phase),sin(phase));
326 0 : phase = arg(Yval); Yval=Complex(cos(phase),sin(phase));
327 :
328 0 : if (tndx(2)==0) ap.aperture->putAt(val*conj(Xval),tndx);
329 0 : else if (tndx(2)==1) ap.aperture->putAt(val*conj(Yval),tndx);
330 0 : else if (tndx(2)==2) ap.aperture->putAt(val*conj(Xval),tndx);
331 0 : else if (tndx(2)==3) ap.aperture->putAt(val*conj(Yval),tndx);
332 : }
333 0 : }
334 0 : tmpAperture += *(ap.aperture);
335 : }
336 0 : (ap.aperture)->copyData( tmpAperture);
337 0 : tmpAperture.resize(IPosition(1,1));//Release temp. store.
338 :
339 0 : StokesCoordinate polnCoord(poln);
340 0 : SpectralCoordinate spectralCoord(MFrequency::TOPO,Freq,1.0,0.0);
341 : // uvCoords.addCoordinate(dirCoord);
342 0 : index = uvCoords.findCoordinate(Coordinate::STOKES);
343 0 : uvCoords.replaceCoordinate(polnCoord,index);
344 0 : index = uvCoords.findCoordinate(Coordinate::SPECTRAL);
345 0 : uvCoords.replaceCoordinate(spectralCoord,index);
346 :
347 0 : ap.aperture->setCoordinateInfo(uvCoords);
348 : // if (doSquint==false)
349 : // {
350 : // String name("apperture.im");
351 : // storeImg(name,*(ap.aperture));
352 : // }
353 :
354 0 : ftAperture(*(ap.aperture));
355 0 : }
356 :
357 :
358 0 : void ALMACalcIlluminationConvFunc::regridAperture(CoordinateSystem& skyCS,
359 : IPosition& skyShape,
360 : TempImage<Complex>& /*uvGrid*/,
361 : const String& telescope,
362 : const MVFrequency& freqQ,
363 : Float pa,
364 : Bool doSquint,
365 : Int bandID)
366 : {
367 0 : LogIO logIO(LogOrigin("ALMACalcIlluminationConvFunc","regridAperture"));
368 0 : CoordinateSystem skyCoords(skyCS);
369 :
370 : Int index;
371 0 : if (bandID != -1) ap.band = bandID;
372 0 : AlwaysAssert(ap.band>=-1, AipsError);
373 :
374 0 : String telescopeName=telescope;
375 :
376 : Float Freq;
377 :
378 0 : if (lastPA == pa){
379 : logIO << "Your CPU is being used to do computations for the same PA as for the previous call. Report this!"
380 0 : << LogIO::NORMAL1;
381 : }
382 :
383 0 : lastPA = pa;
384 :
385 0 : Freq = freqQ.getValue();
386 0 : ap.freq = Freq/1E9;
387 :
388 0 : IPosition imsize(skyShape);
389 :
390 0 : CoordinateSystem uvCoords = makeUVCoords(skyCoords,imsize,freqQ.getValue());
391 :
392 0 : index = uvCoords.findCoordinate(Coordinate::LINEAR);
393 0 : LinearCoordinate lc=uvCoords.linearCoordinate(index);
394 0 : Vector<Double> incr = lc.increment();
395 0 : Double Lambda = C::c/Freq;
396 :
397 0 : index = skyCS.findCoordinate(Coordinate::SPECTRAL);
398 0 : SpectralCoordinate SpC = skyCS.spectralCoordinate(index);
399 0 : Vector<Double> refVal = SpC.referenceValue();
400 0 : refVal(0) = freqQ.getValue();
401 0 : SpC.setReferenceValue(refVal);
402 0 : skyCS.replaceCoordinate(SpC,index);
403 :
404 0 : ap.nx = skyShape(0); ap.ny = skyShape(1);
405 0 : IPosition apertureShape(ap.aperture->shape());
406 0 : apertureShape(0) = ap.nx; apertureShape(1) = ap.ny;
407 0 : ap.aperture->resize(apertureShape);
408 0 : ap.dx = abs(incr(0)*Lambda); ap.dy = abs(incr(1)*Lambda);
409 :
410 : //cout << " dx dy " << ap.dx << " " << ap.dy << endl;
411 :
412 0 : ap.x0 = -(ap.nx/2)*ap.dx;
413 0 : ap.y0 = -(ap.ny/2)*ap.dy;
414 :
415 0 : ap.pa=pa;
416 0 : ap.aperture->set(0.0);
417 :
418 0 : Vector<Int> poln(4);
419 0 : poln(0) = Stokes::XX;
420 0 : poln(1) = Stokes::XY;
421 0 : poln(2) = Stokes::YX;
422 0 : poln(3) = Stokes::YY;
423 :
424 0 : for(uInt i=0; i<4; i++){
425 0 : BeamCalc::Instance()->calculateApertureLinPol(&ap, poln(i));
426 : }
427 :
428 : //
429 : // Set the phase of the aperture function to zero if doSquint==F
430 : // Poln. axis indices
431 : // 0: XX, 1:XY, 2:YX, 3:YY
432 : // This is electic field. 0=> Poln X,
433 : // 1=> Leakage of X->Y
434 : // 2=> Leakage of Y->X
435 : // 3=> Poln Y
436 : //
437 : // The squint is removed in the following code using
438 : // honest-to-god pixel indexing. If this is not the most
439 : // efficient method of doing this in AIPS++ (i.e. instead use
440 : // slices etc.), then this cost will show up in making the
441 : // average PB. Since this goes over each pixel of a full
442 : // stokes (poln. really) complex image, look here (also) for
443 : // optimization (if required).
444 : //
445 0 : if (!doSquint){
446 0 : IPosition PolnXIndex(4,0,0,0,0), PolnYIndex(4,0,0,3,0);
447 0 : IPosition tndx(4,0,0,0,0);
448 0 : for(tndx(3)=0;tndx(3)<apertureShape(3);tndx(3)++) // The freq. axis
449 0 : for(tndx(2)=0;tndx(2)<apertureShape(2);tndx(2)++) // The Poln. axis
450 0 : for(tndx(1)=0;tndx(1)<apertureShape(1);tndx(1)++) // The spatial
451 0 : for(tndx(0)=0;tndx(0)<apertureShape(0);tndx(0)++) // axis.
452 : {
453 0 : PolnXIndex(0)=PolnYIndex(0)=tndx(0);
454 0 : PolnXIndex(1)=PolnYIndex(1)=tndx(1);
455 0 : Complex val, Xval, Yval;
456 : Float phase;
457 0 : val = ap.aperture->getAt(tndx);
458 0 : Yval = ap.aperture->getAt(PolnYIndex);
459 0 : phase = arg(Xval); Xval=Complex(cos(phase),sin(phase));
460 0 : phase = arg(Yval); Yval=Complex(cos(phase),sin(phase));
461 :
462 0 : if (tndx(2)==0) ap.aperture->putAt(val*conj(Xval),tndx);
463 0 : else if (tndx(2)==1) ap.aperture->putAt(val*conj(Yval),tndx);
464 0 : else if (tndx(2)==2) ap.aperture->putAt(val*conj(Xval),tndx);
465 0 : else if (tndx(2)==3) ap.aperture->putAt(val*conj(Yval),tndx);
466 : }
467 0 : }
468 :
469 0 : StokesCoordinate polnCoord(poln);
470 0 : SpectralCoordinate spectralCoord(MFrequency::TOPO,Freq,1.0,0.0);
471 :
472 0 : index = uvCoords.findCoordinate(Coordinate::STOKES);
473 0 : Int inStokes = uvCoords.stokesCoordinate(index).stokes()(0);
474 :
475 0 : uvCoords.replaceCoordinate(polnCoord,index);
476 0 : index = uvCoords.findCoordinate(Coordinate::SPECTRAL);
477 0 : uvCoords.replaceCoordinate(spectralCoord,index);
478 :
479 0 : ap.aperture->setCoordinateInfo(uvCoords);
480 : //
481 : // Now FT the re-gridded Fourier plane to get the primary beam.
482 : //
483 0 : cout << "**Writing ALMA Apertures for Pol " << inStokes << " to disk" << endl;
484 0 : String rname("aperture_pol"+String::toString(inStokes)+".im");
485 0 : storeImg(rname, *(ap.aperture) , true);
486 0 : cout << "Done writing apertures to disk" << endl;
487 :
488 0 : ftAperture(*(ap.aperture));
489 :
490 0 : }
491 :
492 0 : void ALMACalcIlluminationConvFunc::fillPB(ImageInterface<Complex>& inImg,
493 : ImageInterface<Complex>& outImg,
494 : Bool Square)
495 : {
496 0 : IPosition imsize(outImg.shape());
497 0 : IPosition ndx(outImg.shape());
498 0 : IPosition inShape(inImg.shape()),inNdx;
499 0 : Vector<Int> inStokes,outStokes;
500 : Int index,s,index1;
501 0 : index = inImg.coordinates().findCoordinate(Coordinate::STOKES);
502 0 : inStokes = inImg.coordinates().stokesCoordinate(index).stokes();
503 0 : index = outImg.coordinates().findCoordinate(Coordinate::STOKES);
504 0 : outStokes = outImg.coordinates().stokesCoordinate(index).stokes();
505 0 : index = outImg.coordinates().findCoordinate(Coordinate::SPECTRAL);
506 0 : index1 = inImg.coordinates().findCoordinate(Coordinate::SPECTRAL);
507 0 : SpectralCoordinate inSpectralCoords = inImg.coordinates().spectralCoordinate(index1);
508 0 : CoordinateSystem outCS = outImg.coordinates();
509 0 : outCS.replaceCoordinate(inSpectralCoords,index);
510 0 : outImg.setCoordinateInfo(outCS);
511 :
512 0 : ndx(3)=0;
513 0 : for(ndx(2)=0;ndx(2)<imsize(2);ndx(2)++){ // The poln axes
514 0 : Bool found = false;
515 0 : for(s=0;s<inShape(2);s++){
516 0 : if (inStokes(s) == outStokes(ndx(2))){
517 0 : found = true;
518 0 : break;
519 : }
520 : }
521 :
522 0 : if(found){
523 0 : for(ndx(0)=0;ndx(0)<imsize(0);ndx(0)++){
524 0 : for(ndx(1)=0;ndx(1)<imsize(1);ndx(1)++){
525 0 : Complex cval;
526 0 : inNdx = ndx; inNdx(2)=s;
527 0 : cval = inImg.getAt(inNdx);
528 0 : if (Square){
529 0 : cval = cval*conj(cval);
530 : }
531 0 : outImg.putAt(cval*outImg.getAt(ndx),ndx);
532 : }
533 : }
534 : }
535 : else{
536 0 : cerr << "Stokes " << outStokes(ndx(2)) << " not found in response image " << endl;
537 : }
538 :
539 : }
540 :
541 0 : }
542 :
543 0 : void ALMACalcIlluminationConvFunc::fillVP(ImageInterface<Complex>& inImg,
544 : ImageInterface<Complex>& outImg,
545 : Bool Square)
546 : {
547 0 : IPosition imsize(outImg.shape());
548 0 : IPosition ndx(outImg.shape());
549 0 : IPosition inShape(inImg.shape()),inNdx;
550 0 : Vector<Int> inStokes,outStokes;
551 : Int index,s,index1;
552 0 : index = inImg.coordinates().findCoordinate(Coordinate::STOKES);
553 0 : inStokes = inImg.coordinates().stokesCoordinate(index).stokes();
554 0 : index = outImg.coordinates().findCoordinate(Coordinate::STOKES);
555 0 : outStokes = outImg.coordinates().stokesCoordinate(index).stokes();
556 0 : index = outImg.coordinates().findCoordinate(Coordinate::SPECTRAL);
557 0 : index1 = inImg.coordinates().findCoordinate(Coordinate::SPECTRAL);
558 0 : SpectralCoordinate inSpectralCoords = inImg.coordinates().spectralCoordinate(index1);
559 0 : CoordinateSystem outCS = outImg.coordinates();
560 0 : outCS.replaceCoordinate(inSpectralCoords,index);
561 0 : outImg.setCoordinateInfo(outCS);
562 :
563 0 : ndx(3)=0;
564 0 : for(ndx(2)=0;ndx(2)<imsize(2);ndx(2)++) // The poln axes
565 : {
566 0 : for(s=0;s<inShape(2);s++) if (inStokes(s) == outStokes(ndx(2))) break;
567 :
568 0 : for(ndx(0)=0;ndx(0)<imsize(0);ndx(0)++)
569 0 : for(ndx(1)=0;ndx(1)<imsize(1);ndx(1)++)
570 : {
571 0 : Complex cval;
572 0 : inNdx = ndx; inNdx(2)=s;
573 0 : cval = inImg.getAt(inNdx);
574 0 : cval = sqrt(sqrt(real(cval*conj(cval))));
575 0 : if (Square) cval = cval*cval;
576 0 : outImg.putAt(cval*outImg.getAt(ndx),ndx);
577 : }
578 : }
579 0 : }
580 :
581 0 : void ALMACalcIlluminationConvFunc::fillPB(ImageInterface<Complex>& inImg,
582 : ImageInterface<Float>& outImg,
583 : Bool Square)
584 : {
585 0 : IPosition imsize(outImg.shape());
586 0 : IPosition ndx(outImg.shape());
587 0 : IPosition inShape(inImg.shape()),inNdx;
588 :
589 0 : Vector<Int> inStokes,outStokes;
590 : Int index,s;
591 0 : index = inImg.coordinates().findCoordinate(Coordinate::STOKES);
592 0 : inStokes = inImg.coordinates().stokesCoordinate(index).stokes();
593 0 : index = outImg.coordinates().findCoordinate(Coordinate::STOKES);
594 0 : outStokes = outImg.coordinates().stokesCoordinate(index).stokes();
595 0 : ndx(3)=0;
596 :
597 0 : for(ndx(2)=0;ndx(2)<imsize(2);ndx(2)++) // The poln axes
598 : {
599 0 : if (outStokes(ndx(2)) == Stokes::I)
600 : {
601 : //
602 : // Fill the outImg wiht (inImg(XX)+inImage(YY))/2
603 : //
604 0 : for(ndx(0)=0;ndx(0)<imsize(0);ndx(0)++)
605 0 : for(ndx(1)=0;ndx(1)<imsize(1);ndx(1)++)
606 : {
607 0 : Complex cval;
608 0 : inNdx = ndx;
609 0 : inNdx(2)=0; cval = inImg.getAt(inNdx);
610 0 : inNdx(2)=3; cval += inImg.getAt(inNdx);
611 0 : cval/2;
612 :
613 0 : if (Square) cval = cval*conj(cval);
614 :
615 0 : outImg.putAt(abs(cval*outImg.getAt(ndx)),ndx);
616 : }
617 : }
618 0 : else if ((outStokes(ndx(2)) == Stokes::XX) ||
619 0 : (outStokes(ndx(2)) == Stokes::XY) ||
620 0 : (outStokes(ndx(2)) == Stokes::YX) ||
621 0 : (outStokes(ndx(2)) == Stokes::YY))
622 : {
623 0 : for(s=0;s<inShape(2);s++) if (inStokes(s) == outStokes(ndx(2))) break;
624 :
625 0 : for(ndx(0)=0;ndx(0)<imsize(0);ndx(0)++)
626 0 : for(ndx(1)=0;ndx(1)<imsize(1);ndx(1)++)
627 : {
628 0 : Complex cval;
629 0 : inNdx = ndx; inNdx(2)=s;
630 0 : cval = inImg.getAt(inNdx);
631 0 : if (Square) cval = cval*conj(cval);
632 0 : outImg.putAt(abs(cval*outImg.getAt(ndx)),ndx);
633 : }
634 : }
635 0 : else throw(AipsError("Unsupported Stokes found in ALMACalcIlluminationConvFunc::fillPB."));
636 : }
637 0 : }
638 :
639 0 : void ALMACalcIlluminationConvFunc::ftAperture(ImageInterface<Complex>& uvgrid)
640 : {
641 : //
642 : // Make SkyJones
643 : //
644 0 : LatticeFFT::cfft2d(uvgrid);
645 : //
646 : // Now make SkyMuller
647 : //
648 0 : skyMuller(uvgrid);
649 0 : }
650 :
651 0 : void ALMACalcIlluminationConvFunc::loadFromImage(String& /*fileName*/)
652 : {
653 0 : throw(AipsError("ALMACalcIlluminationConvFunc::loadFromImage() not yet supported."));
654 : };
655 :
656 0 : void ALMACalcIlluminationConvFunc::skyMuller(ImageInterface<Complex>& skyJones)
657 : {
658 0 : Array<Complex> buf=skyJones.get(),tmp;
659 :
660 0 : IPosition shape(buf.shape());
661 0 : IPosition sliceStart0(4,0,0,0,0),sliceStart1(4,0,0,0,0),
662 0 : sliceLength(4,shape(0),shape(1),1,1);
663 : //
664 : // Giving up on fancy slicing of arrays etc. (the commented code
665 : // below). Just do pixel-by-pixel multiplications for of the
666 : // Jones planes. For now, computing only the diagonal of the
667 : // SkyMuller.
668 : //
669 0 : IPosition t(4,0,0,0,0),n0(4,0,0,0,0),n1(4,0,0,0,0);
670 :
671 : Float peak;
672 0 : peak=0;
673 0 : for(t(2)=0;t(2)<shape(2);t(2)++)
674 0 : for(t(1)=0;t(1)<shape(1);t(1)++)
675 0 : for(t(0)=0;t(0)<shape(0);t(0)++)
676 0 : if (abs(buf(t)) > peak) peak = abs(buf(t));
677 0 : if (peak > 1E-8)
678 0 : for(t(3)=0;t(3)<shape(3);t(3)++) // Freq axis
679 0 : for(t(2)=0;t(2)<shape(2);t(2)++) // Poln axis
680 0 : for(t(1)=0;t(1)<shape(1);t(1)++) // Y axis
681 0 : for(t(0)=0;t(0)<shape(0);t(0)++) // X axis
682 0 : buf(t) = buf(t)/peak;
683 :
684 0 : tmp = buf;
685 :
686 0 : t(0)=t(1)=t(2)=t(3)=0;
687 :
688 0 : t(2)=0; n0(2)=0; n1(2)=0; //XX
689 0 : for( n0(0)=n1(0)=t(0)=0; n0(0)<shape(0); n0(0)++,n1(0)++,t(0)++)
690 0 : for(n0(1)=n1(1)=t(1)=0; n0(1)<shape(1); n0(1)++,n1(1)++,t(1)++)
691 0 : buf(t) = (tmp(n0)*conj(tmp(n1)));
692 :
693 0 : t(2)=1;n0(2)=1;n1(2)=1; //XY
694 0 : for(n0(0)=n1(0)=t(0)=0; n0(0)<shape(0); n0(0)++,n1(0)++,t(0)++)
695 0 : for(n0(1)=n1(1)=t(1)=0; n0(1)<shape(1); n0(1)++,n1(1)++,t(1)++)
696 0 : buf(t) = (tmp(n0)*conj(tmp(n1)));
697 :
698 0 : t(2)=2;n0(2)=2;n1(2)=2; //YX
699 0 : for(n0(0)=n1(0)=t(0)=0; n0(0)<shape(0); n0(0)++,n1(0)++,t(0)++)
700 0 : for(n0(1)=n1(1)=t(1)=0; n0(1)<shape(1); n0(1)++,n1(1)++,t(1)++)
701 0 : buf(t) = (tmp(n0)*conj(tmp(n1)));
702 :
703 0 : t(2)=3;n0(2)=3;n1(2)=3; //YY
704 0 : for(n0(0)=n1(0)=t(0)=0; n0(0)<shape(0); n0(0)++,n1(0)++,t(0)++)
705 0 : for(n0(1)=n1(1)=t(1)=0; n0(1)<shape(1); n0(1)++,n1(1)++,t(1)++)
706 0 : buf(t) = (tmp(n0)*conj(tmp(n1)));
707 : /*
708 : sliceStart0(3)=0; sliceStart1(3)=0;
709 : Slicer s0(sliceStart0,sliceLength),s1(sliceStart1,sliceLength);
710 :
711 : buf(s0) = tmp(s1);
712 : buf(s0) *= tmp(s1);
713 : //
714 : // Muller[1,1]
715 : //
716 : sliceStart0(3)=0; sliceStart1(3)=1;
717 : buf(Slicer(sliceStart0,sliceLength)) = tmp(Slicer(sliceStart0,sliceLength))
718 : *tmp(Slicer(sliceStart1,sliceLength));
719 : //
720 : // Muller[2,2]
721 : //
722 : sliceStart0(3)=1; sliceStart1(3)=0;
723 : buf(Slicer(sliceStart0,sliceLength)) = tmp(Slicer(sliceStart0,sliceLength))
724 : *tmp(Slicer(sliceStart1,sliceLength));
725 : //
726 : // Muller[3,3]
727 : //
728 : sliceStart0(3)=2; sliceStart1(3)=0;
729 : buf(Slicer(sliceStart0,sliceLength)) = tmp(Slicer(sliceStart0,sliceLength))
730 : *tmp(Slicer(sliceStart1,sliceLength));
731 : */
732 0 : skyJones.put(buf);
733 0 : }
734 :
735 : };
|
