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