Line data Source code
1 : //# Jones.h: Definition of Jones
2 : //# Copyright (C) 1996,1997,2000,2001,2002,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: casa-feedback@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 :
28 :
29 : #include <synthesis/MeasurementComponents/Jones.h>
30 : #include <casacore/casa/aips.h>
31 : #include <casacore/casa/BasicSL/Complex.h>
32 : #include <iostream>
33 : #include <casacore/casa/Exceptions/Error.h>
34 : #include <casacore/casa/namespace.h>
35 :
36 : using namespace casacore;
37 : namespace casa { //# NAMESPACE CASA - BEGIN
38 :
39 :
40 : // Constructor
41 1566 : Jones::Jones() :
42 1566 : j0_(NULL),
43 1566 : ok0_(NULL),
44 1566 : j_(NULL),
45 1566 : ji_(NULL),
46 1566 : ok_(NULL),
47 1566 : oki_(NULL),
48 1566 : cOne_(1.0),
49 1566 : cZero_(0.0),
50 1566 : scalardata_(false),
51 3132 : vtmp_(VisVector::Four,true)
52 1566 : {}
53 :
54 : // In place invert
55 73080 : void Jones::invert() {
56 :
57 73080 : if (!ok_) throw(AipsError("Illegal use of Jones::invert()"));
58 :
59 73080 : if (ok_[0]&&ok_[1]&&ok_[2]&&ok_[3]) {
60 73080 : Complex det,tmp;
61 73080 : tmp=j_[0];
62 73080 : det=j_[0]*j_[3]-j_[1]*j_[2];
63 73080 : if (abs(det)>0.0) {
64 73080 : j_[0]=j_[3]/det;
65 73080 : j_[1]=-j_[1]/det;
66 73080 : j_[2]=-j_[2]/det;
67 73080 : j_[3]=tmp/det;
68 : } else {
69 0 : zero();
70 0 : ok_[0]=ok_[1]=ok_[2]=ok_[3]=false;
71 : }
72 73080 : }
73 : else {
74 0 : ok_[0]=ok_[1]=ok_[2]=ok_[3]=false;
75 0 : zero();
76 : }
77 :
78 73080 : }
79 :
80 : // Set matrix elements according to ok flag
81 : // (so we don't have to check ok flags atomically in apply)
82 121800 : void Jones::setMatByOk() {
83 :
84 : // Not needed if ok_ not set
85 121800 : if (!ok_) return;
86 :
87 121800 : if (!ok_[0] ||
88 121800 : !ok_[1] ||
89 121800 : !ok_[2] ||
90 121800 : !ok_[3]) {
91 0 : j_[0]=j_[3]=cOne_;
92 0 : j_[1]=j_[2]=cZero_;
93 : }
94 :
95 : }
96 :
97 :
98 :
99 : // In-place multipication with another Jones
100 0 : void Jones::operator*=(const Jones& other) {
101 :
102 0 : switch(other.type()) {
103 : case Jones::General:
104 : case Jones::Diagonal:
105 : case Jones::Scalar:
106 : default:
107 0 : throw(AipsError("General multiplication NYI"));
108 : break;
109 : }
110 : }
111 :
112 : // Apply rightward to a VisVector
113 4384800 : void Jones::applyRight(VisVector& v) const {
114 :
115 : // full general right-ward apply
116 : //
117 : // [w x][A B] = [wA+xC wB+xD] [00+12 01+13]
118 : // [y z][C D] [yA+zC yB+ZD] [20+32 21+23]
119 :
120 :
121 : // flag
122 4384800 : if (v.f_) flagRight(v);
123 :
124 4384800 : switch(v.type()) {
125 4384800 : case VisVector::Four: {
126 4384800 : vtmp_.v_[0] = v.v_[2];
127 4384800 : vtmp_.v_[1] = v.v_[3];
128 4384800 : vtmp_.v_[2] = v.v_[0];
129 4384800 : vtmp_.v_[3] = v.v_[1];
130 :
131 4384800 : v.v_[0]*=(j_[0]);
132 4384800 : v.v_[1]*=(j_[0]);
133 4384800 : v.v_[2]*=(j_[3]);
134 4384800 : v.v_[3]*=(j_[3]);
135 :
136 4384800 : vtmp_.v_[0]*=(j_[1]);
137 4384800 : vtmp_.v_[1]*=(j_[1]);
138 4384800 : vtmp_.v_[2]*=(j_[2]);
139 4384800 : vtmp_.v_[3]*=(j_[2]);
140 :
141 21924000 : for (Int i=0;i<4;++i)
142 17539200 : v.v_[i] += vtmp_.v_[i];
143 :
144 4384800 : break;
145 : }
146 0 : default:
147 0 : throw(AipsError("Jones matrix apply (J::aR) incompatible with VisVector."));
148 : }
149 4384800 : }
150 :
151 : // Apply rightward to a VisVector
152 0 : void Jones::applyRight(VisVector& v, Bool& vflag) const {
153 :
154 0 : if (!ok_) throw(AipsError("Illegal use of Jones::applyRight(v,vflag)"));
155 :
156 0 : applyFlag(vflag);
157 0 : applyRight(v);
158 :
159 0 : }
160 :
161 : // Apply leftward (transposed) to a VisVector
162 4384800 : void Jones::applyLeft(VisVector& v) const {
163 :
164 : // full general left-ward (conjugate transpose) apply
165 : //
166 : // [A B][w y*] = [Aw+Bx* Ay*+Bz] [00+11 02+13]
167 : // [C D][x* z] [Cw+Dx* Cy*+Dz] [20+31 22+33]
168 :
169 :
170 : // flag
171 4384800 : if (v.f_) flagLeft(v);
172 :
173 4384800 : switch(v.type()) {
174 4384800 : case VisVector::Four: {
175 4384800 : vtmp_.v_[0] = v.v_[1];
176 4384800 : vtmp_.v_[1] = v.v_[0];
177 4384800 : vtmp_.v_[2] = v.v_[3];
178 4384800 : vtmp_.v_[3] = v.v_[2];
179 :
180 4384800 : Complex c;
181 4384800 : c=conj(j_[0]);
182 4384800 : v.v_[0]*=c;
183 4384800 : v.v_[2]*=c;
184 4384800 : c=conj(j_[3]);
185 4384800 : v.v_[1]*=c;
186 4384800 : v.v_[3]*=c;
187 :
188 4384800 : c=conj(j_[1]);
189 4384800 : vtmp_.v_[0]*=c;
190 4384800 : vtmp_.v_[2]*=c;
191 4384800 : c=conj(j_[2]);
192 4384800 : vtmp_.v_[1]*=c;
193 4384800 : vtmp_.v_[3]*=c;
194 :
195 21924000 : for (Int i=0;i<4;++i)
196 17539200 : v.v_[i] += vtmp_.v_[i];
197 :
198 4384800 : break;
199 : }
200 0 : default:
201 0 : throw(AipsError("Jones matrix apply (J::aL) incompatible with VisVector."));
202 : }
203 4384800 : }
204 :
205 : // Apply leftward to a VisVector
206 0 : void Jones::applyLeft(VisVector& v, Bool& vflag) const {
207 :
208 0 : if (!ok_) throw(AipsError("Illegal use of Jones::applyLeft(v,vflag)"));
209 :
210 0 : applyFlag(vflag);
211 0 : applyLeft(v);
212 :
213 0 : }
214 :
215 : // Set flags according to solution flags
216 : // (non-corr-dep flag version)
217 0 : void Jones::applyFlag(Bool& vflag) const {
218 0 : vflag|=(!(ok_[0]&&ok_[1]&&ok_[2]&&ok_[3]));
219 0 : }
220 : // Corr-dep flag rightward onto a VisVector
221 4384800 : void Jones::flagRight(VisVector& v) const {
222 4384800 : switch(v.type()) {
223 4384800 : case VisVector::Four: {
224 :
225 4384800 : vtmp_.f_[0]=(v.f_[0]||v.f_[2]||(!ok_[0])||(!ok_[1]));
226 4384800 : vtmp_.f_[1]=(v.f_[1]||v.f_[3]||(!ok_[0])||(!ok_[1]));
227 4384800 : vtmp_.f_[2]=(v.f_[0]||v.f_[2]||(!ok_[2])||(!ok_[3]));
228 4384800 : vtmp_.f_[3]=(v.f_[1]||v.f_[3]||(!ok_[2])||(!ok_[3]));
229 :
230 21924000 : for (Int i=0;i<4;++i)
231 17539200 : v.f_[i] |= vtmp_.f_[i];
232 :
233 4384800 : break;
234 : }
235 0 : default:
236 0 : throw(AipsError("Jones matrix apply (J::fR) incompatible with VisVector."));
237 : }
238 4384800 : }
239 : // Corr-dep flag leftward onto a VisVector
240 4384800 : void Jones::flagLeft(VisVector& v) const {
241 4384800 : switch(v.type()) {
242 4384800 : case VisVector::Four: {
243 :
244 4384800 : vtmp_.f_[0]=(v.f_[0]||v.f_[1]||(!ok_[0])||(!ok_[1]));
245 4384800 : vtmp_.f_[1]=(v.f_[0]||v.f_[1]||(!ok_[2])||(!ok_[3]));
246 4384800 : vtmp_.f_[2]=(v.f_[2]||v.f_[3]||(!ok_[0])||(!ok_[1]));
247 4384800 : vtmp_.f_[3]=(v.f_[2]||v.f_[3]||(!ok_[2])||(!ok_[3]));
248 :
249 21924000 : for (Int i=0;i<4;++i)
250 17539200 : v.f_[i] |= vtmp_.f_[i];
251 :
252 4384800 : break;
253 : }
254 0 : default:
255 0 : throw(AipsError("Jones matrix apply (J::fR) incompatible with VisVector."));
256 : }
257 4384800 : }
258 :
259 :
260 0 : void Jones::zero() {
261 0 : ji_=j_;
262 0 : for (Int i=0;i<4;++i,++ji_)
263 0 : (*ji_)=0.0;
264 0 : }
265 :
266 :
267 : // Constructor
268 108 : JonesGenLin::JonesGenLin() : Jones() {}
269 :
270 : // In place invert
271 410 : void JonesGenLin::invert() {
272 :
273 410 : if (!ok_) throw(AipsError("Illegal use of JonesGenLin::invert()"));
274 :
275 : // In linear approx, we merely negate the off-diag terms!
276 :
277 : /*
278 : if (ok_[0]&&ok_[1]) {
279 : j_[0]*=-1.0;
280 : j_[1]*=-1.0;
281 : } else {
282 : zero();
283 : ok_[0]=ok_[1]=false;
284 : }
285 : */
286 :
287 410 : if (ok_[0])
288 410 : j_[0]*=-1.0;
289 : else
290 0 : j_[0]=0.0;
291 :
292 410 : if (ok_[1])
293 410 : j_[1]*=-1.0;
294 : else
295 0 : j_[1]=0.0;
296 :
297 410 : }
298 :
299 : // Set matrix elements according to ok flag
300 : // (so we don't have to check ok flags atomically in apply)
301 410 : void JonesGenLin::setMatByOk() {
302 :
303 : // Not needed if ok_ not set
304 410 : if (!ok_) return;
305 :
306 410 : if (!ok_[0] ||
307 410 : !ok_[1]) {
308 0 : j_[0]=j_[1]=cOne_;
309 : }
310 :
311 : }
312 :
313 :
314 : // In-place multipication with another Jones
315 0 : void JonesGenLin::operator*=(const Jones& other) {
316 :
317 0 : switch(other.type()) {
318 : case Jones::General:
319 : case Jones::GenLinear:
320 : case Jones::Diagonal:
321 : case Jones::Scalar:
322 : default:
323 0 : throw(AipsError("General multiplication NYI"));
324 : break;
325 : }
326 : }
327 :
328 : // Apply rightward to a VisVector
329 13947840 : void JonesGenLin::applyRight(VisVector& v) const {
330 13947840 : switch(v.type()) {
331 13947840 : case VisVector::Four: {
332 :
333 : // Only adjust cross-hands by pars times parallel hands
334 : //
335 : // [1 X][A B] = [A B+XD]
336 : // [Y 1][C D] [YA+C D ]
337 :
338 : // flag
339 13947840 : if (v.f_) flagRight(v);
340 :
341 13947840 : v.v_[1]+=(j_[0]*v.v_[3]);
342 13947840 : v.v_[2]+=(j_[1]*v.v_[0]);
343 :
344 13947840 : break;
345 : }
346 0 : default:
347 0 : throw(AipsError("JonesGenLin matrix apply (J::aR) incompatible with VisVector."));
348 : }
349 13947840 : }
350 :
351 : // Apply rightward to a VisVector
352 0 : void JonesGenLin::applyRight(VisVector& v, Bool& vflag) const {
353 :
354 0 : if (!ok_) throw(AipsError("Illegal use of JonesGenLin::applyRight(v,vflag)"));
355 :
356 0 : applyFlag(vflag);
357 0 : applyRight(v);
358 :
359 0 : }
360 :
361 : // Apply leftward (conjugate transposed) to a VisVector
362 17455680 : void JonesGenLin::applyLeft(VisVector& v) const {
363 17455680 : switch(v.type()) {
364 17455680 : case VisVector::Four: {
365 :
366 : // Only adjust cross-hands by pars times parallel hands
367 : //
368 : // [A B] [1 Y*] = [A B+AY*]
369 : // [C D] [X* 1 ] [C+DX* D ]
370 :
371 : // flag
372 17455680 : if (v.f_) flagLeft(v);
373 :
374 17455680 : v.v_[1]+=(conj(j_[1])*v.v_[0]);
375 17455680 : v.v_[2]+=(conj(j_[0])*v.v_[3]);
376 :
377 17455680 : break;
378 : }
379 0 : default:
380 0 : throw(AipsError("JonesGenLin matrix apply (J::aL) incompatible with VisVector."));
381 : }
382 17455680 : }
383 :
384 : // Apply leftward to a VisVector
385 0 : void JonesGenLin::applyLeft(VisVector& v, Bool& vflag) const {
386 :
387 0 : if (!ok_) throw(AipsError("Illegal use of JonesGenLin::applyLeft(v,vflag)"));
388 :
389 0 : applyFlag(vflag);
390 0 : applyLeft(v);
391 :
392 0 : }
393 :
394 : // Set flags according to solution flags
395 : // (non-corr-dep flag version)
396 0 : void JonesGenLin::applyFlag(Bool& vflag) const {
397 0 : vflag|=(!(ok_[0]&&ok_[1]));
398 0 : }
399 :
400 : // Corr-dep flag rightward to a VisVector
401 3424320 : void JonesGenLin::flagRight(VisVector& v) const {
402 3424320 : switch(v.type()) {
403 3424320 : case VisVector::Four: {
404 :
405 3424320 : v.f_[1] |= ((!ok_[0])||v.f_[3]);
406 3424320 : v.f_[2] |= ((!ok_[1])||v.f_[0]);
407 :
408 3424320 : break;
409 : }
410 0 : default:
411 0 : throw(AipsError("JonesGenLin matrix apply (JGL::fR) incompatible with VisVector."));
412 : }
413 3424320 : }
414 : // Corr-dep flag leftward to a VisVector
415 6932160 : void JonesGenLin::flagLeft(VisVector& v) const {
416 6932160 : switch(v.type()) {
417 6932160 : case VisVector::Four: {
418 :
419 6932160 : v.f_[1] |= ((!ok_[1])||v.f_[0]);
420 6932160 : v.f_[2] |= ((!ok_[0])||v.f_[3]);
421 :
422 6932160 : break;
423 : }
424 0 : default:
425 0 : throw(AipsError("JonesGenLin matrix apply (JGL::fL) incompatible with VisVector."));
426 : }
427 6932160 : }
428 :
429 0 : void JonesGenLin::zero() {
430 0 : ji_=j_;
431 0 : for (Int i=0;i<2;++i,++ji_)
432 0 : (*ji_)=0.0;
433 0 : }
434 :
435 :
436 : // Constructor
437 1408 : JonesDiag::JonesDiag() : Jones() {}
438 :
439 : // In place invert
440 205556 : void JonesDiag::invert() {
441 :
442 205556 : if (!ok_) throw(AipsError("Illegal use of JonesDiag::invert()"));
443 :
444 205556 : ji_=j_;
445 205556 : oki_=ok_;
446 616668 : for (Int i=0;i<2;++i,++ji_,++oki_)
447 411112 : if ((*oki_) && abs(*ji_)>0.0)
448 408510 : (*ji_)=cOne_/(*ji_);
449 : else {
450 2602 : (*ji_)=cZero_;
451 2602 : (*oki_)=false;
452 : }
453 :
454 205556 : }
455 :
456 : // Set matrix elements according to ok flag
457 : // (so we don't have to check ok flags atomically in apply)
458 5314197 : void JonesDiag::setMatByOk() {
459 :
460 : // Not needed if ok_ not set
461 5314197 : if (!ok_) return;
462 :
463 5314197 : if (!ok_[0]) j_[0]=cOne_;
464 5314197 : if (!ok_[1]) j_[1]=cOne_;
465 :
466 : }
467 :
468 :
469 : // In-place multipication with another Jones
470 0 : void JonesDiag::operator*=(const Jones& other) {
471 :
472 0 : if (!ok_) throw(AipsError("Illegal use of JonesDiag::operator*=()"));
473 :
474 0 : switch(other.type()) {
475 0 : case Jones::General:
476 : case Jones::GenLinear:
477 0 : throw(AipsError("Can't multiply Diagonal by General Jones matrix."));
478 : break;
479 0 : case Jones::Diagonal: {
480 0 : if (ok_[0]&=other.ok_[0]) j_[0]*=other.j_[0];
481 0 : if (ok_[1]&=other.ok_[1]) j_[1]*=other.j_[1];
482 0 : break;
483 : }
484 0 : case Jones::Scalar: {
485 0 : if (ok_[0]&=other.ok_[0]) j_[0]*=other.j_[0];
486 0 : if (ok_[1]&=other.ok_[0]) j_[1]*=other.j_[0];
487 0 : break;
488 : }
489 : }
490 0 : }
491 :
492 : // Apply rightward to a VisVector
493 144727443 : void JonesDiag::applyRight(VisVector& v) const {
494 144727443 : if (v.f_) flagRight(v);
495 144727443 : switch(v.type()) {
496 72448248 : case VisVector::Four: {
497 72448248 : v.v_[0]*=j_[0];
498 72448248 : v.v_[1]*=j_[0];
499 72448248 : v.v_[2]*=j_[1];
500 72448248 : v.v_[3]*=j_[1];
501 72448248 : break;
502 : }
503 72271635 : case VisVector::Two: {
504 72271635 : v.v_[0]*=j_[0];
505 72271635 : v.v_[1]*=j_[1];
506 72271635 : break;
507 : }
508 7560 : case VisVector::One: {
509 7560 : v.v_[0]*=(*j_);
510 7560 : break;
511 : }
512 0 : default:
513 0 : throw(AipsError("Jones matrix apply (JD::aR) incompatible with VisVector."));
514 : }
515 144727443 : }
516 :
517 0 : void JonesDiag::applyRight(VisVector& v, Bool& vflag) const {
518 :
519 0 : if (!ok_) throw(AipsError("Illegal use of JonesDiag::applyRight(v,vflag)"));
520 :
521 0 : applyFlag(vflag);
522 0 : applyRight(v);
523 :
524 0 : }
525 :
526 : // Apply leftward (transposed) to a VisVector
527 154092075 : void JonesDiag::applyLeft(VisVector& v) const {
528 154092075 : if (v.f_) flagLeft(v);
529 154092075 : switch(v.type()) {
530 74165632 : case VisVector::Four: {
531 74165632 : Complex c;
532 74165632 : c=conj(j_[0]);
533 74165632 : v.v_[0]*=c;
534 74165632 : v.v_[2]*=c;
535 74165632 : c=conj(j_[1]);
536 74165632 : v.v_[1]*=c;
537 74165632 : v.v_[3]*=c;
538 74165632 : break;
539 : }
540 79916363 : case VisVector::Two: {
541 79916363 : v.v_[0]*=conj(j_[0]);
542 79916363 : v.v_[1]*=conj(j_[1]);
543 79916363 : break;
544 : }
545 10080 : case VisVector::One: {
546 10080 : v.v_[0]*=conj(*j_);
547 10080 : break;
548 : }
549 0 : default:
550 0 : throw(AipsError("Jones matrix apply (JD::aL) incompatible with VisVector."));
551 : }
552 154092075 : }
553 :
554 0 : void JonesDiag::applyLeft(VisVector& v, Bool& vflag) const {
555 :
556 0 : if (!ok_) throw(AipsError("Illegal use of JonesDiag::applyLeft(v,vflag)"));
557 :
558 0 : applyFlag(vflag);
559 0 : applyLeft(v);
560 0 : }
561 :
562 :
563 : // Set flags according to solution flags
564 : // (non-corr-dep flag version)
565 0 : void JonesDiag::applyFlag(Bool& vflag) const {
566 0 : if (scalardata_)
567 0 : vflag|=(!ok_[0]);
568 : else
569 0 : vflag|=(!(ok_[0]&&ok_[1]));
570 0 : }
571 : // Corr-dep flag right-ward onto a VisVector
572 124526187 : void JonesDiag::flagRight(VisVector& v) const {
573 124526187 : switch(v.type()) {
574 75188736 : case VisVector::Four: {
575 75188736 : v.f_[0]|=(!ok_[0]);
576 75188736 : v.f_[1]|=(!ok_[0]);
577 75188736 : v.f_[2]|=(!ok_[1]);
578 75188736 : v.f_[3]|=(!ok_[1]);
579 75188736 : break;
580 : }
581 49337451 : case VisVector::Two: {
582 49337451 : v.f_[0]|=(!ok_[0]);
583 49337451 : v.f_[1]|=(!ok_[1]);
584 49337451 : break;
585 : }
586 0 : case VisVector::One: {
587 0 : v.f_[0]|=(!ok_[0]);
588 0 : break;
589 : }
590 0 : default:
591 0 : throw(AipsError("Jones matrix apply (JD::aF) incompatible with VisVector."));
592 : }
593 124526187 : }
594 :
595 : // Corr-dep flag left-ward onto a VisVector
596 133890819 : void JonesDiag::flagLeft(VisVector& v) const {
597 133890819 : switch(v.type()) {
598 76906120 : case VisVector::Four: {
599 76906120 : v.f_[0]|=(!ok_[0]);
600 76906120 : v.f_[1]|=(!ok_[1]);
601 76906120 : v.f_[2]|=(!ok_[0]);
602 76906120 : v.f_[3]|=(!ok_[1]);
603 76906120 : break;
604 : }
605 56982179 : case VisVector::Two: {
606 56982179 : v.f_[0]|=(!ok_[0]);
607 56982179 : v.f_[1]|=(!ok_[1]);
608 56982179 : break;
609 : }
610 2520 : case VisVector::One: {
611 2520 : v.f_[0]|=(!ok_[0]);
612 2520 : break;
613 : }
614 0 : default:
615 0 : throw(AipsError("Jones matrix apply (JD::aF) incompatible with VisVector."));
616 : }
617 133890819 : }
618 :
619 :
620 :
621 0 : void JonesDiag::zero() {
622 0 : ji_=j_;
623 0 : for (Int i=0;i<2;++i,++ji_)
624 0 : (*ji_)=0.0;
625 0 : }
626 :
627 :
628 : // Constructor
629 132 : JonesScal::JonesScal() : JonesDiag() {}
630 :
631 : // In place invert
632 509675 : void JonesScal::invert() {
633 :
634 509675 : if (!ok_) throw(AipsError("Illegal use of JonesScal::invert()"));
635 :
636 509675 : if ((*ok_) && abs(*j_)>0.0)
637 502367 : (*j_)=cOne_/(*j_);
638 : else {
639 7308 : (*j_)=cZero_;
640 7308 : (*ok_)=false;
641 : }
642 :
643 509675 : }
644 :
645 :
646 : // Set matrix elements according to ok flag
647 : // (so we don't have to check ok flags atomically in apply)
648 624515 : void JonesScal::setMatByOk() {
649 :
650 : // Not needed if ok_ not set
651 624515 : if (!ok_) return;
652 :
653 624515 : if (!ok_[0]) j_[0]=cOne_;
654 :
655 : }
656 :
657 :
658 : // In-place multipication with another Jones
659 0 : void JonesScal::operator*=(const Jones& other) {
660 :
661 0 : if (!ok_) throw(AipsError("Illegal use of JonesScal::operator*="));
662 :
663 0 : switch(other.type()) {
664 0 : case Jones::General:
665 : case Jones::GenLinear:
666 0 : throw(AipsError("Can't multiply Scalar Jones by General Jones matrix."));
667 : break;
668 0 : case Jones::Diagonal:
669 0 : throw(AipsError("Can't multiply Scalar Jones by Diagonal Jones matrix."));
670 : break;
671 0 : case Jones::Scalar: {
672 0 : if ((*ok_)&=(*other.ok_)) (*j_)*=(*other.j_);
673 0 : break;
674 : }
675 : }
676 0 : }
677 :
678 :
679 : // Apply rightward to a VisVector
680 39779845 : void JonesScal::applyRight(VisVector& v) const {
681 39779845 : if (v.f_) flagRight(v);
682 198871275 : for (Int i=0;i<v.vistype_;++i)
683 159091430 : v.v_[i] *= (*j_);
684 39779845 : }
685 :
686 13975 : void JonesScal::applyRight(VisVector& v, Bool& vflag) const {
687 :
688 13975 : if (!ok_) throw(AipsError("Illegal use of JonesScal::applyRight(v,vflag)"));
689 :
690 13975 : applyFlag(vflag);
691 13975 : applyRight(v);
692 13975 : }
693 :
694 : // Apply leftward (transposed) to a VisVector
695 48387580 : void JonesScal::applyLeft(VisVector& v) const {
696 48387580 : if (v.f_) flagLeft(v);
697 48387580 : Complex c;
698 48387580 : c=conj(*j_);
699 241909950 : for (Int i=0;i<v.vistype_;++i)
700 193522370 : v.v_[i] *= c;
701 48387580 : }
702 :
703 13975 : void JonesScal::applyLeft(VisVector& v, Bool& vflag) const {
704 :
705 13975 : if (!ok_) throw(AipsError("Illegal use of JonesScal::applyLeft(v,vflag)"));
706 :
707 13975 : applyFlag(vflag);
708 13975 : applyLeft(v);
709 :
710 13975 : }
711 :
712 : // Set flags according to solution flags
713 : // (non-corr-dep flag version)
714 27950 : void JonesScal::applyFlag(Bool& vflag) const {
715 27950 : vflag|=(!*ok_);
716 27950 : }
717 : // Corr-dep flag rightward onto a VisVector
718 : // (NB: flagLeft call this since flagging commutes for scalars)
719 53708535 : void JonesScal::flagRight(VisVector& v) const {
720 268542675 : for (Int i=0;i<v.vistype_;++i) v.f_[i]|=(!*ok_);
721 53708535 : }
722 :
723 :
724 0 : void JonesScal::zero() {
725 0 : (*j_)=0.0;
726 0 : }
727 :
728 :
729 : // GLOBALS---------------------------
730 :
731 :
732 1566 : Jones* createJones(const Jones::JonesType& jtype) {
733 1566 : Jones *j(NULL);
734 1566 : switch (jtype) {
735 50 : case Jones::General:
736 50 : j = new Jones();
737 50 : break;
738 108 : case Jones::GenLinear:
739 108 : j = new JonesGenLin();
740 108 : break;
741 1276 : case Jones::Diagonal:
742 1276 : j = new JonesDiag();
743 1276 : break;
744 132 : case Jones::Scalar:
745 132 : j = new JonesScal();
746 132 : break;
747 : }
748 1566 : return j;
749 : }
750 :
751 0 : void apply(const Jones& j1, VisVector& v, const Jones& j2) {
752 :
753 : // Apply Jones matrix from right, then left
754 0 : j1.applyRight(v);
755 0 : j2.applyLeft(v);
756 :
757 0 : }
758 :
759 0 : void apply(const Jones& j1, VisVector& v, const Jones& j2, Bool& vflag) {
760 :
761 : // Apply Jones matrix from right, then left
762 0 : j1.applyRight(v,vflag);
763 0 : j2.applyLeft(v,vflag);
764 :
765 0 : }
766 :
767 : // print out a Jones matrix (according to type)
768 0 : ostream& operator<<(ostream& os, const Jones& mat) {
769 :
770 : Complex *ji;
771 0 : ji=mat.j_;
772 :
773 0 : switch (mat.type()) {
774 0 : case Jones::General:
775 0 : cout << "General Jones: " << endl;
776 0 : for (Int i=0;i<2;++i) {
777 0 : cout << " [" << *ji++;
778 0 : cout << ", " << *ji++ << "]";
779 0 : if (i<1) cout << endl;
780 : }
781 0 : break;
782 0 : case Jones::GenLinear:
783 0 : cout << "GenLinear Jones: (off diag) " << endl;
784 0 : cout << " [" << *ji++;
785 0 : cout << ", " << *ji << "]";
786 0 : break;
787 0 : case Jones::Diagonal:
788 0 : cout << "Diagonal Jones: " << endl;
789 0 : cout << " [" << *ji++;
790 0 : cout << ", " << *ji << "]";
791 0 : break;
792 0 : case Jones::Scalar:
793 0 : cout << "Scalar Jones: " << endl;
794 0 : cout << "[" << *ji << "]";
795 0 : break;
796 : }
797 :
798 0 : return os;
799 : }
800 :
801 : // Return the enum for from an int
802 0 : Jones::JonesType jonesType(const Int& n) {
803 0 : switch (n) {
804 0 : case 1:
805 0 : return Jones::Scalar;
806 0 : case 2:
807 0 : return Jones::Diagonal;
808 0 : case 4:
809 0 : return Jones::General;
810 0 : default:
811 0 : throw(AipsError("Bad JonesType."));
812 :
813 : }
814 :
815 : }
816 :
817 : } //# NAMESPACE CASA - END
|
