Line data Source code
1 : //# Mueller.cc: Implementation of Mueller
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 : #include <synthesis/MeasurementComponents/Mueller.h>
29 : #include <synthesis/MeasurementComponents/Jones.h>
30 : #include <synthesis/MeasurementComponents/VisVector.h>
31 : #include <casacore/casa/aips.h>
32 : #include <casacore/casa/BasicSL/Complex.h>
33 : #include <iostream>
34 : #include <casacore/casa/Exceptions/Error.h>
35 :
36 : #include <casacore/casa/namespace.h>
37 :
38 : using namespace casacore;
39 : namespace casa { //# NAMESPACE CASA - BEGIN
40 :
41 : // Constructor
42 0 : Mueller::Mueller() :
43 0 : m0_(NULL),
44 0 : ok0_(NULL),
45 0 : m_(NULL),
46 0 : mi_(NULL),
47 0 : ok_(NULL),
48 0 : oki_(NULL),
49 0 : cOne_(1.0),
50 0 : cZero_(0.0),
51 0 : scalardata_(false),
52 0 : vtmp_(VisVector::Four,true)
53 0 : {}
54 :
55 : // Formation from Jones matrix outer product: General version
56 0 : void Mueller::fromJones(const Jones& jones1, const Jones& jones2) {
57 :
58 0 : if (!ok_ || !jones1.ok_ || !jones2.ok_)
59 0 : throw(AipsError("Illegal use of Mueller::fromJones."));
60 :
61 0 : mi_=m_;
62 0 : oki_=ok_;
63 : Int i,j,io2,jo2,im2,jm2;
64 0 : for (i=0;i<4;++i) {
65 0 : io2=(i/2)*2;
66 0 : im2=(i%2)*2;
67 0 : for (j=0;j<4;++j,++mi_,++oki_){
68 0 : jo2=(j/2);
69 0 : jm2=(j%2);
70 0 : (*mi_) = jones1.j_[io2+jo2]*conj(jones2.j_[im2+jm2]);
71 0 : (*oki_) = (jones1.ok_[io2+jo2] && jones2.ok_[im2+jm2]);
72 : }
73 : }
74 0 : }
75 :
76 : // In-place multiply onto a VisVector: General version
77 0 : void Mueller::apply(VisVector& v) {
78 :
79 0 : switch (v.type()) {
80 0 : case VisVector::Four: {
81 0 : vtmp_ = v; // remember input (necessarily)
82 : // IS VV COPY COMPLETE?
83 0 : mi_=m_; // iteration within the current matrix
84 : Complex *vo,*vi;
85 0 : Complex tc;
86 0 : vo=v.v_;
87 : Int i,j;
88 0 : for (j=0;j<4;++j,++vo) {
89 0 : vi=vtmp_.v_;
90 0 : (*vo) = (*vi);
91 0 : (*vo) *= (*mi_);
92 0 : ++mi_;
93 0 : ++vi;
94 0 : for (i=1;i<4;++i,++mi_,++vi) { //
95 0 : tc = (*vi);
96 0 : tc *= (*mi_);
97 0 : (*vo) += (tc); // ((*mi_)*(*vi));
98 : }
99 : }
100 0 : break;
101 : }
102 0 : default:
103 0 : throw(AipsError("Mueller/VisVector type mismatch in apply"));
104 : break;
105 : }
106 0 : }
107 :
108 0 : void Mueller::apply(VisVector& /*v*/, Bool& /*vflag*/) {
109 :
110 0 : throw(AipsError("Mueller::apply(v,vflag) (general) NYI."));
111 :
112 : }
113 :
114 0 : void Mueller::applyFlag(Bool& /*vflag*/) {
115 0 : throw(AipsError("Mueller::applyFlag(vflag) (general) NYI."));
116 : }
117 :
118 0 : void Mueller::flag(VisVector& /*v*/) {
119 0 : throw(AipsError("Mueller::flag(v) (general) NYI."));
120 : }
121 :
122 :
123 0 : void Mueller::invert() {
124 :
125 0 : throw(AipsError("Invalid attempt to directly invert a general Mueller."));
126 :
127 : }
128 :
129 : // Set matrix elements according to ok flag
130 : // (so we don't have to check ok flags atomically in apply)
131 0 : void Mueller::setMatByOk() {
132 :
133 0 : throw(AipsError("Invalid attempt to setMatByOk."));
134 :
135 : }
136 :
137 : // Multiply onto a vis VisVector, preserving input (copy then in-place apply)
138 0 : void Mueller::apply(VisVector& out, const VisVector& in) {
139 0 : out=in;
140 0 : apply(out);
141 0 : }
142 :
143 0 : void Mueller::zero() {
144 0 : mi_=m_;
145 0 : for (Int i=0;i<16;++i,++mi_)
146 0 : *mi_=0.0;
147 0 : }
148 :
149 : // Constructor
150 0 : MuellerDiag::MuellerDiag() : Mueller() {}
151 :
152 : // Formation from Jones matrix outer product: optimized Diagonal version
153 0 : void MuellerDiag::fromJones(const Jones& jones1, const Jones& jones2) {
154 :
155 0 : if (!ok_ || !jones1.ok_ || !jones2.ok_)
156 0 : throw(AipsError("Illegal use of MuellerDiag::fromJones."));
157 :
158 0 : mi_=m_;
159 0 : oki_=ok_;
160 0 : for (Int i=0;i<4;++i,++mi_) {
161 0 : (*mi_) = jones1.j_[i/2]*conj(jones2.j_[i%2]);
162 0 : (*oki_) = (jones1.ok_[i/2]&&jones2.ok_[i%2]);
163 : }
164 0 : }
165 :
166 0 : void MuellerDiag::invert() {
167 :
168 0 : if (!ok_) throw(AipsError("Illegal use of MuellerDiag::invert."));
169 :
170 0 : mi_=m_;
171 0 : oki_=ok_;
172 0 : for (Int i=0;i<4;++i,++mi_,++oki_)
173 0 : if ((*oki_) && abs(*mi_)>0.0)
174 0 : (*mi_) = cOne_/(*mi_);
175 : else {
176 0 : (*mi_)=cZero_;
177 0 : (*oki_)=false;
178 : }
179 :
180 0 : }
181 :
182 : // Set matrix elements according to ok flag
183 : // (so we don't have to check ok flags atomically in apply)
184 0 : void MuellerDiag::setMatByOk() {
185 :
186 : // Not needed if ok_ not set
187 0 : if (!ok_) return;
188 :
189 0 : if (!ok_[0]) m_[0]=cOne_;
190 0 : if (!ok_[1]) m_[1]=cOne_;
191 0 : if (!ok_[2]) m_[2]=cOne_;
192 0 : if (!ok_[3]) m_[3]=cOne_;
193 :
194 : }
195 :
196 :
197 :
198 : // In-place multiply onto a VisVector: optimized Diagonal version
199 0 : void MuellerDiag::apply(VisVector& v) {
200 :
201 : // Flag
202 0 : if (v.f_) flag(v);
203 :
204 0 : mi_=m_;
205 0 : Complex *vi(v.v_);
206 :
207 0 : switch (v.type()) {
208 0 : case VisVector::Four: {
209 : // element-by-element apply of Mueller diagonal to VisVector 4-vector
210 0 : for (Int i=0;i<4;++i,++mi_,++vi) (*vi)*=(*mi_);
211 0 : break;
212 : }
213 0 : case VisVector::Two: {
214 : // Mueller corner elements apply to VisVector 2-vector
215 0 : for (Int i=0;i<2;++i,++vi,mi_+=3) (*vi)*=(*mi_);
216 0 : break;
217 : }
218 0 : case VisVector::One: {
219 : // Mueller corner element used as scalar (pol-sensitivity TBD)
220 0 : (*vi)*=(*mi_);
221 0 : break;
222 : }
223 : }
224 0 : }
225 :
226 0 : void MuellerDiag::apply(VisVector& v, Bool& vflag) {
227 :
228 0 : if (!ok_) throw(AipsError("Illegal use of MuellerDiag::apply(v,vflag)."));
229 :
230 0 : applyFlag(vflag);
231 0 : apply(v);
232 :
233 0 : }
234 :
235 0 : void MuellerDiag::applyFlag(Bool& vflag) {
236 :
237 0 : if (!ok_) throw(AipsError("Illegal use of MuellerDiag::applyFlag(vflag)."));
238 :
239 0 : if (scalardata_)
240 0 : vflag|=(!ok_[0]);
241 : else
242 0 : vflag|=(!(ok_[0]&&ok_[1]&&ok_[2]&&ok_[3]));
243 0 : }
244 :
245 : // Flag
246 0 : void MuellerDiag::flag(VisVector& v) {
247 :
248 0 : oki_=ok_;
249 0 : Bool *fi(v.f_);
250 :
251 0 : switch (v.type()) {
252 0 : case VisVector::Four: {
253 : // element-by-element flag of Mueller diagonal to VisVector 4-vector
254 0 : for (Int i=0;i<4;++i,++oki_,++fi) (*fi)|=(!(*oki_));
255 0 : break;
256 : }
257 0 : case VisVector::Two: {
258 : // Mueller corner elements apply to VisVector 2-vector
259 0 : for (Int i=0;i<2;++i,++fi,oki_+=3) (*fi)|=(!(*oki_));
260 0 : break;
261 : }
262 0 : case VisVector::One: {
263 : // Mueller corner element used as scalar (pol-sensitivity TBD)
264 0 : (*fi)|=(!(*oki_));
265 0 : break;
266 : }
267 : }
268 0 : }
269 :
270 :
271 0 : void MuellerDiag::zero() {
272 0 : mi_=m_;
273 0 : for (Int i=0;i<4;++i,++mi_)
274 0 : *mi_=0.0;
275 0 : }
276 :
277 : // Constructor
278 0 : MuellerDiag2::MuellerDiag2() : MuellerDiag() {}
279 :
280 : // Formation from Jones matrix outer product: optimized Diag2 version
281 0 : void MuellerDiag2::fromJones(const Jones& jones1, const Jones& jones2) {
282 :
283 0 : if (!ok_ || !jones1.ok_ || !jones2.ok_)
284 0 : throw(AipsError("Illegal use of MuellerDiag2::fromJones."));
285 :
286 0 : mi_=m_;
287 0 : oki_=ok_;
288 0 : for (Int i=0;i<2;++i,++mi_,++oki_) {
289 0 : (*mi_) = jones1.j_[i]*conj(jones2.j_[i]);
290 0 : (*oki_) = (jones1.ok_[i]&&jones2.ok_[i]);
291 : }
292 0 : }
293 :
294 0 : void MuellerDiag2::invert() {
295 :
296 0 : if (!ok_) throw(AipsError("Illegal use of MuellerDiag2::invert."));
297 :
298 0 : mi_=m_;
299 0 : oki_=ok_;
300 0 : for (Int i=0;i<2;++i,++mi_,++oki_)
301 0 : if ((*oki_)&& abs(*mi_)>0.0)
302 0 : (*mi_) = cOne_/(*mi_);
303 : else {
304 0 : (*mi_) = cZero_;
305 0 : (*oki_) = false;
306 : }
307 :
308 0 : }
309 :
310 : // Set matrix elements according to ok flag
311 : // (so we don't have to check ok flags atomically in apply)
312 0 : void MuellerDiag2::setMatByOk() {
313 :
314 : // Not needed if ok_ not set
315 0 : if (!ok_) return;
316 :
317 0 : if (!ok_[0]) m_[0]=cOne_;
318 0 : if (!ok_[1]) m_[1]=cOne_;
319 :
320 : }
321 :
322 :
323 : // In-place multiply onto a VisVector: optimized Diag2 version
324 0 : void MuellerDiag2::apply(VisVector& v) {
325 :
326 : // Flag
327 0 : if (v.f_) flag(v);
328 :
329 0 : switch (v.type()) {
330 0 : case VisVector::Four:
331 : // Apply of "corners-only" Mueller to VisVector 4-vector (x-hands zeroed)
332 0 : v.v_[0]*=m_[0];
333 0 : v.v_[1]*=0.0;
334 0 : v.v_[2]*=0.0;
335 0 : v.v_[3]*=m_[1];
336 0 : break;
337 0 : case VisVector::Two:
338 : // Element-by-element apply of "corners-only" Mueller to VisVector 2-vector
339 0 : for (Int i=0;i<2;++i) v.v_[i]*=m_[i];
340 0 : break;
341 0 : case VisVector::One: {
342 : // Mueller corner element used as scalar (pol-sensitivity TBD)
343 0 : v.v_[0]*=(*m_);
344 0 : break;
345 : }
346 : }
347 :
348 0 : }
349 :
350 0 : void MuellerDiag2::apply(VisVector& v, Bool& vflag) {
351 :
352 0 : if (!ok_) throw(AipsError("Illegal use of MuellerDiag2::apply(v,vflag)."));
353 :
354 0 : applyFlag(vflag);
355 0 : apply(v);
356 :
357 0 : }
358 :
359 0 : void MuellerDiag2::applyFlag(Bool& vflag) {
360 :
361 0 : if (!ok_) throw(AipsError("Illegal use of MuellerDiag2::applyFlag(vflag)."));
362 :
363 0 : if (scalardata_)
364 0 : vflag|=(!ok_[0]);
365 : else
366 0 : vflag|=(!(ok_[0]&&ok_[1]));
367 0 : }
368 :
369 : // Flag a VisVector: optimized Diag2 version
370 0 : void MuellerDiag2::flag(VisVector& v) {
371 :
372 0 : switch (v.type()) {
373 0 : case VisVector::Four:
374 : // Apply of "corners-only" Mueller to VisVector 4-vector
375 0 : v.f_[0]|=(!ok_[0]);
376 0 : v.f_[3]|=(!ok_[1]);
377 0 : break;
378 0 : case VisVector::Two:
379 : // Element-by-element apply of "corners-only" Mueller to VisVector 2-vector
380 0 : for (Int i=0;i<2;++i) v.f_[i]|=(!ok_[i]);
381 0 : break;
382 0 : case VisVector::One: {
383 : // Mueller corner element used as scalar (pol-sensitivity TBD)
384 0 : v.f_[0]|=(!ok_[0]);
385 0 : break;
386 : }
387 : }
388 :
389 0 : }
390 :
391 :
392 :
393 :
394 0 : void MuellerDiag2::zero() {
395 0 : mi_=m_;
396 0 : for (Int i=0;i<2;++i,++mi_)
397 0 : *mi_=0.0;
398 0 : }
399 :
400 : // Constructor
401 0 : MuellerScal::MuellerScal() : MuellerDiag() {}
402 :
403 : // Formation from Jones matrix outer product: optimized Scalar version
404 0 : void MuellerScal::fromJones(const Jones& jones1, const Jones& jones2) {
405 :
406 0 : if (!ok_ || !jones1.ok_ || !jones2.ok_)
407 0 : throw(AipsError("Illegal use of MuellerScal::fromJones."));
408 :
409 0 : (*m_) = (*jones1.j_)*conj(*jones2.j_);
410 0 : (*ok_) = ((*jones1.ok_)&&(*jones2.ok_));
411 0 : }
412 :
413 0 : void MuellerScal::invert() {
414 :
415 0 : if (!ok_) throw(AipsError("Illegal use of MuellerScal::invert."));
416 :
417 0 : if ((*ok_) && abs(*m_)>0.0)
418 0 : (*m_) = cOne_/(*m_);
419 : else {
420 0 : zero();
421 0 : (*ok_)=false;
422 : }
423 :
424 0 : }
425 :
426 : // Set matrix elements according to ok flag
427 : // (so we don't have to check ok flags atomically in apply)
428 0 : void MuellerScal::setMatByOk() {
429 :
430 : // Not needed if ok_ not set
431 0 : if (!ok_) return;
432 :
433 0 : if (!ok_[0]) m_[0]=cOne_;
434 :
435 : }
436 :
437 : // In-place multiply onto a VisVector: optimized Scalar version
438 0 : void MuellerScal::apply(VisVector& v) {
439 :
440 : // Flag
441 0 : if (v.f_) flag(v);
442 :
443 : // Apply single value to all vector elements
444 0 : for (Int i=0;i<v.vistype_;i++) v.v_[i]*=(*m_);
445 0 : }
446 :
447 0 : void MuellerScal::apply(VisVector& v, Bool& vflag) {
448 :
449 0 : if (!ok_) throw(AipsError("Illegal use of MuellerScal::apply(v,vflag)."));
450 :
451 0 : applyFlag(vflag);
452 0 : apply(v);
453 :
454 0 : }
455 :
456 0 : void MuellerScal::applyFlag(Bool& vflag) {
457 :
458 0 : if (!ok_) throw(AipsError("Illegal use of MuellerScal::applyFlag(vflag)."));
459 :
460 0 : vflag|=(!*ok_);
461 0 : }
462 :
463 : // Flag a VisVector: optimized Scalar version
464 0 : void MuellerScal::flag(VisVector& v) {
465 : // Apply single value to all vector elements
466 0 : Bool f=(!ok_[0]);
467 0 : for (Int i=0;i<v.vistype_;i++) v.f_[i]|=f;
468 0 : }
469 :
470 :
471 0 : void MuellerScal::zero() {
472 0 : *m_=0.0;
473 0 : }
474 :
475 :
476 : // Constructor
477 0 : AddMuellerDiag2::AddMuellerDiag2() : MuellerDiag2() {}
478 :
479 0 : void AddMuellerDiag2::invert() {
480 :
481 : // "invert" means "negate" for additive term
482 :
483 0 : if (!ok_) throw(AipsError("Illegal use of AddMuellerDiag2::invert."));
484 :
485 0 : mi_=m_;
486 0 : oki_=ok_;
487 0 : for (Int i=0;i<2;++i,++mi_,++oki_)
488 0 : if ((*oki_))
489 0 : (*mi_)*=-1.0; // negate
490 : else {
491 0 : (*mi_) = cZero_;
492 0 : (*oki_) = false;
493 : }
494 :
495 0 : }
496 :
497 : // Set matrix elements according to ok flag
498 : // (so we don't have to check ok flags atomically in apply)
499 0 : void AddMuellerDiag2::setMatByOk() {
500 :
501 : // Not needed if ok_ not set
502 0 : if (!ok_) return;
503 :
504 0 : if (!ok_[0]) m_[0]=cZero_;
505 0 : if (!ok_[1]) m_[1]=cZero_;
506 :
507 : }
508 :
509 : // In-place multiply onto a VisVector: optimized Diag2 version
510 0 : void AddMuellerDiag2::apply(VisVector& v) {
511 :
512 : // Flag
513 0 : if (v.f_) flag(v);
514 :
515 : // "apply" means "add" for additive term
516 :
517 0 : switch (v.type()) {
518 0 : case VisVector::Four:
519 : // Apply of "corners-only" Mueller to VisVector 4-vector (x-hands zeroed)
520 0 : v.v_[0]+=m_[0]; // add
521 0 : v.v_[1]*=0.0;
522 0 : v.v_[2]*=0.0;
523 0 : v.v_[3]+=m_[1]; // add
524 0 : break;
525 0 : case VisVector::Two:
526 : // Element-by-element apply of "corners-only" Mueller to VisVector 2-vector
527 0 : for (Int i=0;i<2;++i) v.v_[i]+=m_[i]; // add
528 0 : break;
529 0 : case VisVector::One: {
530 : // Mueller corner element used as scalar (pol-sensitivity TBD)
531 0 : v.v_[0]+=(*m_); // add
532 0 : break;
533 : }
534 : }
535 :
536 0 : }
537 :
538 : // Constructor
539 0 : AddMuellerDiag::AddMuellerDiag() : MuellerDiag() {}
540 :
541 0 : void AddMuellerDiag::invert() {
542 :
543 : // "invert" means "negate" for additive term
544 :
545 0 : if (!ok_) throw(AipsError("Illegal use of AddMuellerDiag2::invert."));
546 :
547 0 : mi_=m_;
548 0 : oki_=ok_;
549 0 : for (Int i=0;i<4;++i,++mi_,++oki_)
550 0 : if ((*oki_))
551 0 : (*mi_)*=-1.0; // negate
552 : else {
553 0 : (*mi_) = cZero_;
554 0 : (*oki_) = false;
555 : }
556 :
557 0 : }
558 :
559 : // Set matrix elements according to ok flag
560 : // (so we don't have to check ok flags atomically in apply)
561 0 : void AddMuellerDiag::setMatByOk() {
562 :
563 : // Not needed if ok_ not set
564 0 : if (!ok_) return;
565 :
566 0 : if (!ok_[0]) m_[0]=cZero_;
567 0 : if (!ok_[1]) m_[1]=cZero_;
568 0 : if (!ok_[2]) m_[2]=cZero_;
569 0 : if (!ok_[3]) m_[3]=cZero_;
570 :
571 : }
572 :
573 :
574 : // In-place add onto a VisVector: optimized Diag2 version
575 0 : void AddMuellerDiag::apply(VisVector& v) {
576 :
577 : // Flag
578 0 : if (v.f_) flag(v);
579 :
580 : // "apply" means "add" for additive term
581 :
582 0 : switch (v.type()) {
583 0 : case VisVector::Four:
584 : // Direct apply of Mueller to VisVector 4-vector
585 0 : v.v_[0]+=m_[0]; // add
586 0 : v.v_[1]+=m_[1];
587 0 : v.v_[2]+=m_[2];
588 0 : v.v_[3]+=m_[3];
589 0 : break;
590 0 : case VisVector::Two:
591 : // Element-by-element apply of Mueller corners to VisVector 2-vector
592 0 : v.v_[0]+=m_[0]; // add
593 0 : v.v_[1]+=m_[3]; // add
594 0 : break;
595 0 : case VisVector::One: {
596 : // Mueller corner element used as scalar (pol-sensitivity TBD)
597 0 : v.v_[0]+=m_[0]; // add
598 0 : break;
599 : }
600 : }
601 :
602 0 : }
603 :
604 :
605 :
606 :
607 : // GLOBALS---------------------------
608 :
609 0 : Mueller* createMueller(const Mueller::MuellerType& mtype) {
610 0 : Mueller *m(NULL);
611 0 : switch (mtype) {
612 0 : case Mueller::General:
613 0 : m = new Mueller();
614 0 : break;
615 0 : case Mueller::Diagonal:
616 0 : m = new MuellerDiag();
617 0 : break;
618 0 : case Mueller::Diag2:
619 0 : m = new MuellerDiag2();
620 0 : break;
621 0 : case Mueller::AddDiag2:
622 0 : m = new AddMuellerDiag2();
623 0 : break;
624 0 : case Mueller::AddDiag:
625 0 : m = new AddMuellerDiag();
626 0 : break;
627 0 : case Mueller::Scalar:
628 0 : m = new MuellerScal();
629 0 : break;
630 : }
631 0 : return m;
632 : }
633 :
634 : // Return the enum for from an int
635 : /*
636 : Mueller::MuellerType muellerType(const Int& n) {
637 : switch (n) {
638 : case 1:
639 : return Mueller::Scalar;
640 : case 2:
641 : return Mueller::Diag2;
642 : case 4:
643 : return Mueller::Diagonal;
644 : case 16:
645 : return Mueller::General;
646 : default:
647 : throw(AipsError("Bad MuellerType."));
648 : }
649 : }
650 : */
651 :
652 0 : Mueller::MuellerType muellerType(const Jones::JonesType& jtype,const VisVector::VisType& vtype) {
653 0 : switch(jtype) {
654 0 : case Jones::General:
655 0 : switch(vtype) {
656 0 : case VisVector::Four:
657 0 : return Mueller::General;
658 : break;
659 0 : default:
660 0 : throw(AipsError("Cannot apply General calibration matrices to incomplete visibilities."));
661 : // TBD: Must cater for the case of some spws with differing VisVector VisTypes
662 : }
663 : break;
664 0 : case Jones::Diagonal:
665 0 : switch(vtype) {
666 0 : case VisVector::Four:
667 0 : return Mueller::Diagonal;
668 : break;
669 0 : default:
670 0 : return Mueller::Diag2;
671 : break;
672 : }
673 : break;
674 0 : case Jones::Scalar:
675 0 : return Mueller::Scalar;
676 : break;
677 0 : default:
678 : // can't reach here
679 0 : throw(AipsError("Cannot figure out Mueller algebra from Jones algebra"));
680 :
681 : }
682 : // Signature consistency (can't reach here)
683 : return Mueller::General;
684 : }
685 :
686 :
687 : // print out a Mueller matrix (according to type)
688 0 : ostream& operator<<(ostream& os, const Mueller& mat) {
689 :
690 : Complex *mi;
691 0 : mi=mat.m_;
692 :
693 0 : switch (mat.type()) {
694 0 : case Mueller::General: {
695 0 : cout << "General Mueller: " << endl;;
696 0 : for (Int i=0;i<4;++i) {
697 0 : cout << " [" << *mi++;
698 0 : for (Int j=1;j<4;++j,++mi) cout << ", " << *mi;
699 0 : cout << "]";
700 0 : if (i<3) cout << endl;
701 : }
702 0 : break;
703 : }
704 0 : case Mueller::Diagonal: {
705 0 : cout << "Diagonal Mueller: " << endl;
706 0 : cout << " [";
707 0 : cout << *mi;
708 0 : ++mi;
709 0 : for (Int i=1;i<4;++i,++mi) cout << ", " << *mi;
710 0 : cout << "]";
711 0 : break;
712 : }
713 0 : case Mueller::Diag2: {
714 0 : cout << "Diag2 Mueller: " << endl;
715 0 : cout << " [" << *mi << ", ";
716 0 : ++mi;
717 0 : cout << *mi << "]";
718 0 : break;
719 : }
720 0 : case Mueller::AddDiag: {
721 0 : cout << "AddDiag Mueller: " << endl;
722 0 : cout << " [";
723 0 : cout << *mi;
724 0 : ++mi;
725 0 : for (Int i=1;i<4;++i,++mi) cout << ", " << *mi;
726 0 : cout << "]";
727 0 : break;
728 : }
729 0 : case Mueller::AddDiag2: {
730 0 : cout << "AddDiag2 Mueller: " << endl;
731 0 : cout << " [" << *mi << ", ";
732 0 : ++mi;
733 0 : cout << *mi << "]";
734 0 : break;
735 : }
736 0 : case Mueller::Scalar: {
737 0 : cout << "Scalar Mueller: " << endl;
738 0 : cout << " [ " << *mi << " ]";
739 0 : break;
740 : }
741 : }
742 :
743 0 : return os;
744 : }
745 :
746 : } //# NAMESPACE CASA - END
747 :
748 :
|
