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/TransformMachines2/VLACalcIlluminationConvFunc.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 <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 : #include <casacore/casa/OS/Timer.h>
49 : #include <msvis/MSVis/VisibilityIterator2.h>
50 : using namespace casacore;
51 : namespace casa{
52 : namespace refim{
53 : //
54 : //------------------------------------------------------------------------
55 : //
56 148 : VLACalcIlluminationConvFunc::VLACalcIlluminationConvFunc():
57 148 : IlluminationConvFunc(),convFunc_p(),resolution(),pbRead_p(false),freq_p(0),lastPA(0),ap()
58 : {
59 :
60 296 : LogIO logIO(LogOrigin("VLACalcIlluminationConvFunc","ctor"));
61 148 : ap.oversamp = 3;
62 148 : ap.x0=-13.0; ap.y0=-13.0;
63 148 : ap.dx=0.5; ap.dy=0.5;
64 :
65 148 : ap.nx=ap.ny=104;
66 148 : ap.pa=lastPA=18000000;
67 148 : ap.freq=1.4;
68 148 : ap.band = BeamCalc_VLA_L;
69 : // IPosition shape(4,ap.nx,ap.ny,4,1);
70 148 : IPosition shape(4,ap.nx,ap.ny,1,1);//Set poln. axis to be
71 : //degenerate (len=1). This is
72 : //set to 2 if cross-hand
73 : //functions are requested.
74 148 : ap.aperture = new TempImage<Complex> ();
75 :
76 148 : if (maximumCacheSize() > 0) ap.aperture->setMaximumCacheSize(maximumCacheSize());
77 148 : ap.aperture->resize(shape);
78 :
79 148 : }
80 :
81 :
82 148 : CoordinateSystem VLACalcIlluminationConvFunc::makeUVCoords(const CoordinateSystem& imageCoordSys,
83 : const IPosition& shape,
84 : Double /*refFreq*/)
85 : {
86 148 : CoordinateSystem FTCoords = imageCoordSys;
87 148 : Int dirIndex=FTCoords.findCoordinate(Coordinate::LINEAR);
88 :
89 : // If LINEAR axis is found, assume that the coordsys is alread in FT domain
90 148 : if (dirIndex >= 0) return FTCoords;
91 :
92 64 : dirIndex=FTCoords.findCoordinate(Coordinate::DIRECTION);
93 : //cerr << "DIRIndex " << dirIndex << " shape " << shape << endl;
94 64 : DirectionCoordinate dc=imageCoordSys.directionCoordinate(dirIndex);
95 64 : Vector<Bool> axes(2); axes=true;
96 64 : Vector<Int> dirShape(2); dirShape(0)=shape(0);dirShape(1)=shape(1);
97 64 : Coordinate* FTdc=dc.makeFourierCoordinate(axes,dirShape);
98 : // if (refFreq > 0)
99 : // {
100 : // Int index1 = FTCoords.findCoordinate(Coordinate::SPECTRAL);
101 : // SpectralCoordinate SpC = FTCoords.spectralCoordinate(index1);
102 : // Vector<Double> refVal = SpC.referenceValue();
103 : // refVal = refFreq;
104 : // SpC.setReferenceValue(refVal);
105 : // }
106 :
107 64 : FTCoords.replaceCoordinate(*FTdc,dirIndex);
108 64 : delete FTdc;
109 :
110 64 : return FTCoords;
111 64 : }
112 :
113 : // CoordinateSystem VLACalcIlluminationConvFunc::makeJonesCoords(CoordinateSystem& imageCoordsys)
114 : //
115 : //
116 : //----------------------------------------------------------------------
117 : // Write PB to the pbImage
118 : //
119 0 : void VLACalcIlluminationConvFunc::applyPB(ImageInterface<Float>& /*pbImage*/,
120 : //const VisBuffer2& vb,
121 : Double& /*pa*/,
122 : const Vector<Float>& /*paList*/,
123 : Int /*bandID*/, Bool /*doSquint*/)
124 : {
125 0 : throw(AipsError("applyPB(paList) called!"));
126 : // CoordinateSystem skyCS(pbImage.coordinates());
127 : // IPosition skyShape(pbImage.shape());
128 : // TempImage<Complex> uvGrid;
129 : // if (maximumCacheSize() > 0) uvGrid.setMaximumCacheSize(maximumCacheSize());
130 : // // regridAperture(skyCS, skyShape, uvGrid, vb, paList, false, bandID);
131 : // regridAperture(skyCS, skyShape, uvGrid, pa, paList, doSquint, bandID);
132 :
133 : // fillPB(*(ap.aperture),pbImage);
134 : }
135 0 : void VLACalcIlluminationConvFunc::applyPB(ImageInterface<Float>& pbImage,
136 : //const VisBuffer2& vb,
137 : Double& pa,
138 : Int bandID,
139 : Bool doSquint, Double freqVal)
140 : {
141 0 : CoordinateSystem skyCS(pbImage.coordinates());
142 0 : IPosition skyShape(pbImage.shape());
143 :
144 0 : TempImage<Complex> uvGrid;
145 0 : if (maximumCacheSize() > 0) uvGrid.setMaximumCacheSize(maximumCacheSize());
146 : // regridAperture(skyCS, skyShape, uvGrid, vb,false, bandID);
147 0 : regridAperture(skyCS, skyShape, uvGrid, pa, doSquint, bandID, 0, freqVal);
148 0 : fillPB(*(ap.aperture),pbImage);
149 0 : }
150 148 : void VLACalcIlluminationConvFunc::applyPB(ImageInterface<Complex>& pbImage,
151 : //const VisBuffer2& vb,
152 : Double& pa,
153 : Bool doSquint, Int bandID,Int muellerTerm, Double freqVal)
154 : {
155 148 : CoordinateSystem skyCS(pbImage.coordinates());
156 148 : IPosition skyShape(pbImage.shape());
157 : // cout<<"M = "<<muellerTerm<<"\n";
158 : // cout<<"bandID = "<<bandID<<"\n";
159 : // cout<<"doSquint = "<<doSquint<<"\n";
160 148 : TempImage<Complex> uvGrid;
161 :
162 148 : if (maximumCacheSize() > 0) uvGrid.setMaximumCacheSize(maximumCacheSize());
163 : // regridAperture(skyCS, skyShape, uvGrid, vb, true, bandID);
164 148 : Int index= skyCS.findCoordinate(Coordinate::SPECTRAL);
165 148 : SpectralCoordinate spCS = skyCS.spectralCoordinate(index);
166 : // cout<<"Ref Freq for sky jones is :"<<spCS.referenceValue()(0);
167 : // index=skyCS.findCoordinate(Coordinate::DIRECTION);
168 : // AlwaysAssert(index>=0, AipsError);
169 :
170 148 : regridAperture(skyCS, skyShape, uvGrid, pa, doSquint, bandID, muellerTerm,freqVal);
171 :
172 148 : pbImage.setCoordinateInfo(skyCS);
173 : // {
174 : // string name("aperture.im");
175 : // storeImg(name,*(ap.aperture));
176 : // }
177 148 : fillPB(*(ap.aperture),pbImage);
178 : //{
179 : // string name("apb.im");
180 : // storeImg(name,pbImage);
181 : //}
182 148 : }
183 : //--------------------------------------------------------------------------
184 : // Write PB^2 to the pbImage
185 : //
186 0 : void VLACalcIlluminationConvFunc::applyPBSq(ImageInterface<Float>& /*pbImage*/,
187 : //const VisBuffer2& vb,
188 : Double& /*pa*/,
189 : const Vector<Float>& /*paList*/,
190 : Int /*bandID*/,
191 : Bool /*doSquint*/)
192 : {
193 0 : throw(AipsError("applyPBSq(paList) called!"));
194 : // CoordinateSystem skyCS(pbImage.coordinates());
195 : // IPosition skyShape(pbImage.shape());
196 : // TempImage<Complex> uvGrid;
197 : // if (maximumCacheSize() > 0) uvGrid.setMaximumCacheSize(maximumCacheSize());
198 : // // regridAperture(skyCS, skyShape, uvGrid, vb, paList, false, bandID);
199 : // regridAperture(skyCS, skyShape, uvGrid, pa, paList, doSquint, bandID);
200 :
201 : // fillPB(*(ap.aperture),pbImage, true);
202 : }
203 0 : void VLACalcIlluminationConvFunc::applyPBSq(ImageInterface<Float>& pbImage,
204 : //const VisBuffer2& vb,
205 : Double& pa,
206 : Int bandID,
207 : Bool doSquint)
208 : {
209 0 : CoordinateSystem skyCS(pbImage.coordinates());
210 0 : IPosition skyShape(pbImage.shape());
211 :
212 0 : TempImage<Complex> uvGrid;
213 0 : if (maximumCacheSize() > 0) uvGrid.setMaximumCacheSize(maximumCacheSize());
214 : // regridAperture(skyCS, skyShape, uvGrid, vb,false, bandID);
215 0 : regridAperture(skyCS, skyShape, uvGrid, pa, doSquint, bandID);
216 0 : fillPB(*(ap.aperture),pbImage,true);
217 0 : }
218 0 : void VLACalcIlluminationConvFunc::applyPBSq(ImageInterface<Complex>& pbImage,
219 : //const VisBuffer2& vb,
220 : Double& pa,
221 : Int bandID,
222 : Bool doSquint)
223 : {
224 0 : CoordinateSystem skyCS(pbImage.coordinates());
225 0 : IPosition skyShape(pbImage.shape());
226 :
227 0 : TempImage<Complex> uvGrid;
228 0 : if (maximumCacheSize() > 0) uvGrid.setMaximumCacheSize(maximumCacheSize());
229 : // regridAperture(skyCS, skyShape, uvGrid, vb, true, bandID);
230 0 : regridAperture(skyCS, skyShape, uvGrid, pa, doSquint, bandID);
231 0 : fillPB(*(ap.aperture),pbImage, true);
232 0 : }
233 : //
234 : //--------------------------------------------------------------------------
235 : //
236 148 : void VLACalcIlluminationConvFunc::setApertureParams(ApertureCalcParams& ap,
237 : const Float& Freq, const Float& pa,
238 : const Int& bandID,
239 : const IPosition& skyShape,
240 : const Vector<Double>& uvIncr)
241 : {
242 148 : Double Lambda = C::c/Freq;
243 :
244 148 : ap.pa=pa;
245 148 : ap.band = bandID;
246 148 : ap.freq = Freq/1E9;
247 148 : ap.nx = skyShape(0); ap.ny = skyShape(1);
248 148 : ap.dx = abs(uvIncr(0)*Lambda); ap.dy = abs(uvIncr(1)*Lambda);
249 148 : ap.x0 = -(ap.nx/2)*ap.dx; ap.y0 = -(ap.ny/2)*ap.dy;
250 : //cerr << "pa= " << ap.pa << " band " << ap.band << " freq " << ap.freq << " nx ny " << ap.nx << " " << ap.ny << " dx dy " << ap.dx << " " << ap.dy << endl;
251 : //
252 : // If cross-hand pols. are requested, we need to compute both
253 : // the parallel-hand aperture illuminations.
254 : //
255 : //if ((inStokes == Stokes::RL) || (inStokes == Stokes::LR))
256 : {
257 148 : IPosition apShape(ap.aperture->shape());
258 : //cerr << "APshape " << apShape << endl;
259 148 : apShape(3)=4;
260 148 : ap.aperture->resize(apShape);
261 148 : }
262 148 : }
263 : //
264 : //--------------------------------------------------------------------------
265 : //
266 148 : void VLACalcIlluminationConvFunc::regridApertureEngine(ApertureCalcParams& ap,
267 : const Int& /*inStokes*/)
268 : {
269 148 : IPosition apertureShape(ap.aperture->shape());
270 148 : apertureShape(0) = ap.nx; apertureShape(1) = ap.ny;
271 : //cerr << "new aperture shape " << apertureShape << " old " << (ap.aperture->shape())<< endl;
272 148 : ap.aperture->resize(apertureShape);
273 148 : ap.aperture->set(0.0);
274 : //BeamCalc::Instance()->calculateAperture(&ap,inStokes);
275 : //cerr << ap.aperture->shape() << " " << inStokes << endl;
276 :
277 : // If full-pol. imaging, compute all 4 pols., else only the one given by inStokes.
278 148 : BeamCalc::Instance()->calculateAperture(&ap);// The call in the absence of instokes allows the computation of all
279 : //cerr << "Min max of Aperture " << min(ap.aperture->get()) << " " << max(ap.aperture->get()) << endl;
280 : //BeamCalc::Instance()->calculateAperture(&ap,inStokes);// The call in the absence of instokes allows the computation of all
281 : // the four jones parameters at one time.
282 148 : }
283 : //
284 : //--------------------------------------------------------------------------
285 : //
286 148 : void VLACalcIlluminationConvFunc::regridAperture(CoordinateSystem& skyCS,
287 : IPosition& skyShape,
288 : TempImage<Complex>& /*uvGrid*/,
289 : //const VisBuffer2& vb,
290 : Double& pa,
291 : Bool doSquint, Int bandID,Int muellerTerm ,Double freqVal)
292 : {
293 296 : LogIO logIO(LogOrigin("VLACalcIlluminationConvFunc","regrid"));
294 148 : CoordinateSystem skyCoords(skyCS);
295 :
296 : Int index;
297 : //UNUSED: Double timeValue = getCurrentTimeStamp(vb);
298 148 : AlwaysAssert(bandID>=-1, AipsError);
299 148 : if (bandID != -1) ap.band = bandID;
300 : //Float pa = getPA(vb);
301 : Float Freq, freqHi;
302 :
303 148 : if (lastPA == pa)
304 : {
305 : // LogIO logIO;
306 : logIO << "Your CPU is being used to do computations for the same PA as for the previous call. Report this!"
307 0 : << LogIO::NORMAL1;
308 : }
309 :
310 :
311 148 : if (freqVal > 0)
312 : {
313 148 : Freq=freqHi=freqVal;
314 148 : ap.freq=freqHi/1E09;
315 : }
316 : else
317 : {
318 : //throw(AipsError("Freq. < 0 in VLACICF::regrid"));
319 :
320 : // Vector<Double> chanFreq = vb.frequency();
321 0 : index = skyCS.findCoordinate(Coordinate::SPECTRAL);
322 0 : SpectralCoordinate SpC = skyCS.spectralCoordinate(index);
323 0 : Vector<Double> refVal = SpC.referenceValue();
324 :
325 : // freqHi = max(chanFreq);
326 0 : freqHi = refVal[0];
327 : // freqLo = min(chanFreq);
328 0 : Freq = freqHi ;
329 0 : ap.freq = Freq/1E9;
330 0 : }
331 :
332 148 : IPosition imsize(skyShape);
333 148 : CoordinateSystem uvCoords = makeUVCoords(skyCoords,imsize,freqHi);
334 :
335 148 : index = uvCoords.findCoordinate(Coordinate::LINEAR);
336 148 : LinearCoordinate lc=uvCoords.linearCoordinate(index);
337 148 : Vector<Double> uvIncr = lc.increment();
338 : //Double Lambda = C::c/freqHi;
339 :
340 148 : index = uvCoords.findCoordinate(Coordinate::STOKES);
341 : //cerr << "STOKES index " << index << endl;
342 148 : Int inStokes = uvCoords.stokesCoordinate(index).stokes()(0);
343 :
344 : //Vector<Int> intSkyShape=skyShape.asVector();
345 148 : setApertureParams(ap, Freq, pa, bandID,
346 : skyShape, uvIncr);
347 148 : (ap.aperture)->setCoordinateInfo(uvCoords);
348 148 : regridApertureEngine(ap, inStokes);
349 148 : IPosition apertureShape(ap.aperture->shape());
350 :
351 : // ap.freq = Freq/1E9;
352 : // ap.nx = skyShape(0); ap.ny = skyShape(1);
353 : // ap.dx = abs(uvIncr(0)*Lambda); ap.dy = abs(uvIncr(1)*Lambda);
354 : // ap.x0 = -(ap.nx/2)*ap.dx;
355 : // ap.y0 = -(ap.ny/2)*ap.dy;
356 : // ap.pa=pa;
357 : // if ((inStokes == Stokes::RL) || (inStokes == Stokes::LR))
358 : // ap.aperture->shape()(3)=2;
359 :
360 : // apertureShape(0) = ap.nx; apertureShape(1) = ap.ny;
361 : // ap.aperture->resize(apertureShape);
362 : // ap.aperture->set(0.0);
363 :
364 : // BeamCalc::Instance()->calculateAperture(&ap,inStokes);
365 :
366 148 : if (!doSquint)
367 : {
368 0 : IPosition PolnRIndex(4,0,0,0,0), PolnLIndex(4,0,0,3,0);
369 0 : IPosition tndx(4,0,0,0,0);
370 0 : for(tndx(3)=0;tndx(3)<apertureShape(3);tndx(3)++) // The freq. axis
371 0 : for(tndx(2)=0;tndx(2)<apertureShape(2);tndx(2)++) // The Poln. axis
372 0 : for(tndx(1)=0;tndx(1)<apertureShape(1);tndx(1)++) // The spatial
373 0 : for(tndx(0)=0;tndx(0)<apertureShape(0);tndx(0)++) // axis.
374 : {
375 0 : PolnRIndex(0)=PolnLIndex(0)=tndx(0);
376 0 : PolnRIndex(1)=PolnLIndex(1)=tndx(1);
377 0 : Complex val, Rval;
378 : Float phase;
379 0 : val = ap.aperture->getAt(tndx);
380 0 : Rval = ap.aperture->getAt(PolnRIndex);
381 : //Lval = ap.aperture->getAt(PolnLIndex);
382 0 : phase = arg(Rval); Rval=Complex(cos(phase),sin(phase)); // Isn't this garbage ?
383 : //So when tndx(3)=0 and tndx(2)=0 polnRindex and tndx are the same ...thus then
384 : // first plane is Ae^(-i\Rval)*e^(iRval) which is just A of each pixel
385 : // when tndx(3) and (2) moves to 1 ...now Rval is 0
386 : //After the first plane phase is zeroed i.e done modular...then the next planes are being multiplied by conj (0)
387 : //phase = arg(Lval); Lval=Complex(cos(phase),sin(phase));
388 :
389 : // if (tndx(2)==0) ap.aperture->putAt(val*conj(Rval),tndx);
390 : // else if (tndx(2)==1) ap.aperture->putAt(val*conj(Lval),tndx);
391 : // else if (tndx(2)==2) ap.aperture->putAt(val*conj(Rval),tndx);
392 : // else if (tndx(2)==3) ap.aperture->putAt(val*conj(Lval),tndx);
393 0 : ap.aperture->putAt(val*conj(Rval),tndx);
394 : //cerr << "no squint " << val << " " << Rval << endl;
395 : }
396 0 : }
397 : // cout<<"Completed the regrid Aperture step";
398 : // Vector<Int> poln(4);
399 : // poln(0) = Stokes::RR;
400 : // poln(1) = Stokes::RL;
401 : // poln(2) = Stokes::LR;
402 : // poln(3) = Stokes::LL;
403 148 : Vector<Int> poln(1); poln(0)=inStokes;
404 148 : StokesCoordinate polnCoord(poln);
405 148 : SpectralCoordinate spectralCoord(MFrequency::TOPO,Freq,1.0,0.0);
406 : // uvCoords.addCoordinate(dirCoord);
407 148 : index = uvCoords.findCoordinate(Coordinate::STOKES);
408 148 : uvCoords.replaceCoordinate(polnCoord,index);
409 148 : index = uvCoords.findCoordinate(Coordinate::SPECTRAL);
410 148 : uvCoords.replaceCoordinate(spectralCoord,index);
411 : //logIO << "The Stokes coordinate is", poln(0)<< LogIO::POST;
412 148 : ap.aperture->setCoordinateInfo(uvCoords);
413 148 : if (doSquint==true)
414 : {
415 : // String name("aperture.im");
416 : // storeImg(name,*(ap.aperture));
417 : }
418 :
419 : //
420 : // Now FT the re-gridded Fourier plane to get the primary beam.
421 : //
422 148 : ftAperture(*(ap.aperture),muellerTerm);
423 148 : if (doSquint==true)
424 : {
425 : // String name("ftaperture.im");
426 : // storeImg(name,*(ap.aperture));
427 : }
428 :
429 148 : }
430 0 : void VLACalcIlluminationConvFunc::regridAperture(CoordinateSystem& skyCS,
431 : IPosition& skyShape,
432 : TempImage<Complex>& uvGrid,
433 : const VisBuffer2 &vb,
434 : const Vector<Float>& paList,
435 : Bool doSquint, Int bandID)
436 : {
437 0 : CoordinateSystem skyCoords(skyCS);
438 :
439 : Float pa, Freq;
440 0 : if (bandID != -1) ap.band = bandID;
441 0 : AlwaysAssert(ap.band>=-1, AipsError);
442 0 : Vector<Double> chanFreq = vb.getFrequencies(0);
443 :
444 : const MSSpWindowColumns& spwCol =
445 0 : vb.subtableColumns().spectralWindow();
446 0 : ArrayColumn<Double> chanfreq = spwCol.chanFreq();
447 0 : ScalarColumn<Double> reffreq = spwCol.refFrequency();
448 : // Freq = sum(chanFreq)/chanFreq.nelements();
449 :
450 0 : Freq = max(chanfreq.getColumn());
451 0 : IPosition imsize(skyShape);
452 0 : CoordinateSystem uvCoords = makeUVCoords(skyCoords,imsize,Freq);
453 0 : uvGrid.setCoordinateInfo(uvCoords);
454 :
455 0 : Int index = uvCoords.findCoordinate(Coordinate::LINEAR);
456 0 : LinearCoordinate lc=uvCoords.linearCoordinate(index);
457 0 : Vector<Double> incr = lc.increment();
458 0 : Double Lambda = C::c/Freq;
459 0 : ap.freq = Freq/1E9;
460 :
461 0 : ap.nx = skyShape(0); ap.ny = skyShape(1);
462 0 : IPosition apertureShape(ap.aperture->shape());
463 0 : apertureShape(0) = ap.nx; apertureShape(1) = ap.ny;
464 0 : ap.aperture->resize(apertureShape);
465 :
466 0 : TempImage<Complex> tmpAperture;tmpAperture.resize(apertureShape);
467 0 : if (maximumCacheSize() > 0) tmpAperture.setMaximumCacheSize(maximumCacheSize());
468 0 : tmpAperture.set(0.0);
469 0 : ap.dx = abs(incr(0)*Lambda); ap.dy = abs(incr(1)*Lambda);
470 : // cout << ap.dx << " " << incr(0) << endl;
471 : // ap.x0 = -(25.0/(2*ap.dx)+1)*ap.dx; ap.y0 = -(25.0/(2*ap.dy)+1)*ap.dy;
472 : // Following 3 lines go with the ANT tag in BeamCalc.cc
473 : // Double antRadius=BeamCalcGeometryDefaults[ap.band].Rant;
474 : // ap.x0 = -(antRadius/(ap.dx)+1)*ap.dx;
475 : // ap.y0 = -(antRadius/(ap.dy)+1)*ap.dy;
476 : // Following 2 lines go with the PIX tag in BeamCalc.cc
477 0 : ap.x0 = -(ap.nx/2)*ap.dx;
478 0 : ap.y0 = -(ap.ny/2)*ap.dy;
479 :
480 : //
481 : // Accumulate apertures for a list of PA
482 : //
483 0 : for(uInt ipa=0;ipa<paList.nelements();ipa++)
484 : {
485 0 : pa = paList[ipa];
486 :
487 0 : ap.pa=pa;
488 0 : ap.aperture->set(0.0);
489 0 : BeamCalc::Instance()->calculateAperture(&ap);
490 : //
491 : // Set the phase of the aperture function to zero if doSquint==F
492 : // Poln. axis indices
493 : // 0: RR, 1:RL, 2:LR, 3:LL
494 : // This is electic field. 0=> Poln R,
495 : // 1=> Leakage of R->L
496 : // 2=> Leakage of L->R
497 : // 3=> Poln L
498 : //
499 : // The squint is removed in the following code using
500 : // honest-to-god pixel indexing. If this is not the most
501 : // efficient method of doing this in AIPS++ (i.e. instead use
502 : // slices etc.), then this cost will show up in making the
503 : // average PB. Since this goes over each pixel of a full
504 : // stokes (poln. really) complex image, look here (also) for
505 : // optimization (if required).
506 : //
507 0 : if (!doSquint)
508 : {
509 0 : IPosition PolnRIndex(4,0,0,0,0), PolnLIndex(4,0,0,3,0);
510 0 : IPosition tndx(4,0,0,0,0);
511 0 : for(tndx(3)=0;tndx(3)<apertureShape(3);tndx(3)++) // The freq. axis
512 0 : for(tndx(2)=0;tndx(2)<apertureShape(2);tndx(2)++) // The Poln. axis
513 0 : for(tndx(1)=0;tndx(1)<apertureShape(1);tndx(1)++) // The spatial
514 0 : for(tndx(0)=0;tndx(0)<apertureShape(0);tndx(0)++) // axis.
515 : {
516 0 : PolnRIndex(0)=PolnLIndex(0)=tndx(0);
517 0 : PolnRIndex(1)=PolnLIndex(1)=tndx(1);
518 0 : Complex val, Rval, Lval;
519 : Float phase;
520 0 : val = ap.aperture->getAt(tndx);
521 0 : Rval = ap.aperture->getAt(PolnRIndex);
522 0 : Lval = ap.aperture->getAt(PolnLIndex);
523 0 : phase = arg(Rval); Rval=Complex(cos(phase),sin(phase));
524 0 : phase = arg(Lval); Lval=Complex(cos(phase),sin(phase));
525 :
526 0 : if (tndx(2)==0) ap.aperture->putAt(val*conj(Rval),tndx);
527 0 : else if (tndx(2)==1) ap.aperture->putAt(val*conj(Lval),tndx);
528 0 : else if (tndx(2)==2) ap.aperture->putAt(val*conj(Rval),tndx);
529 0 : else if (tndx(2)==3) ap.aperture->putAt(val*conj(Lval),tndx);
530 : }
531 0 : }
532 0 : tmpAperture += *(ap.aperture);
533 : }
534 0 : (ap.aperture)->copyData(tmpAperture);
535 0 : tmpAperture.resize(IPosition(1,1));//Release temp. store.
536 0 : Vector<Int> poln(4);
537 0 : poln(0) = Stokes::RR;
538 0 : poln(1) = Stokes::RL;
539 0 : poln(2) = Stokes::LR;
540 0 : poln(3) = Stokes::LL;
541 0 : StokesCoordinate polnCoord(poln);
542 0 : SpectralCoordinate spectralCoord(MFrequency::TOPO,Freq,1.0,0.0);
543 : // uvCoords.addCoordinate(dirCoord);
544 0 : index = uvCoords.findCoordinate(Coordinate::STOKES);
545 : //cerr << "STOKES index " << index << endl;
546 0 : uvCoords.replaceCoordinate(polnCoord,index);
547 0 : index = uvCoords.findCoordinate(Coordinate::SPECTRAL);
548 : //cerr << "Spectral index " << index << endl;
549 0 : uvCoords.replaceCoordinate(spectralCoord,index);
550 :
551 0 : ap.aperture->setCoordinateInfo(uvCoords);
552 : //if (doSquint==false)
553 : //{
554 : // String name("aperture.im");
555 : // storeImg(name,*(ap.aperture));
556 : //logIO << "The aperture has been written to aperture.im"<< LogIO::POST;
557 : //}
558 :
559 0 : ftAperture(*(ap.aperture));
560 : //String name("aperture.im");
561 : //storeImg(name,*(ap.aperture));
562 : //logIO << "The fourier transform of the aperture has been written to ftaperture.im"<< LogIO::POST;
563 0 : }
564 :
565 :
566 :
567 148 : void VLACalcIlluminationConvFunc::fillPB(ImageInterface<Complex>& inImg,
568 : ImageInterface<Complex>& outImg,
569 : Bool Square)
570 : {
571 148 : IPosition imsize(outImg.shape());
572 148 : IPosition ndx(outImg.shape());
573 148 : IPosition inShape(inImg.shape()),inNdx;
574 148 : Vector<Int> inStokes,outStokes;
575 : Int index,s,index1;
576 : //cerr << "Complex IMAGENAME " << outImg.name() << " " << max(inImg.get()) << " " << min(inImg.get()) << endl;
577 : // Timer tim;
578 : // tim.mark();
579 148 : index = inImg.coordinates().findCoordinate(Coordinate::STOKES);
580 148 : inStokes = inImg.coordinates().stokesCoordinate(index).stokes();
581 148 : index = outImg.coordinates().findCoordinate(Coordinate::STOKES);
582 148 : outStokes = outImg.coordinates().stokesCoordinate(index).stokes();
583 148 : index = outImg.coordinates().findCoordinate(Coordinate::SPECTRAL);
584 148 : index1 = inImg.coordinates().findCoordinate(Coordinate::SPECTRAL);
585 148 : SpectralCoordinate inSpectralCoords = inImg.coordinates().spectralCoordinate(index1);
586 148 : CoordinateSystem outCS = outImg.coordinates();
587 148 : outCS.replaceCoordinate(inSpectralCoords,index);
588 148 : outImg.setCoordinateInfo(outCS);
589 : //tim.show("fillPB::CSStuff:");
590 148 : ndx(3)=0;
591 : // #ifdef HAS_OMP
592 : // Int Nth=max(omp_get_max_threads()-2,1);
593 : // #endif
594 : //tim.mark();
595 296 : for(ndx(2)=0;ndx(2)<imsize(2);ndx(2)++) // The poln axes
596 : {
597 148 : for(s=0;s<inShape(2);s++) if (inStokes(s) == outStokes(ndx(2))) break;
598 :
599 174228 : for(ndx(0)=0;ndx(0)<imsize(0);ndx(0)++)
600 : //#pragma omp parallel default(none) firstprivate(s,iy) shared(twoPiW,convSize) num_threads(Nth)
601 : {
602 : //#pragma omp for
603 290629632 : for(ndx(1)=0;ndx(1)<imsize(1);ndx(1)++)
604 : {
605 290455552 : Complex cval;
606 290455552 : inNdx = ndx; inNdx(2)=s;
607 290455552 : cval = inImg.getAt(inNdx);
608 290455552 : if (Square) cval = cval*conj(cval);
609 : // cerr << "Complex " << ndx << " val " << cval << endl;
610 290455552 : outImg.putAt(cval*outImg.getAt(ndx),ndx);
611 : }
612 : }
613 : }
614 : //tim.show("fillPB: ");
615 148 : }
616 :
617 0 : void VLACalcIlluminationConvFunc::fillPB(ImageInterface<Complex>& inImg,
618 : ImageInterface<Float>& outImg,
619 : Bool Square)
620 : {
621 :
622 : //cerr << "REAL IMAGENAME " << outImg.name() << " " << max(inImg.get()) << " " << min(inImg.get()) << endl;
623 0 : IPosition imsize(outImg.shape());
624 0 : IPosition ndx(outImg.shape());
625 0 : IPosition inShape(inImg.shape()),inNdx;
626 :
627 0 : Vector<Int> inStokes,outStokes;
628 : Int index,s;
629 :
630 : // Timer tim;
631 : // tim.mark();
632 0 : index = inImg.coordinates().findCoordinate(Coordinate::STOKES);
633 0 : inStokes = inImg.coordinates().stokesCoordinate(index).stokes();
634 0 : index = outImg.coordinates().findCoordinate(Coordinate::STOKES);
635 0 : outStokes = outImg.coordinates().stokesCoordinate(index).stokes();
636 0 : ndx(3)=0;
637 : //tim.show("fillPB::CSStuff2:");
638 :
639 : //tim.mark();
640 0 : for(ndx(2)=0;ndx(2)<imsize(2);ndx(2)++) // The poln axes
641 : {
642 0 : if (outStokes(ndx(2)) == Stokes::I)
643 : {
644 : //
645 : // Fill the outImg wiht (inImg(RR)+inImage(LL))/2
646 : //
647 0 : for(ndx(0)=0;ndx(0)<imsize(0);ndx(0)++)
648 0 : for(ndx(1)=0;ndx(1)<imsize(1);ndx(1)++)
649 : {
650 0 : Complex cval;
651 0 : inNdx = ndx;
652 0 : inNdx(2)=0; cval = inImg.getAt(inNdx);
653 0 : inNdx(2)=3; cval += inImg.getAt(inNdx);
654 0 : cval/2;
655 :
656 0 : if (Square) cval = cval*conj(cval);
657 : //cerr << "I " << ndx << " val " << cval << endl;
658 0 : outImg.putAt(abs(cval*outImg.getAt(ndx)),ndx);
659 : }
660 : }
661 0 : else if ((outStokes(ndx(2)) == Stokes::RR) ||
662 0 : (outStokes(ndx(2)) == Stokes::RL) ||
663 0 : (outStokes(ndx(2)) == Stokes::LR) ||
664 0 : (outStokes(ndx(2)) == Stokes::LL))
665 : {
666 0 : for(s=0;s<inShape(2);s++) if (inStokes(s) == outStokes(ndx(2))) break;
667 :
668 0 : for(ndx(0)=0;ndx(0)<imsize(0);ndx(0)++)
669 0 : for(ndx(1)=0;ndx(1)<imsize(1);ndx(1)++)
670 : {
671 0 : Complex cval;
672 0 : inNdx = ndx; inNdx(2)=s;
673 0 : cval = inImg.getAt(inNdx);
674 0 : if (Square) cval = cval*conj(cval);
675 0 : cerr << "Corr "<< ndx << " val" << cval << endl;
676 0 : outImg.putAt(abs(cval*outImg.getAt(ndx)),ndx);
677 : }
678 : }
679 0 : else throw(AipsError("Unsupported Stokes found in VLACalcIlluminationConvFunc::fillPB."));
680 : }
681 : //tim.show("fillPB2:");
682 0 : }
683 :
684 : /*
685 : void VLACalcIlluminationConvFunc::fillPB(ImageInterface<Complex>& inImg,
686 : ImageInterface<Float>& outImg)
687 : {
688 : IPosition imsize(outImg.shape());
689 : IPosition ndx(outImg.shape());
690 : IPosition inShape(inImg.shape()), inNdx;
691 : Vector<Int> inStokes,outStokes;
692 : Int index;
693 : index = inImg.coordinates().findCoordinate(Coordinate::STOKES);
694 : inStokes = inImg.coordinates().stokesCoordinate(index).stokes();
695 : index = outImg.coordinates().findCoordinate(Coordinate::STOKES);
696 : outStokes = outImg.coordinates().stokesCoordinate(index).stokes();
697 :
698 : ndx(3)=0;
699 : for(ndx(0)=0;ndx(0)<imsize(0);ndx(0)++)
700 : {
701 :
702 : for(ndx(1)=0;ndx(1)<imsize(1);ndx(1)++)
703 : {
704 : //
705 : // Average along the polarization axes and fillin the
706 : // amp. of the average in the output image.
707 : //
708 : Complex cval=0.0;
709 :
710 : ndx(2)=0;cval = inImg.getAt(ndx);
711 : for(ndx(2)=0;ndx(2)<imsize(2);ndx(2)++) // The poln axes
712 : outImg.putAt(abs(cval*outImg.getAt(ndx)),ndx);
713 : }
714 : }
715 : }
716 :
717 : */
718 :
719 : // void VLACalcIlluminationConvFunc::ftAperture(TempImage<Complex>& uvgrid, Bool makeMueller)
720 148 : void VLACalcIlluminationConvFunc::ftAperture(ImageInterface<Complex>& uvgrid, Int muellerTerm)
721 : {
722 : //
723 : // Make SkyJones
724 : //
725 : // {
726 : //String name("uvgrid.im");
727 : //storeImg(name,uvgrid);
728 : // }
729 148 : LatticeFFT::cfft2d(uvgrid);
730 : //cerr <<"FTap mueller = " << muellerTerm << " " << max(uvgrid.get()) << " " << min(uvgrid.get()) << " shape " << uvgrid.shape() << endl;
731 : // {
732 : // Int index = uvgrid.coordinates().findCoordinate(Coordinate::STOKES);
733 : // Int inStokes = uvgrid.coordinates().stokesCoordinate(index).stokes()(0);
734 : // IPosition shape = uvgrid.shape();
735 : //ostringstream tt;
736 : //String name("ftuvgrid.im");
737 : //tt << name << "_"<< inStokes;
738 : //storeImg(String(tt),uvgrid);
739 : // }
740 : // {
741 : //String name("ftuvgrid.im");
742 : //storeImg(name,uvgrid);
743 : // }
744 :
745 : // Now make SkyMuller
746 : //
747 : //skyMuller(uvgrid);
748 148 : Int tempmueller=-1; // Set to -1 to use sanjay's code which pass the regressions.
749 148 : tempmueller=0; // Full-pol. imaging is requested, use PJ's code. Else use SB's.
750 :
751 : //if (uvgrid.shape()(3) > 2) tempmueller=0;
752 :
753 148 : if (tempmueller==0)
754 148 : skyMuller(uvgrid,muellerTerm);
755 : else
756 0 : skyMuller(uvgrid,-1);
757 : //cerr <<"FTap post skymuller " << max(uvgrid.get()) << " " << min(uvgrid.get()) << endl;
758 148 : }
759 :
760 0 : void VLACalcIlluminationConvFunc::loadFromImage(String& /*fileName*/)
761 : {
762 0 : throw(AipsError("VLACalcIlluminationConvFunc::loadFromImage() not yet supported."));
763 : };
764 :
765 0 : void VLACalcIlluminationConvFunc::skyMuller(Array<Complex>& buf,
766 : const IPosition& shape,
767 : const Int& inStokes)
768 : {
769 0 : Array<Complex> tmp;
770 0 : IPosition t(4,0,0,0,0),n0(4,0,0,0,0),n1(4,0,0,0,0);
771 : Float peak;
772 0 : peak=0;
773 0 : for(t(2)=0;t(2)<shape(2);t(2)++)
774 0 : for(t(1)=0;t(1)<shape(1);t(1)++)
775 0 : for(t(0)=0;t(0)<shape(0);t(0)++)
776 0 : if (abs(buf(t)) > peak) peak = abs(buf(t));
777 0 : if (peak > 1E-8)
778 0 : for(t(3)=0;t(3)<shape(3);t(3)++) // Freq axis
779 0 : for(t(2)=0;t(2)<shape(2);t(2)++) // Poln axis
780 0 : for(t(1)=0;t(1)<shape(1);t(1)++) // y axis
781 0 : for(t(0)=0;t(0)<shape(0);t(0)++) // X axis
782 0 : buf(t) = buf(t)/peak;
783 : // {
784 : // skyJones.put(buf);
785 : // String name("skyjones.im");
786 : // storeImg(name,skyJones);
787 : // }
788 :
789 0 : tmp.assign(buf);
790 :
791 0 : t(0)=t(1)=t(2)=t(3)=0;
792 :
793 0 : if ((inStokes == Stokes::RR) || (inStokes == Stokes::LL))
794 : {
795 0 : t(2)=0;n0(2)=0;n1(2)=0; //RR
796 0 : for( n0(0)=n1(0)=t(0)=0;n0(0)<shape(0);n0(0)++,n1(0)++,t(0)++)
797 0 : for(n0(1)=n1(1)=t(1)=0;n0(1)<shape(1);n0(1)++,n1(1)++,t(1)++)
798 0 : buf(t) = (tmp(n0)*conj(tmp(n1)));
799 : }
800 :
801 0 : if ((inStokes == Stokes::LR) || (inStokes == Stokes::RL))
802 : {
803 0 : t(2)=0;n0(3)=1;n1(2)=0; //LR
804 0 : for( n0(0)=n1(0)=t(0)=0;n0(0)<shape(0);n0(0)++,n1(0)++,t(0)++)
805 0 : for(n0(1)=n1(1)=t(1)=0;n0(1)<shape(1);n0(1)++,n1(1)++,t(1)++)
806 0 : buf(t) = (tmp(n0)*conj(tmp(n1)));
807 :
808 : }
809 :
810 : // if (inStokes == Stokes::RL)
811 : // {
812 : // t(2)=0;n0(2)=0;n1(3)=1; //LR
813 : // for( n0(0)=n1(0)=t(0)=0;n0(0)<shape(0);n0(0)++,n1(0)++,t(0)++)
814 : // for(n0(1)=n1(1)=t(1)=0;n0(1)<shape(1);n0(1)++,n1(1)++,t(1)++)
815 : // buf(t) = (tmp(n0)*conj(tmp(n1)));
816 :
817 : // }
818 : //{
819 : //skyJones.put(buf);
820 : // ostringstream tt;
821 : //String name("skyjones.im");
822 : //tt << name << "_" << inStokes;
823 : //storeImg(tt.string(),skyJones);
824 : //}
825 :
826 : //cout<<"Regular skyjones \n";
827 0 : }
828 :
829 : // void VLACalcIlluminationConvFunc::skyMuller(ImageInterface<Complex>& skyJones)
830 : // {
831 : // Int index = skyJones.coordinates().findCoordinate(Coordinate::STOKES);
832 : // Int inStokes = skyJones.coordinates().stokesCoordinate(index).stokes()(0);
833 : // Array<Complex> buf=skyJones.get();
834 : //Vector<Int> shape(buf.shape().asVector());
835 : // IPosition shape=skyJones.shape();
836 : // skyMuller(buf,shape, inStokes);
837 : // skyJones.put(buf);
838 : // }
839 :
840 148 : void VLACalcIlluminationConvFunc::skyMuller(ImageInterface<Complex>& skyJones, Int muellerTerm)
841 : {
842 148 : Int index = skyJones.coordinates().findCoordinate(Coordinate::STOKES);
843 : //cerr << "STOKES index in skymuller " << index << " shape " << skyJones.shape() << endl;
844 148 : Int inStokes = skyJones.coordinates().stokesCoordinate(index).stokes()(0);
845 : //PagedImage<Complex> lala(skyJones.shape(), skyJones.coordinates(), "BEFORESKYMULL.im");
846 : //lala.copyData(skyJones);
847 148 : Array<Complex> buf=skyJones.get(),tmp;
848 148 : IPosition shape=skyJones.shape();
849 148 : if(muellerTerm == -1)
850 : {
851 : //cout<<"####Temp - bypassing the full mueller convolution function generation \n";
852 0 : skyMuller(buf,shape,inStokes);
853 0 : skyJones.put(buf);
854 : }
855 : else
856 : {
857 : //cout<<"####Temp - Proceeding with IQUV convolution function generation\n";
858 148 : Array<Complex> M0,M1,M2,M3;
859 : //Vector<Int> shape(buf.shape().asVector());
860 : // IPosition shape=skyJones.shape();
861 : // cout<<"The shape of sky Jones matrix is "<<shape<<"\n";
862 148 : Array<Complex> Jp,Jq,Jpq,Jqp;
863 : // skyMuller(buf,shape, inStokes);
864 :
865 : // if (muellerTerm == 5)
866 : // {
867 : // skyJones.put(buf);
868 : // String name("skyjones5.im");
869 : // storeImg(name,skyJones);
870 : // }
871 : // Int index = skyJones.coordinates().findCoordinate(Coordinate::STOKES);
872 : // Int inStokes = skyJones.coordinates().stokesCoordinate(index).stokes()(0);
873 : // Array<Complex> buf=skyJones.get(),tmp;
874 :
875 148 : IPosition t(4,0,0,0,0),n0(4,0,0,0,0),n1(4,0,0,0,0);
876 :
877 148 : skyJones.put(buf); ////?????????????
878 : // ostringstream tt;
879 : // String name("skyjones.im");
880 : // tt << name << "_"<< inStokes;
881 : // storeImg(String(tt),skyJones);
882 : // skyJones.put(buf);
883 : // cout<<"Finished writing the initial sky jones \n";
884 : // exit(0);
885 : // skyMuller(buf,shape, inStokes);
886 :
887 148 : Float Normalizesq,pqscale=100.0;
888 : Int midx,midy;
889 148 : Normalizesq=0;
890 148 : midx = shape(0)/2; // This gives the central pixel(not the peak) of RR and LL.
891 148 : midy = shape(1)/2; // This normalization scheme was tested to be correct in the python code.
892 :
893 148 : tmp=buf;
894 296 : for(t(2)=0;t(2)<shape(2);t(2)++) // Start with poln axis as we want to loop in freq for the moment
895 740 : for(t(3)=0;t(3)<shape(3);t(3)++) // The freq axis each one of 4 chans contains a jones element
896 696912 : for(t(1)=0;t(1)<shape(1);t(1)++)
897 1162518528 : for(t(0)=0;t(0)<shape(0);t(0)++)
898 1161822208 : if((t(0)== midx)&&(t(1)==midy)) //?????Seriously quite a bogus looping..this is not a tape !
899 592 : Normalizesq = Normalizesq + abs(buf(t)*buf(t))/2.0; // This needs to be changed so that Normalizesq stays complex
900 :
901 :
902 296 : for(t(2)=0;t(2)<shape(2);t(2)++)
903 740 : for(t(3)=0;t(3)<shape(3);t(3)++)
904 696912 : for(t(1)=0;t(1)<shape(1);t(1)++)
905 1162518528 : for(t(0)=0;t(0)<shape(0);t(0)++)
906 1161822208 : if((t(3)==1)||(t(3)==2))
907 580911104 : tmp(t)=pqscale*conj(tmp(t)/sqrt(Normalizesq)); // This is the crosshand scale factor to be determined and hardcoded, currently is arbitrary.
908 :
909 : else
910 580911104 : tmp(t)=conj(tmp(t)/sqrt(Normalizesq));
911 :
912 : // cout<<"The Jones Matrix has been normalized using:"<< sqrt(Normalizesq)<<"\n";
913 148 : skyJones.put(tmp); //Why ??
914 : //cerr << "after weird normalization " << max(tmp) << " " << min(tmp) << endl;
915 : // ostringstream tt1;
916 : // String name1("skyjones_normalized_conj.im");
917 : // tt1 << name1 << "_"<< inStokes;
918 : // storeImg(String(tt1),skyJones);
919 : // // exit(0);
920 : // // skyJones.put(tmp);
921 : // cout<<"Finished writing the normalized conjugate sky jones \n";
922 :
923 : // cout<<"Begining the compute of Mueller Matrix term images \n";
924 :
925 148 : IPosition sliceStart0(4,0,0,0,0),sliceStart1(4,0,0,0,1),sliceLength0(4,shape(0),shape(1),1,1),sliceLength1(4,shape(0),shape(1),1,1);
926 148 : IPosition sliceStart2(4,0,0,0,2),sliceStart3(4,0,0,0,3),sliceLength2(4,shape(0),shape(1),1,1),sliceLength3(4,shape(0),shape(1),1,1);
927 148 : Slicer s0(sliceStart0,sliceLength0),s1(sliceStart1,sliceLength1);
928 148 : Slicer s2(sliceStart2,sliceLength2),s3(sliceStart3,sliceLength3);
929 :
930 : // For the sake of pixel math we are resizing the four initial jones buffers.
931 148 : IPosition shp=buf.shape();
932 : // shp = buf.shape();
933 : // shp(2)=shp(3);
934 : // shp(3)=1;
935 148 : shp(2)=shp(3)=1;
936 148 : Jp.resize(shp);
937 148 : Jpq.resize(shp);
938 148 : Jqp.resize(shp);
939 148 : Jq.resize(shp);
940 : // cout<<"The Jones buffer shape is :"<< shp <<"\n";
941 :
942 148 : Jp(s0)=tmp(s0);Jpq(s0)=tmp(s1);
943 148 : Jqp(s0)=tmp(s2);Jq(s0)=tmp(s3);
944 148 : M0=M1=M2=M3=tmp; //what an unnecessary mem copy
945 : // We will initialize the Mueller rows to zero as per need and then slice and return with only the first slice with written values
946 :
947 148 : M0(s0)=Jp*conj(Jp); M0(s1)=Jp*conj(Jpq); M0(s2)=Jpq*conj(Jp); M0(s3)=Jpq*conj(Jpq);
948 148 : M1(s0)=Jp*conj(Jqp); M1(s1)=Jp*conj(Jq); M1(s2)=Jpq*conj(Jqp); M1(s3)=Jpq*conj(Jq);
949 148 : M2(s0)=Jqp*conj(Jp); M2(s1)=Jqp*conj(Jpq); M2(s2)=Jq*conj(Jq); M2(s3)=Jq*conj(Jpq);
950 148 : M3(s0)=Jqp*conj(Jqp); M3(s1)=Jqp*conj(Jq); M3(s2)=Jq*conj(Jqp); M3(s3)=Jq*conj(Jq);
951 : // cout<<"Slicing Complete"<<"\n";
952 : // cout<<"Mueller row selection in place is :" << muellerTerm << "\n";
953 148 : if (muellerTerm <= 3)
954 : {
955 74 : M0=0; // ??????? all these multiplication up there is just to warm the CPU !
956 74 : if (muellerTerm==0) {
957 74 : M0(s0)=Jp*conj(Jp);
958 : }
959 0 : else if (muellerTerm==1){
960 0 : M0(s0)=Jp*conj(Jpq);
961 : }
962 0 : else if (muellerTerm==2){
963 0 : M0(s0)=Jpq*conj(Jp);
964 : }
965 : else {
966 0 : M0(s0)=Jpq*conj(Jpq);
967 : }
968 : //cerr << "after M0 " << max(M0) << " " << min(M0) << endl;
969 74 : skyJones.put(M0);
970 : //String name2("M0.im");
971 : //storeImg(name2,skyJones);
972 : // cout<<"Writing M0 to disk, muellerTerm : "<< muellerTerm <<"\n";
973 74 : M0.resize();
974 : }
975 74 : else if ((4<=muellerTerm)&&(muellerTerm<8))
976 : {
977 0 : M1=0;
978 0 : if (muellerTerm==4) {
979 0 : M1(s0)=Jp*conj(Jqp);
980 : }
981 0 : else if (muellerTerm==5){
982 0 : M1(s0)=Jp*conj(Jq);
983 : }
984 0 : else if (muellerTerm==6){
985 0 : M1(s0)=Jpq*conj(Jqp);
986 : }
987 : else {
988 0 : M1(s0)=Jpq*conj(Jq);
989 : }
990 0 : skyJones.put(M1);
991 : // String name3("M1.im");
992 : // storeImg(name3,skyJones);
993 : // cout<<"Writing M1 to disk, muellerTerm : "<< muellerTerm <<"\n";
994 0 : M1.resize();
995 : }
996 74 : else if ((8<=muellerTerm)&&(muellerTerm<12))
997 : {
998 0 : M2=0;
999 0 : if (muellerTerm==8) {
1000 0 : M2(s0)=Jqp*conj(Jp);
1001 : }
1002 0 : else if (muellerTerm==9){
1003 0 : M2(s0)=Jqp*conj(Jpq);
1004 : }
1005 0 : else if (muellerTerm==10){
1006 0 : M2(s0)=Jq*conj(Jp);
1007 : }
1008 : else {
1009 0 : M2(s0)=Jq*conj(Jpq);
1010 : }
1011 0 : skyJones.put(M2);
1012 : // String name4("M2.im");
1013 : // storeImg(name4,skyJones);
1014 : // cout<<"Writing M2 to disk, muellerTerm : "<< muellerTerm <<"\n";
1015 0 : M2.resize();
1016 : }
1017 74 : else if ((12<=muellerTerm)&&(muellerTerm<16))
1018 : {
1019 74 : M3=0;
1020 74 : if (muellerTerm==12) {
1021 0 : M3(s0)=Jqp*conj(Jqp);
1022 : }
1023 74 : else if (muellerTerm==13){
1024 0 : M3(s0)=Jqp*conj(Jq);
1025 : }
1026 74 : else if (muellerTerm==14){
1027 0 : M3(s0)=Jq*conj(Jqp);
1028 : }
1029 : else {
1030 74 : M3(s0)=Jq*conj(Jq);
1031 : }
1032 74 : skyJones.put(M3);
1033 : // String name5("M3.im");
1034 : // storeImg(name5,skyJones);
1035 : // cout<<"Writing M3 to disk, muellerTerm : "<< muellerTerm <<"\n";
1036 74 : M3.resize();
1037 : }
1038 : // cout<<"Mueller Matrix row computation complete \n";
1039 148 : }
1040 : // skyJones.put(buf);
1041 :
1042 148 : }
1043 :
1044 :
1045 :
1046 0 : void VLACalcIlluminationConvFunc::makeFullJones(ImageInterface<Complex>& /*pbImage*/,
1047 : const VisBuffer2& /*vb*/,
1048 : Bool /*doSquint*/,
1049 : Int /*bandID*/,
1050 : Double /*freqVal*/)
1051 0 : {}
1052 :
1053 : };
1054 : };
|
