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