Line data Source code
1 : //# PBMath2DImage.cc: Implementation for PBMath2DImage
2 : //# Copyright (C) 1996,1997,1998,1999,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be adressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //#
27 : //# $Id$
28 :
29 : #include <casacore/casa/aips.h>
30 : #include <casacore/casa/Arrays/ArrayMath.h>
31 : #include <casacore/casa/Exceptions/Error.h>
32 : #include <casacore/casa/BasicSL/Complex.h>
33 : #include <casacore/casa/BasicMath/Math.h>
34 : #include <casacore/scimath/Mathematics/FFTServer.h>
35 : #include <synthesis/TransformMachines/PBMath2DImage.h>
36 : #include <casacore/images/Images/PagedImage.h>
37 : #include <casacore/images/Images/TempImage.h>
38 : #include <casacore/images/Images/ImageRegrid.h>
39 : #include <casacore/images/Images/ImageSummary.h>
40 : #include <msvis/MSVis/StokesVector.h>
41 : #include <synthesis/TransformMachines/StokesImageUtil.h>
42 : #include <casacore/lattices/Lattices/LatticeStepper.h>
43 : #include <casacore/lattices/Lattices/LatticeIterator.h>
44 : #include <casacore/lattices/LEL/LatticeExpr.h>
45 : #include <casacore/lattices/LEL/LatticeExprNode.h>
46 : #include <casacore/casa/Utilities/Assert.h>
47 : #include <components/ComponentModels/ComponentType.h>
48 : #include <casacore/casa/Quanta.h>
49 : #include <casacore/measures/Measures.h>
50 : #ifdef _OPENMP
51 : #include <omp.h>
52 : #endif
53 :
54 : using namespace casacore;
55 : typedef unsigned long long ooLong;
56 : namespace casa {
57 :
58 0 : PBMath2DImage::PBMath2DImage(ImageInterface<Float>& reJonesImage):
59 : PBMath2D(), reJonesImage_p(0), reRegridJonesImage_p(0),
60 : imJonesImage_p(0), imRegridJonesImage_p(0),
61 : incrementsReJones_p(0), incrementsImJones_p(0),
62 : referencePixelReJones_p(0), referencePixelImJones_p(0),
63 0 : pa_p(0.0)
64 : {
65 0 : LogIO os(LogOrigin("PBMath2DImage", "PBMath2DImage"));
66 :
67 0 : os << "Using two-dimensional purely real image model for antenna voltage pattern" << LogIO::POST;
68 :
69 0 : reJonesImage_p = new TempImage<Float>(reJonesImage.shape(), reJonesImage.coordinates());
70 0 : reJonesImage_p->copyData(reJonesImage);
71 0 : incrementsReJones_p=new Vector<Double>(reJonesImage_p->coordinates().directionCoordinate(0).increment());
72 0 : referencePixelReJones_p=new Vector<Double>(reJonesImage_p->coordinates().directionCoordinate(0).referencePixel());
73 :
74 0 : };
75 :
76 0 : PBMath2DImage::PBMath2DImage(ImageInterface<Float>& reJonesImage,
77 0 : ImageInterface<Float>& imJonesImage) :
78 : PBMath2D(), reJonesImage_p(0), reRegridJonesImage_p(0),
79 : imJonesImage_p(0), imRegridJonesImage_p(0),
80 : incrementsReJones_p(0), incrementsImJones_p(0),
81 : referencePixelReJones_p(0), referencePixelImJones_p(0),
82 0 : pa_p(0.0)
83 : {
84 :
85 0 : LogIO os(LogOrigin("PBMath2DImage", "PBMath2DImage"));
86 :
87 0 : os << "Using two-dimensional real and imaginary image models for antenna voltage pattern" << LogIO::POST;
88 :
89 0 : checkJonesCongruent(reJonesImage, imJonesImage);
90 :
91 : // Save images and useful information
92 0 : reJonesImage_p = new TempImage<Float>(reJonesImage.shape(), reJonesImage.coordinates());
93 0 : reJonesImage_p->copyData(reJonesImage);
94 0 : incrementsReJones_p=new Vector<Double>(reJonesImage_p->coordinates().directionCoordinate(0).increment());
95 0 : referencePixelReJones_p=new Vector<Double>(reJonesImage_p->coordinates().directionCoordinate(0).referencePixel());
96 :
97 0 : imJonesImage_p = new TempImage<Float>(imJonesImage.shape(), imJonesImage.coordinates());
98 0 : imJonesImage_p->copyData(imJonesImage);
99 0 : incrementsImJones_p=new Vector<Double>(imJonesImage_p->coordinates().directionCoordinate(0).increment());
100 0 : referencePixelImJones_p=new Vector<Double>(imJonesImage_p->coordinates().directionCoordinate(0).referencePixel());
101 0 : };
102 0 : PBMath2DImage::PBMath2DImage(const ImageInterface<Complex>& jonesImage) :
103 : PBMath2D(), reJonesImage_p(0), reRegridJonesImage_p(0),
104 : imJonesImage_p(0), imRegridJonesImage_p(0),
105 : incrementsReJones_p(0), incrementsImJones_p(0),
106 : referencePixelReJones_p(0), referencePixelImJones_p(0),
107 0 : pa_p(0.0){
108 0 : reJonesImage_p = new TempImage<Float>(jonesImage.shape(), jonesImage.coordinates());
109 0 : reJonesImage_p->copyData(LatticeExpr<float> (real(jonesImage)));
110 0 : imJonesImage_p = new TempImage<Float>(jonesImage.shape(), jonesImage.coordinates());
111 0 : imJonesImage_p->copyData(LatticeExpr<float> (imag(jonesImage)));
112 0 : incrementsReJones_p=new Vector<Double>(reJonesImage_p->coordinates().directionCoordinate(0).increment());
113 0 : referencePixelReJones_p=new Vector<Double>(reJonesImage_p->coordinates().directionCoordinate(0).referencePixel());
114 0 : incrementsImJones_p=new Vector<Double>(imJonesImage_p->coordinates().directionCoordinate(0).increment());
115 0 : referencePixelImJones_p=new Vector<Double>(imJonesImage_p->coordinates().directionCoordinate(0).referencePixel());
116 0 : }
117 :
118 :
119 :
120 0 : PBMath2DImage::~PBMath2DImage()
121 : {
122 0 : if(reJonesImage_p) delete reJonesImage_p; reJonesImage_p=0;
123 0 : if(imJonesImage_p) delete imJonesImage_p; imJonesImage_p=0;
124 0 : if(incrementsReJones_p) delete incrementsReJones_p; incrementsReJones_p=0;
125 0 : if(incrementsImJones_p) delete incrementsImJones_p; incrementsImJones_p=0;
126 0 : if(referencePixelReJones_p) delete referencePixelReJones_p; referencePixelReJones_p=0;
127 0 : if(referencePixelImJones_p) delete referencePixelImJones_p; referencePixelImJones_p=0;
128 0 : };
129 :
130 0 : PBMath2DImage& PBMath2DImage::operator=(const PBMath2DImage& other)
131 : {
132 0 : if (this == &other)
133 0 : return *this;
134 :
135 0 : PBMath2D::operator=(other);
136 0 : reJonesImage_p = other.reJonesImage_p;
137 0 : imJonesImage_p = other.imJonesImage_p;
138 0 : incrementsReJones_p = other.incrementsReJones_p;
139 0 : incrementsImJones_p = other.incrementsImJones_p;
140 0 : referencePixelReJones_p = other.referencePixelReJones_p;
141 0 : referencePixelImJones_p = other.referencePixelImJones_p;
142 :
143 0 : return *this;
144 : };
145 :
146 :
147 :
148 : void
149 0 : PBMath2DImage::summary(Int nValues)
150 : {
151 0 : PBMath2D::summary(nValues);
152 0 : LogIO os(LogOrigin("PBMath2DImage", "summary"));
153 :
154 : {
155 0 : os << "Original image of the real part of Jones matrix:" << LogIO::POST;
156 0 : ImageSummary<Float> is(*reJonesImage_p);
157 0 : is.list(os);
158 : }
159 :
160 0 : if(imJonesImage_p) {
161 : os << "Original image of the imaginary part of Jones matrix:"
162 0 : << LogIO::POST;
163 0 : ImageSummary<Float> is(*imJonesImage_p);
164 0 : is.list(os);
165 : }
166 :
167 0 : };
168 :
169 : // Apply the Jones matrices to the input cube (x,y,pol)
170 : // giving an output cube (x,y,pol). These two methods
171 : // should be moved to Fortran for optimization
172 :
173 : // Complex to Complex
174 : void
175 0 : PBMath2DImage::applyJones(const Array<Float>* reJones,
176 : const Array<Float>* imJones,
177 : const Array<Complex>& in,
178 : Array<Complex>& out,
179 : Vector<Int>& polmap,
180 : Bool /*inverse*/,
181 : Bool /*conjugate*/,
182 : Int ipower, // ie, 1=VP, 2=PB
183 : Float /*cutoff*/,
184 : Bool circular,
185 : Bool forward)
186 : {
187 :
188 0 : LogIO os(LogOrigin("PBMath2DImage", "applyJones"));
189 : // This should never be called
190 0 : if((ipower!=2)&&(ipower!=1)) {
191 : os << "Logic error - trying to apply illegal power of PB"
192 0 : << LogIO::EXCEPTION;
193 : }
194 :
195 0 : Int nx=in.shape()(0);
196 0 : Int ny=in.shape()(1);
197 0 : Int npol=in.shape()(2);
198 :
199 : // Loop through x, y coordinates of this cube
200 0 : Matrix<Complex> cmat(2,2);
201 0 : IPosition sp0(4, 0, 0, polmap(3), 0);
202 0 : IPosition sp1(4, 0, 0, polmap(2), 0);
203 0 : IPosition sp2(4, 0, 0, polmap(1), 0);
204 0 : IPosition sp3(4, 0, 0, polmap(0), 0);
205 0 : for (Int ix=0; ix<nx; ix++) {
206 0 : sp0(0)=ix;
207 0 : sp1(0)=ix;
208 0 : sp2(0)=ix;
209 0 : sp3(0)=ix;
210 :
211 0 : for (Int iy=0; iy<ny; iy++) {
212 :
213 0 : sp0(1)=iy;
214 0 : sp1(1)=iy;
215 0 : sp2(1)=iy;
216 0 : sp3(1)=iy;
217 :
218 : // E Jones for this pixel
219 0 : mjJones4 j4;
220 0 : if(imJones) {
221 0 : cmat(0,0)=Complex((*reJones)(sp0), (*imJones)(sp0));
222 0 : cmat(1,0)=Complex((*reJones)(sp1), (*imJones)(sp1));
223 0 : cmat(0,1)=Complex((*reJones)(sp2), (*imJones)(sp2));
224 0 : cmat(1,1)=Complex((*reJones)(sp3), (*imJones)(sp3));
225 : }
226 : else {
227 0 : cmat(0,0)=Complex((*reJones)(sp0), 0.0);
228 0 : cmat(1,0)=Complex((*reJones)(sp1), 0.0);
229 0 : cmat(0,1)=Complex((*reJones)(sp2), 0.0);
230 0 : cmat(1,1)=Complex((*reJones)(sp3), 0.0);
231 : }
232 0 : mjJones2 j2(cmat);
233 :
234 : // Make the relevant Jones matrix
235 : // E
236 0 : if(ipower==1) { // VP
237 0 : mjJones2 j2unit(Complex(1.0, 0.0));
238 0 : directProduct(j4, j2, j2unit);
239 : }
240 : // Primary beam = E . conj(E)
241 0 : else if(ipower==2) { // PB
242 : // Make the conjugate before constructing the
243 : // mjJones since otherwise reference semantics
244 : // get us
245 0 : mjJones2 j2conj(conj(cmat));
246 0 : directProduct(j4, j2, j2conj);
247 : }
248 :
249 : // Subtlety - we have to distinguish between applying
250 : // the PB and applying the adjoint. The former is needed
251 : // for predictions (sky->UV) and the latter is needed
252 : // for inversion (UV->sky). For circular polarization
253 : // this affects only the cross hand terms
254 0 : if(!forward) {
255 0 : j4=mjJones4(adjoint(j4.matrix()));
256 : }
257 :
258 : // Now apply the Jones matrix
259 0 : if(npol==1) {
260 0 : IPosition ip0(4, ix, iy, 0, 0);
261 0 : CStokesVector outCS(in(ip0), 0.0, 0.0, 0.0);
262 0 : outCS*=j4;
263 0 : out(ip0)=outCS(0);
264 : }
265 0 : else if(npol==2) {
266 0 : IPosition ip0(4, ix, iy, 0, 0);
267 0 : IPosition ip1(4, ix, iy, 1, 0);
268 0 : if(circular) {
269 0 : CStokesVector outCS(in(ip0), 0.0, 0.0, in(ip1));
270 0 : outCS*=j4;
271 0 : out(ip0)=outCS(0);
272 0 : out(ip1)=outCS(3);
273 : }
274 : else {
275 0 : CStokesVector outCS(in(ip0), in(ip1), 0.0, 0.0);
276 0 : outCS*=j4;
277 0 : out(ip0)=outCS(0);
278 0 : out(ip1)=outCS(1);
279 : }
280 : }
281 0 : else if(npol==4) {
282 0 : IPosition ip0(4, ix, iy, 0, 0);
283 0 : IPosition ip1(4, ix, iy, 1, 0);
284 0 : IPosition ip2(4, ix, iy, 2, 0);
285 0 : IPosition ip3(4, ix, iy, 3, 0);
286 0 : CStokesVector outCS(in(ip0), in(ip1), in(ip2), in(ip3));
287 0 : outCS*=j4;
288 0 : out(ip0)=outCS(0);
289 0 : out(ip1)=outCS(1);
290 0 : out(ip2)=outCS(2);
291 0 : out(ip3)=outCS(3);
292 : }
293 : }
294 : }
295 0 : }
296 :
297 :
298 : // Complex to Complex
299 : void
300 0 : PBMath2DImage::applyJonesFast(const Float*& reJones,
301 : const Float*& imJones,
302 : const Array<Complex>& in,
303 : Array<Complex>& out,
304 : Vector<Int>& polmap,
305 : Bool /*inverse*/,
306 : Bool /*conjugate*/,
307 : Int ipower, // ie, 1=VP, 2=PB
308 : Float /*cutoff*/,
309 : Bool circular,
310 : Bool forward)
311 : {
312 :
313 0 : LogIO os(LogOrigin("PBMath2DImage", "applyJones"));
314 : // This should never be called
315 0 : if((ipower!=2)&&(ipower!=1)) {
316 : os << "Logic error - trying to apply illegal power of PB"
317 0 : << LogIO::EXCEPTION;
318 : }
319 :
320 : Bool delpolmap;
321 0 : const Int *polmap1=polmap.getStorage(delpolmap);
322 0 : const Float* reJones1=reJones;
323 0 : const Float* imJones1=imJones;
324 :
325 0 : Int nx=in.shape()(0);
326 0 : Int ny=in.shape()(1);
327 0 : Int npol=in.shape()(2);
328 0 : Bool lala=false;
329 0 : Float laloo=0.0;
330 :
331 : Bool delin, delout;
332 0 : Complex *outstor=out.getStorage(delout);
333 0 : const Complex *instor=in.getStorage(delin);
334 : //Int ind0, ind1, ind2, ind3;
335 : // Loop through x, y coordinates of this cube
336 0 : Matrix<Complex> cmat(2,2);
337 0 : IPosition sp0(4, 0, 0, polmap(3), 0);
338 0 : IPosition sp1(4, 0, 0, polmap(2), 0);
339 0 : IPosition sp2(4, 0, 0, polmap(1), 0);
340 0 : IPosition sp3(4, 0, 0, polmap(0), 0);
341 :
342 :
343 0 : #pragma omp parallel default(none) shared(outstor) firstprivate(reJones1, imJones1 , instor, polmap1, lala, laloo, ipower, circular, forward, nx,ny, npol)
344 :
345 : {
346 : #pragma omp for
347 : for (Int ix=0; ix<nx; ix++) {
348 : //sp0(0)=ix;
349 : //sp1(0)=ix;
350 : //sp2(0)=ix;
351 : //sp3(0)=ix;
352 :
353 : applyJonesFastX(reJones1, imJones1, instor, outstor, polmap1,
354 : lala, lala,ipower, // ie, 1=VP, 2=PB
355 : laloo,circular,forward,ix, nx, ny, npol);
356 : /* for (Int iy=0; iy<ny; iy++) {
357 :
358 : //sp0(1)=iy;
359 : //sp1(1)=iy;
360 : //sp2(1)=iy;
361 : //sp3(1)=iy;
362 :
363 : ind0=ix+nx*iy +(nx*ny)*polmap(3);
364 : ind1=ix+nx*iy +(nx*ny)*polmap(2);
365 : ind2=ix+nx*iy +(nx*ny)*polmap(1);
366 : ind3=ix+nx*iy +(nx*ny)*polmap(0);
367 : // E Jones for this pixel
368 : mjJones4 j4;
369 : if(imJones) {
370 : cmat(0,0)=Complex(reJones[ind0], imJones[ind0]);
371 : cmat(1,0)=Complex(reJones[ind1], imJones[ind1]);
372 : cmat(0,1)=Complex(reJones[ind2], imJones[ind2]);
373 : cmat(1,1)=Complex(reJones[ind3], imJones[ind3]);
374 : }
375 : else {
376 : cmat(0,0)=Complex(reJones[ind0], 0.0);
377 : cmat(1,0)=Complex(reJones[ind1], 0.0);
378 : cmat(0,1)=Complex(reJones[ind2], 0.0);
379 : cmat(1,1)=Complex(reJones[ind3], 0.0);
380 : }
381 : mjJones2 j2(cmat);
382 :
383 : // Make the relevant Jones matrix
384 : // E
385 : if(ipower==1) { // VP
386 : mjJones2 j2unit(Complex(1.0, 0.0));
387 : directProduct(j4, j2, j2unit);
388 : }
389 : // Primary beam = E . conj(E)
390 : else if(ipower==2) { // PB
391 : // Make the conjugate before constructing the
392 : // mjJones since otherwise reference semantics
393 : // get us
394 : mjJones2 j2conj(conj(cmat));
395 : directProduct(j4, j2, j2conj);
396 : }
397 :
398 : // Subtlety - we have to distinguish between applying
399 : // the PB and applying the adjoint. The former is needed
400 : // for predictions (sky->UV) and the latter is needed
401 : // for inversion (UV->sky). For circular polarization
402 : // this affects only the cross hand terms
403 : if(!forward) {
404 : j4=mjJones4(adjoint(j4.matrix()));
405 : }
406 :
407 : // Now apply the Jones matrix
408 : if(npol==1) {
409 : //IPosition ip0(4, ix, iy, 0, 0);
410 : ind0=ix+iy*nx;
411 : CStokesVector outCS(instor[ind0], 0.0, 0.0, 0.0);
412 : outCS*=j4;
413 : outstor[ind0]=outCS(0);
414 : }
415 : else if(npol==2) {
416 : //IPosition ip0(4, ix, iy, 0, 0);
417 : //IPosition ip1(4, ix, iy, 1, 0);
418 : ind0=ix+iy*nx;
419 : ind1=ix+iy*nx+nx*ny;
420 : if(circular) {
421 : CStokesVector outCS(instor[ind0], 0.0, 0.0, instor[ind1]);
422 : outCS*=j4;
423 : outstor[ind0]=outCS(0);
424 : outstor[ind1]=outCS(3);
425 : }
426 : else {
427 : CStokesVector outCS(instor[ind0], instor[ind1], 0.0, 0.0);
428 : outCS*=j4;
429 : outstor[ind0]=outCS(0);
430 : outstor[ind1]=outCS(1);
431 : }
432 : }
433 : else if(npol==4) {
434 : //IPosition ip0(4, ix, iy, 0, 0);
435 : //IPosition ip1(4, ix, iy, 1, 0);
436 : //IPosition ip2(4, ix, iy, 2, 0);
437 : //IPosition ip3(4, ix, iy, 3, 0);
438 : ind0=ix+iy*nx;
439 : ind1=ix+iy*nx+nx*ny;
440 : ind2=ix+iy*nx+2*nx*ny;
441 : ind3=ix+iy*nx+3*nx*ny;
442 : CStokesVector outCS(instor[ind0], instor[ind1], instor[ind2], instor[ind3]);
443 : outCS*=j4;
444 : outstor[ind0]=outCS(0);
445 : outstor[ind1]=outCS(1);
446 : outstor[ind2]=outCS(2);
447 : outstor[ind3]=outCS(3);
448 : }
449 : }
450 : */
451 :
452 : }
453 : } //OMP
454 0 : out.putStorage(outstor, delout);
455 0 : in.freeStorage(instor, delin);
456 0 : polmap.freeStorage(polmap1, delpolmap);
457 0 : }
458 :
459 :
460 :
461 :
462 : void
463 0 : PBMath2DImage::applyJonesFastX(const Float*& reJones,
464 : const Float*& imJones,
465 : const Complex*& instor,
466 : Complex*& outstor,
467 : const Int*& polmap,
468 : Bool /*inverse*/,
469 : Bool /*conjugate*/,
470 : Int ipower, // ie, 1=VP, 2=PB
471 : Float /*cutoff*/,
472 : Bool circular,
473 : Bool forward,
474 : const Int ix, const Int nx, const Int ny, const Int npol){
475 :
476 : ooLong ind0, ind1, ind2, ind3;
477 : // Loop through x, y coordinates of this cube
478 0 : Matrix<Complex> cmat(2,2);
479 :
480 0 : for (ooLong iy=0; iy< ooLong(ny) ; iy++) {
481 :
482 : //sp0(1)=iy;
483 : //sp1(1)=iy;
484 : //sp2(1)=iy;
485 : //sp3(1)=iy;
486 :
487 0 : ind0=ooLong(ix)+ooLong(nx)*iy +(ooLong(nx)*ooLong(ny))*ooLong(polmap[3]);
488 0 : ind1=ooLong(ix)+ooLong(nx)*iy +(ooLong(nx)*ooLong(ny))*ooLong(polmap[2]);
489 0 : ind2=ooLong(ix)+ooLong(nx)*iy +(ooLong(nx)*ooLong(ny))*ooLong(polmap[1]);
490 0 : ind3=ooLong(ix)+ooLong(nx)*iy +(ooLong(nx)*ooLong(ny))*ooLong(polmap[0]);
491 : // E Jones for this pixel
492 0 : mjJones4 j4;
493 0 : if(imJones) {
494 0 : cmat(0,0)=Complex(reJones[ind0], imJones[ind0]);
495 0 : cmat(1,0)=Complex(reJones[ind1], imJones[ind1]);
496 0 : cmat(0,1)=Complex(reJones[ind2], imJones[ind2]);
497 0 : cmat(1,1)=Complex(reJones[ind3], imJones[ind3]);
498 : }
499 : else {
500 0 : cmat(0,0)=Complex(reJones[ind0], 0.0);
501 0 : cmat(1,0)=Complex(reJones[ind1], 0.0);
502 0 : cmat(0,1)=Complex(reJones[ind2], 0.0);
503 0 : cmat(1,1)=Complex(reJones[ind3], 0.0);
504 : }
505 0 : mjJones2 j2(cmat);
506 :
507 : // Make the relevant Jones matrix
508 : // E
509 0 : if(ipower==1) { // VP
510 0 : mjJones2 j2unit(Complex(1.0, 0.0));
511 0 : directProduct(j4, j2, j2unit);
512 : }
513 : // Primary beam = E . conj(E)
514 0 : else if(ipower==2) { // PB
515 : // Make the conjugate before constructing the
516 : // mjJones since otherwise reference semantics
517 : // get us
518 0 : mjJones2 j2conj(conj(cmat));
519 0 : directProduct(j4, j2, j2conj);
520 : }
521 :
522 : // Subtlety - we have to distinguish between applying
523 : // the PB and applying the adjoint. The former is needed
524 : // for predictions (sky->UV) and the latter is needed
525 : // for inversion (UV->sky). For circular polarization
526 : // this affects only the cross hand terms
527 0 : if(!forward) {
528 0 : j4=mjJones4(adjoint(j4.matrix()));
529 : }
530 :
531 : // Now apply the Jones matrix
532 0 : if(npol==1) {
533 : //IPosition ip0(4, ix, iy, 0, 0);
534 0 : ind0=ix+iy*nx;
535 0 : CStokesVector outCS(instor[ind0], 0.0, 0.0, 0.0);
536 0 : outCS*=j4;
537 0 : outstor[ind0]=outCS(0);
538 : }
539 0 : else if(npol==2) {
540 : //IPosition ip0(4, ix, iy, 0, 0);
541 : //IPosition ip1(4, ix, iy, 1, 0);
542 0 : ind0=ix+iy*nx;
543 0 : ind1=ix+iy*nx+nx*ny;
544 0 : if(circular) {
545 0 : CStokesVector outCS(instor[ind0], 0.0, 0.0, instor[ind1]);
546 0 : outCS*=j4;
547 0 : outstor[ind0]=outCS(0);
548 0 : outstor[ind1]=outCS(3);
549 : }
550 : else {
551 0 : CStokesVector outCS(instor[ind0], instor[ind1], 0.0, 0.0);
552 0 : outCS*=j4;
553 0 : outstor[ind0]=outCS(0);
554 0 : outstor[ind1]=outCS(1);
555 : }
556 : }
557 0 : else if(npol==4) {
558 : //IPosition ip0(4, ix, iy, 0, 0);
559 : //IPosition ip1(4, ix, iy, 1, 0);
560 : //IPosition ip2(4, ix, iy, 2, 0);
561 : //IPosition ip3(4, ix, iy, 3, 0);
562 0 : ind0=ix+iy*nx;
563 0 : ind1=ix+iy*nx+nx*ny;
564 0 : ind2=ix+iy*nx+2*nx*ny;
565 0 : ind3=ix+iy*nx+3*nx*ny;
566 0 : CStokesVector outCS(instor[ind0], instor[ind1], instor[ind2], instor[ind3]);
567 0 : outCS*=j4;
568 0 : outstor[ind0]=outCS(0);
569 0 : outstor[ind1]=outCS(1);
570 0 : outstor[ind2]=outCS(2);
571 0 : outstor[ind3]=outCS(3);
572 : }
573 : }
574 :
575 :
576 :
577 :
578 :
579 0 : }
580 :
581 :
582 : // Float to Float on real part of complex stokes
583 : // This is only really useful for the weights image
584 : // from SkyEquation - perhaps make SkyEquation a
585 : // friend?
586 : void
587 0 : PBMath2DImage::applyJones(const Array<Float>* reJones,
588 : const Array<Float>* imJones,
589 : const Array<Float>& in,
590 : Array<Float>& out,
591 : Vector<Int>& polmap,
592 : Float /*cutoff*/,
593 : Bool circular)
594 : {
595 0 : LogIO os(LogOrigin("PBMath2DImage", "applyJones"));
596 :
597 :
598 0 : Matrix<Complex> s(4,4), sAdjoint(4,4);
599 0 : Matrix<Complex> Lambda(4,4);
600 0 : s.set(Complex(0.0,0.0));
601 0 : sAdjoint.set(Complex(0.0));
602 0 : Lambda.set(Complex(0.0,0.0));
603 0 : s=0;
604 0 : if(!circular) {
605 0 : s(0,0)=Complex(0.5);
606 0 : s(0,1)=Complex(0.5);
607 0 : s(1,2)=Complex(0.5);
608 0 : s(1,3)=Complex(0.0,0.5);
609 0 : s(2,2)=Complex(0.5);
610 0 : s(2,3)=Complex(0.0,-0.5);
611 0 : s(3,0)=Complex(0.5);
612 0 : s(3,1)=Complex(-0.5);
613 : }
614 : else {
615 0 : s(0,0)=Complex(0.5);
616 0 : s(0,3)=Complex(0.5);
617 0 : s(1,1)=Complex(0.5);
618 0 : s(1,2)=Complex(0.0,0.5);
619 0 : s(2,1)=Complex(0.5);
620 0 : s(2,2)=Complex(0.0,-0.5);
621 0 : s(3,0)=Complex(0.5);
622 0 : s(3,3)=Complex(-0.5);
623 : }
624 0 : sAdjoint=adjoint(s);
625 :
626 0 : Int nx=in.shape()(0);
627 0 : Int ny=in.shape()(1);
628 0 : Int npol=in.shape()(2);
629 :
630 0 : Matrix<Complex> cmat(2,2);
631 0 : Matrix<Complex> prod(4,4);
632 0 : prod.set(Complex(0.0));
633 0 : IPosition sp0(4, 0, 0, polmap(3), 0);
634 0 : IPosition sp1(4, 0, 0, polmap(2), 0);
635 0 : IPosition sp2(4, 0, 0, polmap(1), 0);
636 0 : IPosition sp3(4, 0, 0, polmap(0), 0);
637 :
638 0 : IPosition ip0(4, 0, 0, 0, 0);
639 0 : IPosition ip1(4, 0, 0, 1, 0);
640 0 : IPosition ip2(4, 0, 0, 2, 0);
641 0 : IPosition ip3(4, 0, 0, 3, 0);
642 0 : for (Int ix=0; ix<nx; ++ix) {
643 0 : sp0(0)=ix;
644 0 : sp1(0)=ix;
645 0 : sp2(0)=ix;
646 0 : sp3(0)=ix;
647 0 : ip0(0)=ix;
648 0 : ip1(0)=ix;
649 0 : ip2(0)=ix;
650 0 : ip3(0)=ix;
651 0 : for (Int iy=0; iy<ny; ++iy) {
652 :
653 0 : sp0(1)=iy;
654 0 : sp1(1)=iy;
655 0 : sp2(1)=iy;
656 0 : sp3(1)=iy;
657 :
658 0 : ip0(1)=iy;
659 0 : ip1(1)=iy;
660 0 : ip2(1)=iy;
661 0 : ip3(1)=iy;
662 :
663 : // Lambda is the covariance matrix of the noise in the
664 : // image plane for each complex stokes e.g. RR, LR, etc.
665 0 : Lambda=0.0;
666 0 : if(npol==1) {
667 0 : Lambda(0,0)=in(ip0);
668 0 : Lambda(3,3)=in(ip0);
669 : }
670 0 : else if(npol==2) {
671 0 : Lambda(0,0)=in(ip0);
672 0 : Lambda(3,3)=in(ip1);
673 : }
674 0 : else if(npol==4) {
675 0 : Lambda(0,0)=in(ip0);
676 0 : Lambda(1,1)=in(ip1);
677 0 : Lambda(2,2)=in(ip2);
678 0 : Lambda(3,3)=in(ip3);
679 : }
680 :
681 : // E Jones for this pixel
682 0 : mjJones4 j4(Matrix<Complex>(4, 4, Complex(0.0, 0.0)));
683 0 : if(imJones) {
684 0 : cmat(0,0)=Complex((*reJones)(sp0), (*imJones)(sp0));
685 0 : cmat(1,0)=Complex((*reJones)(sp1), (*imJones)(sp1));
686 0 : cmat(0,1)=Complex((*reJones)(sp2), (*imJones)(sp2));
687 0 : cmat(1,1)=Complex((*reJones)(sp3), (*imJones)(sp3));
688 : }
689 : else {
690 0 : cmat(0,0)=Complex((*reJones)(sp0), 0.0);
691 0 : cmat(1,0)=Complex((*reJones)(sp1), 0.0);
692 0 : cmat(0,1)=Complex((*reJones)(sp2), 0.0);
693 0 : cmat(1,1)=Complex((*reJones)(sp3), 0.0);
694 : }
695 0 : mjJones2 j2(cmat);
696 :
697 : // Direct product of E with conj(E) ~ PB
698 0 : mjJones2 j2conj(conj(cmat));
699 0 : directProduct(j4, j2, j2conj);
700 : // We need the element by element multiplication
701 0 : Matrix<Complex> matj4(j4.matrix());
702 0 : Matrix<Complex> amatj4(adjoint(matj4));
703 : // Noise covar in complex stokes = Adjoint(PB) . Lambda . PB
704 0 : Lambda=product(amatj4, product(Lambda, matj4));
705 0 : prod=product(sAdjoint,product(Lambda,s));
706 : // Noise covar in real stokes
707 0 : if(npol==1) {
708 0 : out(ip0)=real(prod(0,0))+real(prod(3,3));
709 : }
710 0 : else if (npol==2) {
711 0 : out(ip0)=2.0*real(prod(0,0));
712 0 : out(ip1)=2.0*real(prod(3,3));
713 : }
714 0 : else if (npol==4) {
715 0 : out(ip0)=2.0*real(prod(0,0));
716 0 : out(ip1)=2.0*real(prod(1,1));
717 0 : out(ip2)=2.0*real(prod(2,2));
718 0 : out(ip3)=2.0*real(prod(3,3));
719 : }
720 : }
721 : }
722 :
723 0 : }
724 :
725 0 : ImageInterface<Complex>& PBMath2DImage::apply(const ImageInterface<Complex>& in,
726 : ImageInterface<Complex>& out,
727 : const MDirection& sp,
728 : const Quantity parAngle,
729 : const BeamSquint::SquintType /*doSquint*/,
730 : Bool inverse,
731 : Bool conjugate,
732 : Int ipower, // ie, 1=VP, 2=PB
733 : Float cutoff,
734 : Bool forward)
735 : {
736 0 : LogIO os(LogOrigin("PBMath2DImage", "apply"));
737 :
738 0 : if(!in.shape().isEqual(out.shape(), true)) {
739 : os << "Input and output images have different shapes"
740 0 : << LogIO::EXCEPTION;
741 :
742 : }
743 :
744 : // {
745 : // PagedImage<Float> inImage(in.shape(), in.coordinates(), "RealBeforeApplyJones");
746 : // LatticeExpr<Float> le(real(in));
747 : // inImage.copyData(le);
748 : // }
749 : // {
750 : // PagedImage<Float> inImage(in.shape(), in.coordinates(), "ImagBeforeApplyJones");
751 : // LatticeExpr<Float> le(imag(in));
752 : // inImage.copyData(le);
753 : // }
754 :
755 :
756 0 : const IPosition oshape(out.shape());
757 0 : updateJones(out.coordinates(), oshape, sp, parAngle);
758 :
759 : // Read through the images in pol, x, y, and frequency
760 0 : IPosition iShape(in.shape());
761 0 : IPosition oShape(out.shape());
762 0 : IPosition jShape(reRegridJonesImage_p->shape());
763 0 : Int nx=iShape(0);
764 0 : Int ny=iShape(1);
765 0 : Int npol=iShape(2);
766 :
767 : // Find out if these are circular images
768 0 : Vector<Int> polmap(4);
769 : StokesImageUtil::PolRep polframe;
770 0 : Int insm=StokesImageUtil::CStokesPolMap(polmap, polframe, in.coordinates());
771 0 : Int outsm=StokesImageUtil::CStokesPolMap(polmap, polframe, out.coordinates());
772 0 : if(insm!=outsm) {
773 : os << "Input and Output images have different polarization frames"
774 0 : << LogIO::EXCEPTION;
775 : }
776 0 : Bool circular=(polframe==StokesImageUtil::CIRCULAR);
777 :
778 : // Now get the polarization remapping for the Jones image
779 : //Int jsm= // CStokesPolMap sets polmap and polframe.
780 0 : StokesImageUtil::CStokesPolMap(polmap, polframe,
781 0 : reJonesImage_p->coordinates());
782 :
783 : // For the input and output images, get all polarizations for x, y plane
784 0 : IPosition inPlane(4, nx, ny, npol, 1);
785 0 : LatticeStepper inls(iShape, inPlane, IPosition(4, 0, 1, 2, 3));
786 0 : RO_LatticeIterator<Complex> inli(in, inls);
787 :
788 0 : IPosition outPlane(4, nx, ny, npol, 1);
789 0 : LatticeStepper outls(oShape, outPlane, IPosition(4, 0, 1, 2, 3));
790 0 : LatticeIterator<Complex> outli(out, outls);
791 :
792 : // For the Jones image, get all 4 polarizations for one x, y plane
793 0 : IPosition jPlane(4, nx, ny, 4, 1);
794 0 : LatticeStepper jls(jShape, jPlane, IPosition(4, 0, 1, 2, 3));
795 0 : RO_LatticeIterator<Float> reJonesli(*reRegridJonesImage_p, jls);
796 0 : reJonesli.reset();
797 :
798 : Bool delreal, delimag;
799 0 : if(imRegridJonesImage_p) {
800 0 : RO_LatticeIterator<Float> imJonesli(*imRegridJonesImage_p, jls);
801 :
802 0 : imJonesli.reset();
803 0 : for (inli.reset(), outli.reset(); !inli.atEnd(); inli++, outli++) {
804 0 : const Float *restor=reJonesli.cursor().getStorage(delreal);
805 0 : const Float *imstor=imJonesli.cursor().getStorage(delimag);
806 0 : applyJonesFast(restor, imstor,
807 : inli.cursor(), outli.rwCursor(), polmap,
808 : inverse, conjugate,
809 : ipower, cutoff, circular, forward);
810 0 : reJonesli.cursor().freeStorage(restor, delreal);
811 0 : imJonesli.cursor().freeStorage(imstor, delimag);
812 : }
813 : }
814 : else {
815 0 : for (inli.reset(), outli.reset(); !inli.atEnd(); inli++, outli++) {
816 0 : const Float *restor=reJonesli.cursor().getStorage(delreal);
817 0 : const Float *imstor=NULL;
818 0 : applyJonesFast(restor,
819 : imstor,
820 : inli.cursor(),
821 : outli.rwCursor(), polmap, inverse, conjugate,
822 : ipower, cutoff, circular, forward);
823 0 : reJonesli.cursor().freeStorage(restor, delreal);
824 : }
825 : }
826 :
827 : // {
828 : // PagedImage<Float> outImage(out.shape(), out.coordinates(), "RealAfterApplyJones");
829 : // LatticeExpr<Float> le(real(out));
830 : // outImage.copyData(le);
831 : // }
832 : // {
833 : // PagedImage<Float> outImage(out.shape(), out.coordinates(), "ImagAfterApplyJones");
834 : // LatticeExpr<Float> le(imag(out));
835 : // outImage.copyData(le);
836 : // }
837 0 : return out;
838 : }
839 :
840 0 : ImageInterface<Float>& PBMath2DImage::apply(const ImageInterface<Float>& in,
841 : ImageInterface<Float>& out,
842 : const MDirection& sp,
843 : const Quantity parAngle,
844 : const BeamSquint::SquintType /*doSquint*/,
845 : Float cutoff,
846 : Int /*ipower*/)
847 : {
848 0 : LogIO os(LogOrigin("PBMath2DImage", "apply"));
849 :
850 0 : const IPosition oshape(out.shape());
851 0 : if(!in.shape().isEqual(out.shape(), true)) {
852 : os << "Input and output images have different shapes"
853 0 : << LogIO::EXCEPTION;
854 :
855 : }
856 :
857 0 : updateJones(out.coordinates(), oshape, sp, parAngle);
858 :
859 : // Read through the images in pol, x, y, and frequency
860 0 : IPosition iShape(in.shape());
861 0 : IPosition oShape(out.shape());
862 0 : IPosition jShape(reRegridJonesImage_p->shape());
863 0 : Int nx=iShape(0);
864 0 : Int ny=iShape(1);
865 0 : Int npol=iShape(2);
866 :
867 : // Find out if these are circular images
868 0 : Vector<Int> polmap(4);
869 : StokesImageUtil::PolRep polframe;
870 0 : Int insm=StokesImageUtil::CStokesPolMap(polmap, polframe, in.coordinates());
871 0 : Int outsm=StokesImageUtil::CStokesPolMap(polmap, polframe, out.coordinates());
872 0 : if(insm!=outsm) {
873 : os << "Input and Output images have different polarization frames"
874 0 : << LogIO::EXCEPTION;
875 : }
876 0 : Bool circular=(insm<1);
877 :
878 : // Now get the polarization remapping for the Jones image
879 : //Int jsm=
880 0 : StokesImageUtil::CStokesPolMap(polmap, polframe,
881 0 : reJonesImage_p->coordinates());
882 : // For the input and output images, get all polarizations for x, y plane
883 0 : IPosition inPlane(4, nx, ny, npol, 1);
884 0 : LatticeStepper inls(iShape, inPlane, IPosition(4, 0, 1, 2, 3));
885 0 : RO_LatticeIterator<Float> inli(in, inls);
886 :
887 0 : IPosition outPlane(4, nx, ny, npol, 1);
888 0 : LatticeStepper outls(oShape, outPlane, IPosition(4, 0, 1, 2, 3));
889 0 : LatticeIterator<Float> outli(out, outls);
890 :
891 : // For the Jones image, get all 4 polarizations for one x, y plane
892 0 : IPosition jPlane(4, nx, ny, 4, 1);
893 0 : LatticeStepper jls(jShape, jPlane, IPosition(4, 0, 1, 2, 3));
894 0 : RO_LatticeIterator<Float> reJonesli(*reRegridJonesImage_p, jls);
895 0 : reJonesli.reset();
896 :
897 0 : if(imRegridJonesImage_p) {
898 0 : RO_LatticeIterator<Float> imJonesli(*imRegridJonesImage_p, jls);
899 :
900 0 : imJonesli.reset();
901 0 : for (inli.reset(), outli.reset(); !inli.atEnd(); inli++, outli++) {
902 0 : applyJones(&reJonesli.cursor(), &imJonesli.cursor(),
903 : inli.cursor(), outli.rwCursor(), polmap, cutoff, circular);
904 : }
905 : }
906 : else {
907 0 : for (inli.reset(), outli.reset(); !inli.atEnd(); inli++, outli++) {
908 0 : applyJones(&(reJonesli.cursor()),
909 : (const Array<Float>*)0,
910 : inli.cursor(),
911 : outli.rwCursor(), polmap, cutoff, circular);
912 : }
913 : }
914 :
915 0 : return out;
916 : }
917 :
918 :
919 0 : SkyComponent& PBMath2DImage::apply(SkyComponent& in,
920 : SkyComponent& out,
921 : const MDirection& sp,
922 : const Quantity frequency,
923 : const Quantity parAngle,
924 : const BeamSquint::SquintType /*doSquint*/,
925 : Bool /*inverse*/,
926 : Bool /*conjugate*/,
927 : Int ipower, // ie, 1=VP, 2=PB, 4=PB^2
928 : Float /*cutoff*/,
929 : Bool forward)
930 : {
931 0 : LogIO os(LogOrigin("PBMath2DImage", "apply"));
932 :
933 : // Now get the polarization remapping for the Jones image
934 0 : Vector<Int> polmap(4);
935 : StokesImageUtil::PolRep polframe;
936 :
937 : // jsm (= circular or linear) is not used, but CStokesPolMap also sets polframe.
938 : //Int jsm=
939 0 : StokesImageUtil::CStokesPolMap(polmap, polframe, reJonesImage_p->coordinates());
940 :
941 : // First get the frequency of the output image
942 0 : Double desiredFrequency=frequency.getValue("Hz");
943 :
944 0 : CoordinateSystem reCoords(reJonesImage_p->coordinates());
945 0 : Int spectralIndex=reCoords.findCoordinate(Coordinate::SPECTRAL);
946 0 : AlwaysAssert(spectralIndex>=0, AipsError);
947 : SpectralCoordinate
948 0 : imageSpectralCoord=reCoords.spectralCoordinate(spectralIndex);
949 0 : Double reFrequency=imageSpectralCoord.referenceValue()(0);
950 0 : Double freqScale = desiredFrequency/reFrequency;
951 :
952 : // Find the coordinates of the Sky Component in the
953 : // real and imag Jones images, remembering about the
954 : // pointing position sp. We do the opposite from the
955 : // image case - we convert the component position
956 : // back to an offset from the pointing position
957 : // and then calculate the location in Az,El
958 0 : MDirection compLoc=in.shape().refDirection();
959 0 : MVDirection spmvd(sp.getAngle());
960 0 : MVDirection compmvd(compLoc.getAngle());
961 0 : Vector<Double> world(2);
962 : // Scale for the observing frequency: when the desired frequency
963 : // is higher than the frequency for the Jones image, the separation
964 : // should be LARGER.
965 0 : Double sep=spmvd.separation(compmvd,"rad").getValue();
966 0 : sep*=freqScale;
967 :
968 0 : Double pa=+spmvd.positionAngle(compmvd,"rad").getValue();
969 0 : pa-=parAngle.getValue("rad");
970 0 : world(0)=-sep*sin(pa);
971 0 : world(1)=sep*cos(pa);
972 :
973 0 : Vector<Double> rePix(2), imPix(2);
974 : {
975 0 : Int directionIndex=reCoords.findCoordinate(Coordinate::DIRECTION);
976 : DirectionCoordinate
977 0 : reDirectionCoord=reCoords.directionCoordinate(directionIndex);
978 0 : reDirectionCoord.toPixel(rePix, world);
979 : }
980 0 : if(imJonesImage_p) {
981 0 : CoordinateSystem imCoords(imJonesImage_p->coordinates());
982 0 : Int directionIndex=imCoords.findCoordinate(Coordinate::DIRECTION);
983 : DirectionCoordinate
984 0 : imDirectionCoord=imCoords.directionCoordinate(directionIndex);
985 0 : imDirectionCoord.toPixel(imPix, world);
986 : }
987 :
988 : // E Jones for this pixel
989 0 : mjJones4 j4;
990 :
991 0 : Bool offImage=false;
992 0 : Matrix<Complex> cmat(2,2);
993 0 : Array<Float> reJones;
994 0 : Int nx=reJonesImage_p->shape()(0);
995 0 : Int ny=reJonesImage_p->shape()(1);
996 0 : Int npol=reJonesImage_p->shape()(2);
997 :
998 0 : Int px=Int(rePix(0)+0.5);
999 0 : Int py=Int(rePix(1)+0.5);
1000 0 : if(imJonesImage_p) {
1001 0 : Array<Float> imJones;
1002 0 : if((px>-1)&&(px<nx)&&(py>-1)&&(py<ny)) {
1003 0 : IPosition blc(4, px, py, 0, 0);
1004 0 : IPosition len(4, 1, 1, npol, 1);
1005 0 : IPosition ip0(4, 0, 0, polmap(3), 0);
1006 0 : IPosition ip1(4, 0, 0, polmap(2), 0);
1007 0 : IPosition ip2(4, 0, 0, polmap(1), 0);
1008 0 : IPosition ip3(4, 0, 0, polmap(0), 0);
1009 0 : reJonesImage_p->getSlice(reJones, blc, len);
1010 0 : imJonesImage_p->getSlice(imJones, blc, len);
1011 0 : cmat(0,0)=Complex(reJones(ip0), imJones(ip0));
1012 0 : cmat(1,0)=Complex(reJones(ip1), imJones(ip1));
1013 0 : cmat(0,1)=Complex(reJones(ip2), imJones(ip2));
1014 0 : cmat(1,1)=Complex(reJones(ip3), imJones(ip3));
1015 : }
1016 : else {
1017 : // cerr << freqScale << " " << px << " " << py << " " << world << " " << sep << " " << pa << endl;
1018 0 : offImage = true;
1019 : }
1020 : }
1021 : else {
1022 0 : if((px>-1)&&(px<nx)&&(py>-1)&&(py<ny)) {
1023 0 : IPosition blc(4, px, py, 0, 0);
1024 0 : IPosition len(4, 1, 1, npol, 1);
1025 0 : IPosition ip0(4, 0, 0, polmap(3), 0);
1026 0 : IPosition ip1(4, 0, 0, polmap(2), 0);
1027 0 : IPosition ip2(4, 0, 0, polmap(1), 0);
1028 0 : IPosition ip3(4, 0, 0, polmap(0), 0);
1029 0 : reJonesImage_p->getSlice(reJones, blc, len);
1030 0 : cmat(0,0)=Complex(reJones(ip0), 0.0);
1031 0 : cmat(1,0)=Complex(reJones(ip1), 0.0);
1032 0 : cmat(0,1)=Complex(reJones(ip2), 0.0);
1033 0 : cmat(1,1)=Complex(reJones(ip3), 0.0);
1034 : }
1035 : else {
1036 0 : offImage = true;
1037 : }
1038 :
1039 : }
1040 :
1041 : // Check for the component being off source
1042 0 : if(offImage) {
1043 0 : Vector<Double> zero(4);
1044 0 : zero=0.0;
1045 0 : out=in.copy();
1046 0 : out.flux().setValue(zero);
1047 : }
1048 : else {
1049 0 : mjJones2 j2(cmat);
1050 :
1051 : // Make the relevant Jones matrix
1052 : // E
1053 0 : if(ipower==1) { // VP
1054 0 : mjJones2 j2unit(Complex(1.0, 0.0));
1055 0 : directProduct(j4, j2, j2unit);
1056 : }
1057 : // Primary beam = E . conj(E)
1058 0 : else if(ipower==2) { // PB
1059 : // Make the conjugate before constructing the
1060 : // mjJones since otherwise reference semantics
1061 : // get us
1062 0 : mjJones2 j2conj(conj(cmat));
1063 0 : directProduct(j4, j2, j2conj);
1064 : }
1065 :
1066 : // Subtlety - we have to distinguish between applying
1067 : // the PB and applying the adjoint. The former is needed
1068 : // for predictions (sky->UV) and the latter is needed
1069 : // for inversion (UV->sky)
1070 0 : if(!forward) {
1071 0 : j4=mjJones4(adjoint(j4.matrix()));
1072 : }
1073 :
1074 : // Now apply the Jones matrix
1075 0 : Vector<Complex> compFluxIn(4);
1076 0 : if(polframe==StokesImageUtil::CIRCULAR) {
1077 0 : in.flux().convertPol(ComponentType::CIRCULAR);
1078 : }
1079 : else {
1080 0 : in.flux().convertPol(ComponentType::LINEAR);
1081 : }
1082 : // Reverse index convention
1083 0 : for (Int pol=0;pol<4;pol++) {
1084 0 : compFluxIn(3-pol)=in.flux().value()(pol);
1085 : }
1086 0 : CStokesVector outCS;
1087 0 : outCS=compFluxIn;
1088 0 : outCS*=j4;
1089 0 : Vector<DComplex> compFlux(4);
1090 : // Reverse index convention
1091 0 : for (Int pol=0;pol<4;pol++) {
1092 0 : compFlux(3-pol)=outCS(pol);
1093 : }
1094 0 : out=in.copy();
1095 0 : out.flux().setValue(compFlux);
1096 : }
1097 0 : return out;
1098 : }
1099 :
1100 0 : void PBMath2DImage::updateJones(const CoordinateSystem& coords,
1101 : const IPosition& shape,
1102 : const MDirection& pc,
1103 : const Quantity& paAngle)
1104 : {
1105 :
1106 0 : LogIO os(LogOrigin("PBMath2DImage", "updateJones"));
1107 :
1108 0 : if (!StokesImageUtil::standardImageCoordinates(coords)) {
1109 : os << "Image to be PB corrected does not have standard axes order"
1110 0 : << LogIO::EXCEPTION;
1111 :
1112 : }
1113 : // First get the frequency of the output image
1114 0 : Int spectralIndex=coords.findCoordinate(Coordinate::SPECTRAL);
1115 0 : AlwaysAssert(spectralIndex>=0, AipsError);
1116 : SpectralCoordinate
1117 0 : imageSpectralCoord=coords.spectralCoordinate(spectralIndex);
1118 0 : Double desiredFrequency=imageSpectralCoord.referenceValue()(0);
1119 :
1120 : // Next get the frequency of the Jones image
1121 0 : spectralIndex=reJonesImage_p->coordinates().findCoordinate(Coordinate::SPECTRAL);
1122 0 : AlwaysAssert(spectralIndex>=0, AipsError);
1123 0 : imageSpectralCoord=reJonesImage_p->coordinates().spectralCoordinate(spectralIndex);
1124 0 : Double reFrequency=imageSpectralCoord.referenceValue()(0);
1125 :
1126 0 : Double reFreqScale = reFrequency/desiredFrequency;
1127 : // os << "Desired frequency = " << desiredFrequency << "Hz" << LogIO::POST;
1128 : // os << "Real PB model frequency = " << reFrequency << "Hz" << LogIO::POST;
1129 : // os << "Scaling real Jones image cell size by " << reFreqScale << LogIO::POST;
1130 :
1131 0 : Double imFreqScale = 1.0;
1132 0 : if(imJonesImage_p) {
1133 0 : spectralIndex=imJonesImage_p->coordinates().findCoordinate(Coordinate::SPECTRAL);
1134 0 : AlwaysAssert(spectralIndex>=0, AipsError);
1135 0 : imageSpectralCoord=imJonesImage_p->coordinates().spectralCoordinate(spectralIndex);
1136 0 : Double imFrequency=imageSpectralCoord.referenceValue()(0);
1137 0 : imFreqScale = imFrequency/desiredFrequency;
1138 : // os << "Imag PB model frequency = " << imFrequency << "Hz" << LogIO::POST;
1139 : // os << "Scaling imag Jones image cell size by " << imFreqScale << LogIO::POST;
1140 : }
1141 :
1142 0 : Float pa=-paAngle.getValue("rad");
1143 0 : Matrix<Double> xform(2,2);
1144 0 : xform(0,0) = cos(pa);
1145 0 : xform(1,1) = cos(pa);
1146 0 : xform(0,1) = -sin(pa);
1147 0 : xform(1,0) = +sin(pa);
1148 :
1149 : // cout << "updateJones: position angle = " << 180.0*pa/C::pi << endl;
1150 :
1151 0 : IPosition desiredShape(shape);
1152 0 : desiredShape(2)=reJonesImage_p->shape()(2);
1153 0 : desiredShape(3)=reJonesImage_p->shape()(3);
1154 0 : Int jonesdirectionindex=reJonesImage_p->coordinates().findCoordinate(Coordinate::DIRECTION);
1155 0 : Int imdirectionindex=coords.findCoordinate(Coordinate::DIRECTION);
1156 : DirectionCoordinate
1157 0 : imageDirectionCoord=coords.directionCoordinate(imdirectionindex);
1158 0 : imageDirectionCoord.setWorldAxisUnits(reJonesImage_p->coordinates().directionCoordinate(jonesdirectionindex).worldAxisUnits());
1159 : // Now set the desired coordinates for the regridded Jones images
1160 : // The desired coordinates should have the same direction axis
1161 : // as the image to input image and the same stokes and frequency
1162 : // axes as the input Jones images.
1163 : {
1164 :
1165 :
1166 0 : AlwaysAssert(imdirectionindex>=0, AipsError);
1167 0 : Int spectralIndex=coords.findCoordinate(Coordinate::SPECTRAL);
1168 0 : AlwaysAssert(spectralIndex>=0, AipsError);
1169 : SpectralCoordinate
1170 0 : imageSpectralCoord=coords.spectralCoordinate(spectralIndex);
1171 0 : Vector<Double> freq(1);
1172 0 : freq(0)=desiredFrequency;
1173 0 : imageSpectralCoord.setReferenceValue(freq);
1174 :
1175 0 : CoordinateSystem desiredCoords(reJonesImage_p->coordinates());
1176 0 : desiredCoords.replaceCoordinate(imageDirectionCoord, jonesdirectionindex);
1177 0 : desiredCoords.replaceCoordinate(imageSpectralCoord, spectralIndex);
1178 :
1179 0 : if(reRegridJonesImage_p && desiredCoords.near(reRegridJonesImage_p->coordinates(), 1.0e-4 ) && (desiredShape == (reRegridJonesImage_p->shape()))) {
1180 0 : return;
1181 : }
1182 :
1183 : // Delete any old images
1184 0 : if(reRegridJonesImage_p) delete reRegridJonesImage_p; reRegridJonesImage_p=0;
1185 0 : if(imRegridJonesImage_p) delete imRegridJonesImage_p; imRegridJonesImage_p=0;
1186 :
1187 :
1188 0 : reRegridJonesImage_p = new TempImage<Float>(desiredShape,
1189 0 : desiredCoords);
1190 0 : reRegridJonesImage_p->set(0.0);
1191 0 : if(imJonesImage_p) {
1192 0 : imRegridJonesImage_p = new TempImage<Float>(desiredShape,
1193 0 : desiredCoords);
1194 0 : imRegridJonesImage_p->set(0.0);
1195 : }
1196 : if(0) {
1197 : os << "Regridded image of the real part of Jones matrix:" << LogIO::POST;
1198 : ImageSummary<Float> is(*reRegridJonesImage_p);
1199 : is.list(os);
1200 : }
1201 :
1202 : }
1203 :
1204 : // Now fake the direction coordinates of the input Jones images to have
1205 : // the same RA, DEC, etc as the pointing center and the xform matrix
1206 : // needed to rotate +pa in angle.
1207 0 : CoordinateSystem originalReJonesCoords(reJonesImage_p->coordinates());
1208 : {
1209 0 : AlwaysAssert(imdirectionindex>=0, AipsError);
1210 : DirectionCoordinate
1211 0 : reJonesDirectionCoord=coords.directionCoordinate(imdirectionindex);
1212 : //////////////set the right units
1213 0 : reJonesDirectionCoord.setWorldAxisUnits(reJonesImage_p->coordinates().directionCoordinate(jonesdirectionindex).worldAxisUnits());
1214 : // Set the reference value to the pointing center
1215 0 : reJonesDirectionCoord.setReferenceValue(pc.getAngle().getValue(reJonesDirectionCoord.worldAxisUnits()(0)));
1216 : // Set the xform for a rotation
1217 0 : reJonesDirectionCoord.setLinearTransform(xform);
1218 : // Set the increments correctly to the original values, scaled
1219 : // as needed
1220 0 : Vector<Double> increments(2);
1221 0 : increments=reFreqScale*(*incrementsReJones_p);
1222 0 : increments(0)=-increments(0);
1223 0 : reJonesDirectionCoord.setIncrement(increments);
1224 0 : reJonesDirectionCoord.setReferencePixel(*referencePixelReJones_p);
1225 :
1226 : // Now reset the coordinates of the Real Jones image
1227 0 : CoordinateSystem reJonesCoords(reJonesImage_p->coordinates());
1228 0 : reJonesCoords.replaceCoordinate(reJonesDirectionCoord, jonesdirectionindex);
1229 0 : reJonesImage_p->setCoordinateInfo(reJonesCoords);
1230 :
1231 : if(0) {
1232 : os << "Munged image of the real part of Jones matrix:" << LogIO::POST;
1233 : ImageSummary<Float> is(*reJonesImage_p);
1234 : is.list(os);
1235 : }
1236 :
1237 : }
1238 :
1239 0 : if(imJonesImage_p) {
1240 0 : AlwaysAssert(imdirectionindex>=0, AipsError);
1241 : DirectionCoordinate
1242 0 : imJonesDirectionCoord=coords.directionCoordinate(imdirectionindex);
1243 : // Set the reference value to the pointing center
1244 0 : imJonesDirectionCoord.setWorldAxisUnits(imJonesImage_p->coordinates().directionCoordinate(jonesdirectionindex).worldAxisUnits());
1245 0 : imJonesDirectionCoord.setReferenceValue(pc.getAngle().getValue(imJonesDirectionCoord.worldAxisUnits()(0)));
1246 : // Set the xform for a rotation
1247 0 : imJonesDirectionCoord.setLinearTransform(xform);
1248 : // Set the increments correctly to the original values
1249 0 : Vector<Double> increments(2);
1250 0 : increments=imFreqScale*(*incrementsImJones_p);
1251 0 : increments(0)=-increments(0);
1252 0 : imJonesDirectionCoord.setIncrement(increments);
1253 0 : imJonesDirectionCoord.setReferencePixel(*referencePixelImJones_p);
1254 : // Now reset the coordinates of the Imag Jones image
1255 0 : CoordinateSystem imJonesCoords(imJonesImage_p->coordinates());
1256 0 : imJonesCoords.replaceCoordinate(imJonesDirectionCoord, jonesdirectionindex);
1257 0 : imJonesImage_p->setCoordinateInfo(imJonesCoords);
1258 : }
1259 :
1260 : // Now do the regridding
1261 0 : IPosition whichAxes(2, 0, 1);
1262 0 : ImageRegrid<Float> ir;
1263 0 : ir.regrid(*reRegridJonesImage_p, Interpolate2D::LINEAR, whichAxes,
1264 0 : *reJonesImage_p);
1265 0 : if(imJonesImage_p) {
1266 0 : ir.regrid(*imRegridJonesImage_p, Interpolate2D::LINEAR, whichAxes,
1267 0 : *imJonesImage_p);
1268 : }
1269 :
1270 : // Check for empty PB
1271 0 : LatticeExprNode maxRePB=max(*reRegridJonesImage_p);
1272 0 : LatticeExprNode maxImPB=max(*reRegridJonesImage_p);
1273 0 : if(maxRePB.getFloat()==0.0) {
1274 0 : throw(AipsError("PBMath2DImage: regridded real Jones image is empty"));
1275 : }
1276 0 : if(maxImPB.getFloat()==0.0) {
1277 0 : throw(AipsError("PBMath2DImage: regridded imag Jones image is empty"));
1278 : }
1279 :
1280 : // For debugging purposes
1281 : // if(0) {
1282 : // PagedImage<Float> reRegridJones(reRegridJonesImage_p->shape(),
1283 : // reRegridJonesImage_p->coordinates(),
1284 : // "RealRegridJones");
1285 : // reRegridJones.copyData(*reRegridJonesImage_p);
1286 :
1287 : // if(imJonesImage_p) {
1288 : // PagedImage<Float> imRegridJones(imRegridJonesImage_p->shape(),
1289 : // imRegridJonesImage_p->coordinates(),
1290 : // "ImagRegridJones");
1291 : // imRegridJones.copyData(*imRegridJonesImage_p);
1292 : // }
1293 : // }
1294 : // Reset the coordinates
1295 0 : reJonesImage_p->setCoordinateInfo(originalReJonesCoords);
1296 0 : if(imJonesImage_p) {
1297 0 : imJonesImage_p->setCoordinateInfo(originalReJonesCoords);
1298 : }
1299 : // Reset the coordinates
1300 0 : reJonesImage_p->setCoordinateInfo(originalReJonesCoords);
1301 0 : if(imJonesImage_p) {
1302 0 : imJonesImage_p->setCoordinateInfo(originalReJonesCoords);
1303 : }
1304 : }
1305 :
1306 0 : void PBMath2DImage::checkJonesCongruent(ImageInterface<Float>& reJones,
1307 : ImageInterface<Float>& imJones)
1308 : {
1309 0 : LogIO os(LogOrigin("PBMath2DImage", "checkJonesCongruent"));
1310 0 : if(!reJones.shape().isEqual(imJones.shape(), true)) {
1311 : os << "Real and imaginary primary beam images have different shapes"
1312 0 : << reJones.shape().asVector()
1313 0 : << imJones.shape().asVector()
1314 0 : << LogIO::EXCEPTION;
1315 :
1316 : }
1317 0 : if(reJones.shape().asVector()(2)!=4) {
1318 : os << "Primary beam images must have 4 polarization values: number is "
1319 0 : << reJones.shape().asVector()(2)
1320 0 : << LogIO::EXCEPTION;
1321 :
1322 : }
1323 0 : if (!StokesImageUtil::standardImageCoordinates(reJones)) {
1324 : os << "Primary beam image (real part) does not have standard axes order"
1325 0 : << LogIO::EXCEPTION;
1326 :
1327 : }
1328 0 : if (!StokesImageUtil::standardImageCoordinates(imJones)) {
1329 : os << "Primary beam image (imaginary part) does not have standard axes order"
1330 0 : << LogIO::EXCEPTION;
1331 :
1332 : }
1333 0 : }
1334 :
1335 0 : void PBMath2DImage::checkImageCongruent(ImageInterface<Float>& image)
1336 : {
1337 0 : LogIO os(LogOrigin("PBMath2DImage", "checkImageCongruent"));
1338 0 : if (!StokesImageUtil::standardImageCoordinates(image)) {
1339 : os << "Image does not have standard axes order"
1340 0 : << LogIO::EXCEPTION;
1341 :
1342 : }
1343 0 : }
1344 :
1345 0 : Int PBMath2DImage::support(const CoordinateSystem& cs){
1346 :
1347 0 : Int spectralIndex=cs.findCoordinate(Coordinate::SPECTRAL);
1348 0 : AlwaysAssert(spectralIndex>=0, AipsError);
1349 : SpectralCoordinate
1350 0 : imageSpectralCoord=cs.spectralCoordinate(spectralIndex);
1351 0 : Double desiredFrequency=imageSpectralCoord.referenceValue()(0);
1352 :
1353 : // Next get the frequency of the Jones image
1354 0 : spectralIndex=reJonesImage_p->coordinates().findCoordinate(Coordinate::SPECTRAL);
1355 0 : AlwaysAssert(spectralIndex>=0, AipsError);
1356 0 : imageSpectralCoord=reJonesImage_p->coordinates().spectralCoordinate(spectralIndex);
1357 0 : Double reFrequency=imageSpectralCoord.referenceValue()(0);
1358 :
1359 0 : Double reFreqScale = reFrequency/desiredFrequency;
1360 :
1361 0 : Double npixels=reJonesImage_p->shape()(0)*reFreqScale;
1362 :
1363 0 : Int dirIndex=cs.findCoordinate(Coordinate::DIRECTION);
1364 : DirectionCoordinate
1365 0 : directionCoord=cs.directionCoordinate(dirIndex);
1366 :
1367 0 : Vector<String> dirunit=directionCoord.worldAxisUnits();
1368 :
1369 0 : npixels=npixels/fabs(directionCoord.increment()(0));
1370 0 : dirIndex=reJonesImage_p->coordinates().findCoordinate(Coordinate::DIRECTION);
1371 :
1372 0 : directionCoord=(reJonesImage_p->coordinates()).directionCoordinate(dirIndex);
1373 0 : directionCoord.setWorldAxisUnits(dirunit);
1374 0 : npixels=npixels*fabs(directionCoord.increment()(0));
1375 :
1376 0 : return Int(floor(npixels));
1377 :
1378 : }
1379 :
1380 : } //#End casa namespace
|
