Line data Source code
1 : //# VisCal.cc: Implementation of VisCal classes
2 : //# Copyright (C) 1996,1997,1998,1999,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 addressed 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 : //# $Id: VisCal.cc,v 1.15 2006/02/06 19:23:11 gmoellen Exp $
27 :
28 : #include <synthesis/MeasurementComponents/VisCal.h>
29 :
30 : #include <synthesis/MeasurementComponents/MSMetaInfoForCal.h>
31 : #include <msvis/MSVis/VisBuffer.h>
32 :
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/casa/Arrays/MaskArrMath.h>
35 : #include <casacore/casa/Arrays/ArrayLogical.h>
36 : #include <casacore/casa/Arrays/ArrayIter.h>
37 : #include <casacore/casa/Arrays/MaskedArray.h>
38 : #include <casacore/scimath/Mathematics/MatrixMathLA.h>
39 : #include <casacore/casa/BasicSL/String.h>
40 : #include <casacore/casa/Utilities/Assert.h>
41 : #include <casacore/casa/Quanta/MVTime.h>
42 : #include <casacore/casa/Quanta/Quantum.h>
43 : #include <casacore/casa/Quanta/QuantumHolder.h>
44 : #include <casacore/casa/Exceptions/Error.h>
45 : #include <casacore/casa/OS/Memory.h>
46 :
47 : #include <sstream>
48 :
49 : #include <casacore/casa/Logging/LogMessage.h>
50 : #include <casacore/casa/Logging/LogSink.h>
51 :
52 : #include <casacore/casa/Quanta/MVTime.h>
53 :
54 : #define PRTLEV 0
55 :
56 : using namespace casacore;
57 : namespace casa { //# NAMESPACE CASA - BEGIN
58 :
59 : // **********************************************************
60 : // VisCal Implementations
61 : //
62 :
63 0 : VisCal::VisCal(VisSet& vs) :
64 : msName_(vs.msName()),
65 0 : msmc_( new MSMetaInfoForCal(vs.msName()) ), // a local one
66 : delmsmc_(True), // must be delete in dtor
67 0 : nSpw_(vs.numberSpw()),
68 0 : nAnt_(vs.numberAnt()),
69 : nBln_(0),
70 : currSpw_(0),
71 0 : currTime_(vs.numberSpw(),0.0),
72 : currScan_(vs.numberSpw(),-1),
73 : currObs_(vs.numberSpw(),-1),
74 : currField_(vs.numberSpw(),-1),
75 : currIntent_(vs.numberSpw(),-1),
76 : currFreq_(vs.numberSpw(),-1),
77 0 : lastTime_(vs.numberSpw(),0.0),
78 : refTime_(0.0),
79 : refFreq_(0.0),
80 : nChanPar_(vs.numberSpw(),1),
81 : nChanMat_(vs.numberSpw(),1),
82 : startChan_(vs.numberSpw(),0),
83 : interval_(0.0),
84 : applied_(false),
85 0 : V_(vs.numberSpw(),NULL),
86 0 : currCPar_(vs.numberSpw(),NULL),
87 0 : currRPar_(vs.numberSpw(),NULL),
88 0 : currParOK_(vs.numberSpw(),NULL),
89 0 : PValid_(vs.numberSpw(),false),
90 : calWt_(false),
91 0 : currWtScale_(vs.numberSpw(),NULL),
92 : prtlev_(PRTLEV),
93 0 : extratag_("")
94 : {
95 0 : if (prtlev()>2) cout << "VC::VC(vs)" << endl;
96 :
97 : // Initialize
98 0 : initVisCal();
99 0 : }
100 :
101 0 : VisCal::VisCal(String msname,Int MSnAnt,Int MSnSpw) :
102 : msName_(msname),
103 0 : msmc_( new MSMetaInfoForCal(msname) ), // a local one
104 : delmsmc_(True), // must be delete in dtor
105 : nSpw_(MSnSpw),
106 : nAnt_(MSnAnt),
107 : nBln_(0),
108 : currSpw_(0),
109 : currTime_(MSnSpw,0.0),
110 : currScan_(MSnSpw,-1),
111 : currObs_(MSnSpw,-1),
112 : currField_(MSnSpw,-1),
113 : currIntent_(MSnSpw,-1),
114 : currFreq_(MSnSpw,-1),
115 : lastTime_(MSnSpw,0.0),
116 : refTime_(0.0),
117 : refFreq_(0.0),
118 : nChanPar_(MSnSpw,1),
119 : nChanMat_(MSnSpw,1),
120 : startChan_(MSnSpw,0),
121 : interval_(0.0),
122 : applied_(false),
123 : V_(MSnSpw,NULL),
124 : currCPar_(MSnSpw,NULL),
125 : currRPar_(MSnSpw,NULL),
126 : currParOK_(MSnSpw,NULL),
127 : PValid_(MSnSpw,false),
128 : calWt_(false),
129 : currWtScale_(MSnSpw,NULL),
130 : prtlev_(PRTLEV),
131 0 : extratag_("")
132 : {
133 0 : if (prtlev()>2) cout << "VC::VC(msname,MSnAnt,MSnSpw)" << endl;
134 :
135 : // Initialize
136 0 : initVisCal();
137 0 : }
138 :
139 0 : VisCal::VisCal(const MSMetaInfoForCal& msmc) :
140 : msName_(msmc.msname()), // from the msmc (it is still used in some places)
141 : msmc_(&msmc),
142 : delmsmc_(False), // not local, don't delete
143 0 : nSpw_(msmc.nSpw()),
144 0 : nAnt_(msmc.nAnt()),
145 : nBln_(0),
146 : currSpw_(0),
147 0 : currTime_(nSpw_,0.0),
148 : currScan_(nSpw_,-1),
149 : currObs_(nSpw_,-1),
150 : currField_(nSpw_,-1),
151 : currIntent_(nSpw_,-1),
152 : currFreq_(nSpw_,-1),
153 0 : lastTime_(nSpw_,0.0),
154 : refTime_(0.0),
155 : refFreq_(0.0),
156 : nChanPar_(nSpw_,1),
157 : nChanMat_(nSpw_,1),
158 : startChan_(nSpw_,0),
159 : interval_(0.0),
160 : applied_(False),
161 0 : V_(nSpw_,NULL),
162 0 : currCPar_(nSpw_,NULL),
163 0 : currRPar_(nSpw_,NULL),
164 0 : currParOK_(nSpw_,NULL),
165 0 : PValid_(nSpw_,False),
166 : calWt_(False),
167 0 : currWtScale_(nSpw_,NULL),
168 : prtlev_(PRTLEV),
169 0 : extratag_("")
170 : {
171 0 : if (prtlev()>2) cout << "VC::VC(MSMetaInfoForCal)" << endl;
172 :
173 : // Initialize
174 0 : initVisCal();
175 0 : }
176 :
177 0 : VisCal::VisCal(const Int& nAnt) :
178 : msName_(""),
179 : msmc_(NULL),
180 : delmsmc_(False),
181 : nSpw_(1),
182 0 : nAnt_(nAnt),
183 : nBln_(0),
184 : currSpw_(0),
185 : currTime_(1,0.0),
186 : currField_(1,-1),
187 : currIntent_(1,-1),
188 : currFreq_(1,-1),
189 : lastTime_(1,0.0),
190 : refTime_(0.0),
191 : refFreq_(0.0),
192 : nChanPar_(1,1),
193 : nChanMat_(1,1),
194 : startChan_(1,0),
195 : interval_(0.0),
196 : applied_(false),
197 : V_(1,NULL),
198 : currCPar_(1,NULL),
199 : currRPar_(1,NULL),
200 : currParOK_(1,NULL),
201 : PValid_(1,false),
202 : calWt_(false),
203 : currWtScale_(1,NULL),
204 : prtlev_(PRTLEV),
205 0 : extratag_("")
206 : {
207 0 : if (prtlev()>2) cout << "VC::VC(i,j,k)" << endl;
208 :
209 : // Initialize
210 0 : initVisCal();
211 0 : }
212 :
213 0 : VisCal::~VisCal() {
214 :
215 0 : if (prtlev()>2) cout << "VC::~VC()" << endl;
216 :
217 0 : deleteVisCal();
218 :
219 0 : if (delmsmc_)
220 0 : delete msmc_;
221 :
222 0 : }
223 :
224 0 : void VisCal::setApply() {
225 :
226 0 : if (prtlev()>2) cout << "VC::setApply()" << endl;
227 :
228 : // This is the apply context
229 0 : setApplied(true);
230 :
231 : // Establish non-trivial paramter arrays
232 0 : for (Int ispw=0;ispw<nSpw();ispw++) {
233 0 : currSpw()=ispw;
234 :
235 0 : if (nChanPar()>0) {
236 0 : switch (parType()) {
237 0 : case VisCalEnum::COMPLEX: {
238 0 : currCPar().resize(nPar(),nChanPar(),nElem());
239 0 : currCPar()=Complex(1.0);
240 0 : break;
241 : }
242 0 : case VisCalEnum::REAL: {
243 0 : currRPar().resize(nPar(),nChanPar(),nElem());
244 0 : currRPar()=1.0;
245 0 : break;
246 : }
247 0 : default:
248 0 : throw(AipsError("Parameters must be entirely Real or entirely Complex for now!"));
249 : }
250 0 : currParOK().resize(nPar(),nChanPar(),nElem());
251 0 : currParOK()=true;
252 : }
253 : }
254 :
255 0 : }
256 :
257 0 : void VisCal::setApply(const Record& apply) {
258 : // Inputs:
259 : // apply Record& Contains application params
260 : //
261 :
262 0 : if (prtlev()>2) cout << "VC::setApply(apply)" << endl;
263 :
264 : // Extract calWt
265 0 : if (apply.isDefined("calwt"))
266 0 : calWt()=apply.asBool("calwt");
267 :
268 : // This is apply context
269 0 : setApplied(true);
270 :
271 : // Initialize flag counting
272 0 : initCalFlagCount();
273 :
274 :
275 0 : }
276 :
277 0 : void VisCal::setCallib(const Record& callib,const MeasurementSet& /*selms*/) {
278 :
279 : // Forward to setApply
280 0 : VisCal::setApply(callib);
281 :
282 0 : }
283 :
284 :
285 0 : String VisCal::applyinfo() {
286 :
287 0 : ostringstream o;
288 0 : o << typeName() << " <pre-computed>";
289 :
290 0 : return String(o);
291 :
292 : }
293 :
294 0 : void VisCal::correct(VisBuffer& vb, Bool trial) {
295 :
296 0 : if (prtlev()>3) cout << "VC::correct(vb)" << endl;
297 :
298 : // Count pre-cal flags
299 0 : countInFlag(vb);
300 :
301 : // Call non-in-place version, in-place-wise:
302 0 : correct(vb,vb.visCube(),trial);
303 :
304 : // Count post-cal flags
305 0 : countOutFlag(vb);
306 :
307 0 : }
308 :
309 0 : void VisCal::correct2(vi::VisBuffer2& vb, Bool trial, Bool doWtSp, Bool dosync) {
310 :
311 0 : if (prtlev()>3) cout << "VC::correct2(vb)" << endl;
312 :
313 :
314 : // If calibration is unavailable, just return!
315 0 : if (!calAvailable(vb)) {
316 0 : return;
317 : }
318 :
319 : // Reaching here, calibration is available, so do it!
320 :
321 : // Count pre-cal flags
322 0 : countInFlag2(vb);
323 :
324 : // Bring calibration up-to-date w/ the vb
325 0 : if (dosync)
326 0 : syncCal2(vb,true);
327 :
328 : // Organize for weight calibration
329 0 : Cube<Float> wtcube;
330 0 : if (doWtSp) {
331 : // refer directly to weightSpectrum
332 0 : wtcube.reference(vb.weightSpectrum());
333 : }
334 : else {
335 : // refer to unchan'd weight matrix, rendered as Cube
336 0 : IPosition sh2=vb.weight().shape();
337 0 : IPosition sh3(3,sh2[0],1,sh2[1]);
338 0 : wtcube.reference(vb.weight().reform(sh3));
339 : }
340 :
341 : // Refer to corrected data
342 0 : Cube<Complex> vCC;
343 0 : vCC.reference(vb.visCubeCorrected());
344 :
345 : // Apply the calibration
346 0 : applyCal2(vb,vCC,wtcube,trial);
347 :
348 : // Count post-cal flags
349 0 : countOutFlag2(vb);
350 :
351 : }
352 :
353 0 : void VisCal::corrupt(VisBuffer& vb) {
354 :
355 0 : if (prtlev()>3) cout << "VC::corrupt(vb)" << endl;
356 :
357 : // Call non-in-place version, in-place-wise:
358 0 : corrupt(vb,vb.modelVisCube());
359 0 : }
360 :
361 0 : void VisCal::corrupt2(vi::VisBuffer2& vb) {
362 :
363 0 : if (prtlev()>3) cout << "VC::corrupt(vb2)" << endl;
364 :
365 : // If calibration is unavailable, just return!
366 0 : if (!calAvailable(vb)) {
367 0 : return;
368 : }
369 :
370 : // Call non-in-place version, in-place-wise:
371 0 : Cube<Complex> m;
372 0 : m.reference(vb.visCubeModel()); // access model non-const
373 0 : corrupt2(vb,m);
374 : }
375 :
376 0 : void VisCal::correct(VisBuffer& vb, Cube<Complex>& Vout, Bool trial) {
377 :
378 0 : if (prtlev()>3) cout << " VC::correct(vb,Vout)" << endl;
379 :
380 : // Prepare output Cube<Complex> for its own in-place apply:
381 : // (this is a no-op if referencing same storage)
382 0 : Vout = vb.visCube();
383 :
384 : // Bring calibration up-to-date with the vb,
385 : // with inversion turned ON
386 0 : syncCal(vb,true);
387 :
388 : // Call generic row-by-row apply, with inversion turned ON
389 0 : applyCal(vb,Vout,trial);
390 :
391 0 : }
392 :
393 : // void VisCal::corrupt(VisBuffer& vb, Cube<Complex>& Mout) {
394 0 : void VisCal::corrupt(VisBuffer& vb, Cube<Complex>& Mout) {
395 :
396 0 : if (prtlev()>3) cout << " VC::corrupt()" << endl;
397 :
398 :
399 : // Ensure weight calibration off internally for corrupt
400 : // (corruption doesn't re-scale the data!)
401 0 : Bool userCalWt=calWt();
402 0 : calWt()=false;
403 :
404 : // Prepare output Cube<Complex> for its own in-place apply:
405 : // (this is a no-op if referencing same storage)
406 0 : Mout = vb.modelVisCube();
407 :
408 : // Bring calibration up-to-date with the vb,
409 : // with inversion turned OFF
410 0 : syncCal(vb,false);
411 :
412 : // Call generic row-by-row apply, with inversion turned OFF
413 0 : applyCal(vb,Mout);
414 :
415 : // Restore user's calWt()
416 0 : calWt()=userCalWt;
417 :
418 0 : }
419 :
420 0 : void VisCal::corrupt2(vi::VisBuffer2& vb, Cube<Complex>& Mout) {
421 :
422 0 : if (prtlev()>3) cout << " VC::corrupt2(VB2)" << endl;
423 :
424 :
425 : // Ensure weight calibration off internally for corrupt
426 : // (corruption doesn't re-scale the data!)
427 0 : Bool userCalWt=calWt();
428 0 : calWt()=False;
429 :
430 : // Prepare output Cube<Complex> for its own in-place apply:
431 : // (this is a no-op if referencing same storage)
432 0 : Mout = vb.visCubeModel();
433 :
434 : // Bring calibration up-to-date with the vb,
435 : // with inversion turned OFF
436 0 : syncCal2(vb,False);
437 :
438 : // Call generic row-by-row apply, with inversion turned OFF
439 0 : Cube<Float> w0; // place-holder
440 0 : applyCal2(vb,Mout,w0);
441 :
442 : // Restore user's calWt()
443 0 : calWt()=userCalWt;
444 :
445 0 : }
446 :
447 : // VI2-related refactorings-----------
448 :
449 0 : void VisCal::setMeta(Int obs, Int scan, Double time,
450 : Int spw, const Vector<Double>& freq,
451 : Int fld) {
452 :
453 : // Call internal method
454 0 : syncMeta(spw,time,fld,freq,freq.nelements());
455 0 : currScan()=scan;
456 0 : currObs()=obs;
457 :
458 0 : }
459 :
460 : // Setup solvePar shape for a spw
461 0 : void VisCal::sizeApplyParCurrSpw(Int nVisChan) {
462 :
463 : // Sizes the solvePar arrays for the currSpw()
464 :
465 0 : if (prtlev()>3) cout << " VC::sizeApplyParCurrSpw()" << endl;
466 :
467 : // Use nVisChan only for freqDepPar() types
468 0 : Int nChan = ( freqDepPar() ? nVisChan : 1);
469 :
470 : // Keep old way informed (needed?)
471 0 : nChanPar()=nChan;
472 :
473 : // Now, size the arrays:
474 :
475 0 : IPosition parsh(3,nPar(),nChan,nElem()); // multi-chan
476 0 : switch (parType()) {
477 0 : case VisCalEnum::COMPLEX: {
478 0 : currCPar().resize(parsh);
479 0 : currCPar()=Complex(1.0);
480 0 : break;
481 : }
482 0 : case VisCalEnum::REAL: {
483 0 : currRPar().resize(parsh);
484 0 : currRPar()=0.0;
485 0 : break;
486 : }
487 0 : default:
488 : throw(AipsError("Internal error(Calibrater Module): Unsupported parameter type "
489 0 : "COMPLEXREAL found in VC::sizeApplyParCurrSpw()"));
490 : }
491 :
492 0 : currParOK().resize(parsh);
493 0 : currParOK()=True;
494 :
495 0 : }
496 :
497 0 : void VisCal::setDefApplyParCurrSpw(Bool sync, Bool doInv) {
498 :
499 : // TBD: generalize for type-dep def values, etc.
500 0 : switch (parType()) {
501 0 : case VisCalEnum::COMPLEX: {
502 0 : AlwaysAssert(currCPar().nelements()>0,AipsError);
503 0 : currCPar().set(1.0);
504 0 : break;
505 : }
506 0 : case VisCalEnum::REAL: {
507 0 : AlwaysAssert(currRPar().nelements()>0,AipsError);
508 0 : currRPar().set(0.0);
509 0 : break;
510 : }
511 0 : default:
512 : throw(AipsError("Internal error(Calibrater Module): Unsupported parameter type "
513 0 : "COMPLEXREAL found in VC::setDefApplyParCurrSpw()"));
514 : }
515 :
516 0 : currParOK().set(True);
517 0 : validateP();
518 :
519 0 : if (sync)
520 0 : syncCal(doInv);
521 :
522 0 : }
523 :
524 : // Set parameters to specified values in the currSpw(),
525 : // and optionally sync matrices
526 0 : void VisCal::setApplyParCurrSpw(const Cube<Complex> cpar,
527 : bool sync, bool doInv) {
528 :
529 0 : switch (parType()) {
530 0 : case VisCalEnum::COMPLEX: {
531 0 : AlwaysAssert(currCPar().nelements()>0,AipsError);
532 0 : AlwaysAssert(currCPar().shape()==cpar.shape(),AipsError); // shapes must match
533 0 : currCPar()=cpar;
534 0 : break;
535 : }
536 0 : case VisCalEnum::REAL: {
537 0 : throw(AipsError("Cannot set real parameters with complext values!"));
538 : break;
539 : }
540 0 : default:
541 : throw(AipsError("Internal error(Calibrater Module): Unsupported parameter type "
542 0 : "COMPLEXREAL found in VC::setDefApplyParCurrSpw()"));
543 : }
544 :
545 0 : currParOK().set(True);
546 0 : validateP();
547 :
548 : // Synchronize matrices
549 0 : if (sync)
550 0 : syncCal(doInv);
551 :
552 0 : }
553 :
554 0 : void VisCal::setApplyParCurrSpw(const Cube<float> rpar,
555 : bool sync, bool doInv) {
556 :
557 0 : switch (parType()) {
558 0 : case VisCalEnum::REAL: {
559 0 : AlwaysAssert(currRPar().nelements()>0,AipsError);
560 0 : AlwaysAssert(currRPar().shape()==rpar.shape(),AipsError); // shapes must match
561 0 : currRPar()=rpar;
562 0 : break;
563 : }
564 0 : case VisCalEnum::COMPLEX: {
565 0 : throw(AipsError("Cannot set complex parameters with real values!"));
566 : break;
567 : }
568 0 : default:
569 : throw(AipsError("Internal error(Calibrater Module): Unsupported parameter type "
570 0 : "COMPLEXREAL found in VC::setDefApplyParCurrSpw()"));
571 : }
572 :
573 0 : currParOK().set(True);
574 0 : validateP();
575 :
576 : // Synchronize matrices
577 0 : if (sync)
578 0 : syncCal(doInv);
579 :
580 0 : }
581 :
582 :
583 0 : void VisCal:: initCalFlagCount() {
584 0 : ndataIn_=nflagIn_=nflagOut_=0;
585 0 : }
586 :
587 0 : Record VisCal::actionRec() {
588 0 : Record cf;
589 0 : cf.define("type",typeName());
590 0 : if (isApplied()) {
591 0 : cf.define("ndata",ndataIn_);
592 0 : cf.define("nflagIn",nflagIn_);
593 0 : cf.define("nflagOut",nflagOut_);
594 : }
595 0 : return cf;
596 : }
597 :
598 :
599 0 : void VisCal::state() {
600 :
601 0 : cout << "===========================================================" << endl;
602 0 : if (prtlev()>3) cout << "VC::state():" << endl;
603 0 : cout << boolalpha;
604 0 : cout << " Type=" << typeName() << " Enum=" << type() << endl;
605 0 : cout << " nSpw= " << nSpw()
606 0 : << "; nAnt= " << nAnt()
607 0 : << "; nBln= " << nBln()
608 0 : << "; nPar= " << nPar()
609 0 : << "; nElem= " << nElem()
610 0 : << "; nCalMat= " << nCalMat() << endl;
611 0 : cout << " isApplied= " << isApplied()
612 0 : << "; isSolvable= " << isSolvable() << endl;
613 0 : cout << " freqDepPar= " << freqDepPar()
614 0 : << "; freqDepMat= " << freqDepMat()
615 0 : << "; timeDepMat= " << timeDepMat() << endl;
616 0 : cout << " interval= " << interval() << endl;
617 0 : cout << " currSpw= " << currSpw() << endl;
618 0 : cout << " nChanPar= " << nChanPar()
619 0 : << "; nChanMat= " << nChanMat()
620 0 : << "; startChan= " << startChan() << endl;
621 0 : cout << " lastTime= " << lastTime()
622 0 : << "; currTime= " << currTime()
623 0 : << "; refTime= " << refTime() << endl;
624 0 : cout << " refFreq= " << refFreq() << endl;
625 0 : cout << " currField=" << currField() << endl;
626 0 : cout << " currIntent=" << currIntent() << endl;
627 0 : cout << " PValid() = " << PValid() << endl;
628 0 : cout << " currCPar().shape() = " << currCPar().shape()
629 0 : << " (" << currCPar().data() << ")" << endl;
630 0 : cout << " currRPar().shape() = " << currRPar().shape()
631 0 : << " (" << currRPar().data() << ")" << endl;
632 0 : cout << " currParOK().shape() = " << currParOK().shape()
633 0 : << " (" << currParOK().data() << ") "
634 0 : << " (ntrue=" << ntrue(currParOK()) << ")" << endl;
635 :
636 0 : cout << "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" << endl;
637 :
638 0 : }
639 :
640 0 : void VisCal::currMetaNote() {
641 :
642 : cout << " ("
643 0 : << "time=" << MVTime(refTime()/C::day).string(MVTime::YMD,7)
644 0 : << " field=" << currField()
645 0 : << " spw=" << currSpw()
646 0 : << ")"
647 0 : << endl;
648 :
649 0 : }
650 :
651 : // VisCal PROTECTED:
652 :
653 :
654 0 : void VisCal::countInFlag(const VisBuffer& vb) {
655 0 : Int ncorr=vb.nCorr();
656 0 : ndataIn_+=(vb.flag().nelements()*ncorr);
657 0 : nflagIn_+=(ntrue(vb.flag())*ncorr);
658 0 : }
659 0 : void VisCal::countInFlag2(const vi::VisBuffer2& vb) {
660 0 : ndataIn_+=(vb.flagCube().nelements());
661 0 : nflagIn_+=(ntrue(vb.flagCube()));
662 0 : }
663 :
664 0 : void VisCal::countOutFlag(const VisBuffer& vb){
665 0 : nflagOut_+=ntrue(vb.flag())*vb.nCorr();
666 0 : }
667 :
668 0 : void VisCal::countOutFlag2(const vi::VisBuffer2& vb){
669 0 : nflagOut_+=ntrue(vb.flagCube());
670 0 : }
671 :
672 0 : void VisCal::syncCal(const VisBuffer& vb,
673 : const Bool& doInv) {
674 :
675 0 : if (prtlev()>4) cout << " VC::syncCal(vb)" << endl;
676 :
677 : // Set current meta date from the VisBuffer
678 0 : syncMeta(vb);
679 :
680 : // Check the current calibration for validity
681 0 : checkCurrCal();
682 :
683 : // Procede with generalized sync of calibration
684 0 : syncCal(doInv);
685 :
686 0 : }
687 :
688 0 : void VisCal::syncCal2(const vi::VisBuffer2& vb,
689 : const Bool& doInv) {
690 :
691 0 : if (prtlev()>4) cout << " VC::syncCal2(vb)" << endl;
692 :
693 : // Set current meta date from the VisBuffer
694 0 : syncMeta2(vb);
695 :
696 : // Check the current calibration for validity
697 0 : checkCurrCal();
698 :
699 : // Procede with generalized sync of calibration
700 0 : syncCal(doInv);
701 :
702 :
703 0 : }
704 :
705 :
706 0 : void VisCal::syncCal(VisCal& vc) {
707 :
708 0 : if (prtlev()>4) cout << " VC::syncCal(vc)" << endl;
709 :
710 : // cout << " VisCal::syncCal(VisCal): " << currCPar().data()
711 : // << endl;
712 :
713 : // Set current meta date from the VisCal
714 0 : syncMeta(vc);
715 :
716 : // Check the current calibration for validity
717 0 : checkCurrCal();
718 :
719 : // Procede with generalized sync of calibration
720 0 : syncCal(false);
721 :
722 : // cout << " VisCal::syncCal(VisCal): " << currCPar().data()
723 : // << endl;
724 :
725 0 : }
726 :
727 :
728 0 : void VisCal::syncMeta(const VisBuffer& vb) {
729 :
730 0 : if (prtlev()>4) cout << " VC::syncMeta(vb)" << endl;
731 :
732 0 : syncMeta(vb.spectralWindow(),
733 0 : vb.time()(0),
734 0 : vb.fieldId(),
735 0 : vb.frequency(),
736 0 : vb.nChannel());
737 :
738 0 : currObs()=vb.observationId()(0);
739 :
740 :
741 : // Ensure VisVector for data acces has correct form
742 0 : Int ncorr(vb.corrType().nelements());
743 0 : if (V().type() != ncorr)
744 0 : V().setType(visType(ncorr));
745 :
746 0 : }
747 :
748 0 : void VisCal::syncMeta2(const vi::VisBuffer2& vb) {
749 :
750 0 : if (prtlev()>4) cout << " VC::syncMeta2(vb)" << endl;
751 :
752 0 : syncMeta(vb.spectralWindows()(0),
753 0 : vb.time()(0),
754 0 : vb.fieldId()(0),
755 0 : vb.getFrequencies(0),
756 0 : vb.nChannels());
757 :
758 0 : currScan()=vb.scan()(0);
759 0 : currObs()=vb.observationId()(0);
760 0 : currIntent()=vb.stateId()(0);
761 :
762 : // Ensure VisVector for data acces has correct form
763 0 : Int ncorr(vb.nCorrelations());
764 0 : if (V().type() != ncorr)
765 0 : V().setType(visType(ncorr));
766 :
767 0 : }
768 :
769 0 : void VisCal::syncMeta(VisCal& vc) {
770 :
771 0 : if (prtlev()>4) cout << " VC::syncMeta(vc)" << endl;
772 :
773 :
774 0 : syncMeta(vc.currSpw(),
775 0 : vc.currTime(),
776 0 : vc.currField(),
777 0 : vc.currFreq(),
778 0 : vc.nChanMat()); // the number of channels to be applied to
779 :
780 0 : currObs()=vc.currObs();
781 :
782 :
783 : // TBD: Ensure nElem matches (here?)
784 :
785 0 : }
786 :
787 0 : void VisCal::syncMeta(const Int& spw,
788 : const Double& time,
789 : const Int& field,
790 : const Vector<Double>& freq,
791 : const Int& nchan) {
792 :
793 0 : if (prtlev()>4) cout << " VC::syncMeta(...)" << endl;
794 :
795 : // Remember which spw this is---EVERYTHING below pivots on this
796 0 : currSpw()=spw;
797 :
798 : // Remember other meta-data
799 0 : currTime()=time;
800 0 : currField()=field;
801 :
802 0 : currFreq().resize();
803 0 : currFreq()=freq;
804 0 : currFreq()/=1.0e9; // In GHz
805 :
806 : // Use center channel for now
807 0 : refFreq()=currFreq()(nchan/2);
808 :
809 : // Ensure cal matrix channelization is correct vis-a-vis the target channelization
810 : // (in solve context, this is handled differently, outside this method)
811 0 : if (isApplied()) setCalChannelization(nchan);
812 :
813 0 : if (prtlev()>5) cout << " meta: t,fld,freq=" << time << field << freq << endl;
814 0 : }
815 :
816 0 : void VisCal::setCalChannelization(const Int& nChanDat) {
817 :
818 0 : if (prtlev()>4) cout << " VC::setCalChannelization()" << endl;
819 :
820 : // This method sets up cal Matrix channelization for the current
821 : // spw according to the channel count supplied by the caller.
822 : // (Presumably the data channelization)
823 :
824 : // Ensure matrix channel axis length is correct
825 0 : if (freqDepMat())
826 :
827 : // If parameters are channelized
828 0 : if (freqDepPar())
829 : // follow parameter channelization
830 : // (e.g., B)
831 : // nChanMat() = nChanPar();
832 0 : nChanMat() = nChanDat; // NEWCALTABLE
833 : else
834 : // cal is f(freq) for each data channel:
835 : // (e.g., fringe-fitting)
836 0 : nChanMat() = nChanDat;
837 :
838 : else {
839 : // Matrices are single channel
840 : // (e.g., G)
841 0 : nChanMat()=1;
842 :
843 : // cout << "nChanPar() = " << nChanPar() << endl;
844 :
845 : // So must be parameters:
846 0 : AlwaysAssert((!freqDepPar()),AipsError);
847 : // NCT AlwaysAssert((nChanPar()==1),AipsError);
848 :
849 : }
850 :
851 0 : }
852 :
853 0 : void VisCal::checkCurrCal() {
854 : // Based on meta-data changes, determine mainly if
855 : // new calibration PARAMETERS should be sought
856 : // NB: A finding of "true" does not necessarily
857 : // mean new pars will be found by CalInterp (we are
858 : // here only forcing it to try)
859 : // NB: this method may also invalidateCalMat(), to
860 : // trigger matrix reformation, independent of whether
861 : // new pars are found ultimately found by CalInterp
862 :
863 0 : if (prtlev()>4) cout << " VC::checkCurrCal: " << endl;
864 :
865 : // We potentially need new calibration (for currSpw) IF...
866 :
867 : // ...timestamp has changed
868 0 : if (currTime() != lastTime() ) {
869 0 : lastTime()=currTime();
870 0 : invalidateP(); // merely signals need to check (doesn't break reference)
871 :
872 : // If matrices are time-dependent within solution interval,
873 : // a time change ALWAYS invalidates them
874 : // (i.e., PValid=F may not yield _new_ parameters in
875 : // calcPar, but this forces matrix re-sync anyway)
876 0 : if (timeDepMat()) invalidateCalMat();
877 :
878 : //if (prtlev()>5 and timeDepMat()) cout << " invalidateCalMat() " << endl;
879 :
880 : }
881 :
882 0 : }
883 :
884 0 : void VisCal::syncCal(const Bool& doInv) {
885 :
886 0 : if (prtlev()>4) cout << " VC::syncCal(doInv)" << endl;
887 :
888 : // cout << " VisCal::syncCal(doInv): " << currCPar().data()
889 : // << endl;
890 :
891 : // Synchronize the parameters
892 0 : syncPar();
893 :
894 : // cout << " VisCal::syncCal(doInv): " << currCPar().data()
895 : // << endl;
896 :
897 0 : if (!PValid())
898 : // TBD: Improve this message (with some meta data)
899 0 : throw(AipsError("Could not find any calibration parameters."));
900 :
901 : // Synchronize the matrices
902 0 : syncCalMat(doInv);
903 :
904 : // cout << " VisCal::syncCal(doInv): " << currCPar().data()
905 : // << endl;
906 :
907 0 : }
908 :
909 :
910 :
911 0 : void VisCal::syncPar() {
912 :
913 0 : if (prtlev()>5) cout << " VC::syncPar() (PValid()="
914 0 : << PValid() << ")" << endl;
915 :
916 : // If we require new parameters, calculate them by some means (e.g., interpolation):
917 0 : if (!PValid()) calcPar();
918 :
919 0 : }
920 :
921 0 : void VisCal::calcPar() {
922 :
923 0 : if (prtlev()>6) cout << " VC::calcPar()" << endl;
924 :
925 : // Ensure we have some parameters
926 0 : if (parType()==VisCalEnum::COMPLEX && (currCPar().shape()!=IPosition(3,nPar(),nChanPar(),nElem()) ||
927 0 : currParOK().shape()!=IPosition(3,nPar(),nChanPar(),nElem())) ) {
928 0 : cout << "currCPar() = " << currCPar() << endl;
929 0 : cout << "currParOK() = " << currParOK() << endl;
930 0 : throw(AipsError("No (complex) parameters available!"));
931 : }
932 0 : else if (parType()==VisCalEnum::REAL && (currRPar().shape()!=IPosition(3,nPar(),nChanPar(),nElem()) ||
933 0 : currParOK().shape()!=IPosition(3,nPar(),nChanPar(),nElem())) ) {
934 0 : cout << "currRPar() = " << currRPar() << endl;
935 0 : cout << "currParOK() = " << currParOK() << endl;
936 0 : throw(AipsError("No (real) parameters available!"));
937 : }
938 0 : else if (parType()==VisCalEnum::COMPLEXREAL)
939 0 : throw(AipsError("We can't handle combined Real and Complex parameters yet."));
940 : else
941 : // Parameters appear to be available, so signal validation
942 0 : validateP();
943 :
944 : // Force subsequent matrix calculation
945 : // (Specializations will only do this if actual _new_ parameters
946 : // have been found)
947 0 : invalidateCalMat();
948 :
949 0 : }
950 :
951 :
952 : // VisCal PRIVATE method implementations...........................
953 :
954 0 : void VisCal::initVisCal() {
955 :
956 0 : if (prtlev()>2) cout << " VC::initVisCal()" << endl;
957 :
958 : // Number of baselines (including ACs)
959 0 : nBln_ = nAnt_*(nAnt_+1)/2;
960 :
961 0 : for (Int ispw=0;ispw<nSpw(); ispw++) {
962 0 : currCPar_[ispw] = new Cube<Complex>();
963 0 : currRPar_[ispw] = new Cube<Float>();
964 0 : currParOK_[ispw] = new Cube<Bool>();
965 0 : V_[ispw] = new VisVector(VisVector::Four);
966 0 : currWtScale_[ispw] = new Cube<Float>();
967 : }
968 :
969 0 : }
970 :
971 0 : void VisCal::deleteVisCal() {
972 :
973 0 : if (prtlev()>2) cout << " VC::deleteVisCal()" << endl;
974 :
975 0 : for (Int ispw=0; ispw<nSpw(); ispw++) {
976 0 : if (currCPar_[ispw]) delete currCPar_[ispw];
977 0 : if (currRPar_[ispw]) delete currRPar_[ispw];
978 0 : if (currParOK_[ispw]) delete currParOK_[ispw];
979 0 : if (V_[ispw]) delete V_[ispw];
980 0 : if (currWtScale_[ispw]) delete currWtScale_[ispw];
981 : }
982 0 : currCPar_=NULL;
983 0 : currRPar_=NULL;
984 0 : currParOK_=NULL;
985 0 : V_=NULL;
986 0 : currWtScale_=NULL;
987 0 : }
988 :
989 0 : void VisCal::setCurrField(const Int& ifld) {
990 0 : currField_.set(ifld);
991 0 : }
992 :
993 : // **********************************************************
994 : // VisMueller Implementations
995 : //
996 :
997 0 : VisMueller::VisMueller(VisSet& vs) :
998 : VisCal(vs),
999 0 : M_(vs.numberSpw(),NULL),
1000 0 : currMElem_(vs.numberSpw(),NULL),
1001 0 : currMElemOK_(vs.numberSpw(),NULL),
1002 0 : MValid_(vs.numberSpw(),false)
1003 : {
1004 :
1005 0 : if (prtlev()>2) cout << "VM::VM(vs)" << endl;
1006 :
1007 0 : initVisMueller();
1008 0 : }
1009 :
1010 0 : VisMueller::VisMueller(String msname,Int MSnAnt,Int MSnSpw) :
1011 : VisCal(msname,MSnAnt,MSnSpw),
1012 : M_(MSnSpw,NULL),
1013 : currMElem_(MSnSpw,NULL),
1014 : currMElemOK_(MSnSpw,NULL),
1015 0 : MValid_(MSnSpw,false)
1016 : {
1017 :
1018 0 : if (prtlev()>2) cout << "VM::VM(msname,MSnAnt,MSnSpw)" << endl;
1019 :
1020 0 : initVisMueller();
1021 0 : }
1022 :
1023 0 : VisMueller::VisMueller(const MSMetaInfoForCal& msmc) :
1024 : VisCal(msmc),
1025 0 : M_(nSpw(),NULL),
1026 0 : currMElem_(nSpw(),NULL),
1027 0 : currMElemOK_(nSpw(),NULL),
1028 0 : MValid_(nSpw(),False)
1029 : {
1030 :
1031 0 : if (prtlev()>2) cout << "VM::VM(MSMetaInfoForCal)" << endl;
1032 :
1033 0 : initVisMueller();
1034 0 : }
1035 :
1036 :
1037 0 : VisMueller::VisMueller(const Int& nAnt) :
1038 : VisCal(nAnt),
1039 : M_(1,NULL),
1040 : currMElem_(1,NULL),
1041 : currMElemOK_(1,NULL),
1042 0 : MValid_(1,false)
1043 : {
1044 0 : if (prtlev()>2) cout << "VM::VM(i,j,k)" << endl;
1045 :
1046 0 : initVisMueller();
1047 0 : }
1048 :
1049 0 : VisMueller::~VisMueller() {
1050 0 : if (prtlev()>2) cout << "VM::~VM()" << endl;
1051 :
1052 0 : deleteVisMueller();
1053 0 : }
1054 :
1055 : // Report the VisMueller portion of the state
1056 0 : void VisMueller::state() {
1057 :
1058 : // Call parent
1059 0 : VisCal::state();
1060 :
1061 0 : if (applyByMueller()) {
1062 0 : if (prtlev()>3) cout << "VM::state()" << endl;
1063 0 : cout << boolalpha;
1064 : // cout << "This is an baseline-based (Mueller) calibration." << endl;
1065 0 : cout << " applyByMueller=" << applyByMueller() << endl;
1066 0 : cout << " muellerType= " << muellerType() << endl;
1067 0 : cout << " trivialMuellerElem= " << trivialMuellerElem() << endl;
1068 0 : cout << " MValid() = " << MValid() << endl;
1069 :
1070 0 : cout << " currMElem().shape() = " << currMElem().shape()
1071 0 : << " (" << currMElem().data() << ")" << endl;
1072 0 : cout << " currMElemOK().shape() = " << currMElemOK().shape()
1073 0 : << " (" << currMElemOK().data() << ") "
1074 0 : << " (ntrue=" << ntrue(currMElemOK()) << ")" << endl;
1075 :
1076 :
1077 0 : cout << "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" << endl;
1078 : }
1079 0 : }
1080 :
1081 : // Apply this calibration to VisBuffer visibilities
1082 0 : void VisMueller::applyCal(VisBuffer& vb, Cube<Complex>& Vout,
1083 : Bool trial) {
1084 :
1085 0 : if (prtlev()>3) cout << " VM::applyCal()" << endl;
1086 :
1087 : // CURRENTLY ASSUMES ONLY *ONE* TIMESTAMP IN VISBUFFER
1088 : /*
1089 : cout << "VM::applyCal: type= " << typeName() << " trial = "
1090 : << boolalpha << trial << " calWt = " << calWt()
1091 : << " freqDepPar() = " << freqDepPar()
1092 : << " freqDepMat() = " << freqDepMat()
1093 : << endl;
1094 : */
1095 :
1096 : // Data info/indices
1097 : Int* dataChan;
1098 0 : Bool* flagR=&vb.flagRow()(0);
1099 0 : Bool* flag=&vb.flag()(0,0);
1100 0 : Int* a1=&vb.antenna1()(0);
1101 0 : Int* a2=&vb.antenna2()(0);
1102 0 : Matrix<Float> wtmat;
1103 :
1104 : // Access to weights
1105 0 : if (calWt() && !trial)
1106 0 : wtmat.reference(vb.weightMat());
1107 :
1108 0 : ArrayIterator<Float> wt(wtmat,1);
1109 0 : Vector<Float> wtvec;
1110 :
1111 0 : if (V().type()==VisVector::One) {
1112 0 : M().setScalarData(true);
1113 : }
1114 : else
1115 0 : M().setScalarData(false);
1116 :
1117 : // iterate rows
1118 0 : Int& nRow(vb.nRow());
1119 0 : Int& nChanDat(vb.nChannel());
1120 : Int ibln;
1121 0 : for (Int row=0; row<nRow; row++,flagR++,a1++,a2++) {
1122 :
1123 : // The basline number
1124 0 : ibln=blnidx(*a1,*a2);
1125 :
1126 : // Solution channel registration
1127 0 : Int solCh0(0);
1128 0 : dataChan=&vb.channel()(0);
1129 :
1130 : // If cal _parameters_ are not freqDep (e.g., a delay)
1131 : // the startChan() should be the same as the first data channel
1132 0 : if (freqDepMat() && !freqDepPar())
1133 0 : startChan()=(*dataChan);
1134 :
1135 : // Solution and data array registration
1136 0 : M().sync(currMElem()(0,solCh0,ibln),currMElemOK()(0,solCh0,ibln));
1137 0 : if (!trial)
1138 0 : V().sync(Vout(0,0,row));
1139 :
1140 :
1141 0 : for (Int chn=0; chn<nChanDat; chn++,flag++,V()++,dataChan++) {
1142 :
1143 0 : if (trial)
1144 0 : M().applyFlag(*flag);
1145 : else
1146 : // data and solution ok, do the apply
1147 0 : M().apply(V(),*flag);
1148 :
1149 : // inc soln ch axis if freq-dependent
1150 0 : if (freqDepMat())
1151 0 : M()++;
1152 :
1153 : } // chn
1154 :
1155 : // If requested update the weights
1156 0 : if (calWt() && !trial) {
1157 0 : wtvec.reference(wt.array());
1158 0 : updateWt(wtvec,*a1,*a2);
1159 : }
1160 :
1161 0 : if (calWt() && !trial)
1162 0 : wt.next();
1163 :
1164 : }
1165 :
1166 0 : }
1167 :
1168 : // Apply this calibration to VisBuffer visibilities
1169 0 : void VisMueller::applyCal2(vi::VisBuffer2& vb,
1170 : Cube<Complex>& Vout, Cube<Float>& Wout,
1171 : Bool trial) {
1172 :
1173 0 : if (prtlev()>3) cout << " VM::applyCal()" << endl;
1174 :
1175 : // CURRENTLY ASSUMES ONLY *ONE* TIMESTAMP IN VISBUFFER
1176 : /*
1177 : cout << "VM::applyCal: type= " << typeName() << " trial = "
1178 : << boolalpha << trial << " calWt = " << calWt()
1179 : << " freqDepPar() = " << freqDepPar()
1180 : << " freqDepMat() = " << freqDepMat()
1181 : << endl;
1182 : */
1183 :
1184 0 : Vector<Bool> flagRv(vb.flagRow());
1185 0 : Vector<Int> a1v(vb.antenna1());
1186 0 : Vector<Int> a2v(vb.antenna2());
1187 0 : Cube<Bool> flagCube(vb.flagCube());
1188 0 : Cube<Complex> visCube(Vout);
1189 0 : ArrayIterator<Float> wt(Wout,2);
1190 0 : Matrix<Float> wtmat;
1191 :
1192 : // Data info/indices
1193 : Int* dataChan;
1194 0 : Bool* flagR=&flagRv(0);
1195 0 : Int* a1=&a1v(0);
1196 0 : Int* a2=&a2v(0);
1197 :
1198 : // Ensure VisVector for data acces has correct form
1199 0 : Int ncorr(vb.nCorrelations());
1200 0 : if (V().type() != ncorr)
1201 0 : V().setType(visType(ncorr));
1202 :
1203 0 : if (V().type()==VisVector::One) {
1204 0 : M().setScalarData(true);
1205 : }
1206 : else
1207 0 : M().setScalarData(false);
1208 :
1209 : // iterate rows
1210 0 : Int nRow=vb.nRows();
1211 0 : Int nChanDat=vb.nChannels();
1212 : Int ibln;
1213 0 : Vector<Int> dataChanv(vb.getChannelNumbers(0)); // All rows have same chans
1214 0 : for (Int row=0; row<nRow; row++,flagR++,a1++,a2++) {
1215 :
1216 : // The basline number
1217 0 : ibln=blnidx(*a1,*a2);
1218 :
1219 : // Solution channel registration
1220 0 : Int solCh0(0);
1221 0 : dataChan=&dataChanv(0);
1222 :
1223 : // If cal _parameters_ are not freqDep (e.g., a delay)
1224 : // the startChan() should be the same as the first data channel
1225 0 : if (freqDepMat() && !freqDepPar())
1226 0 : startChan()=(*dataChan);
1227 :
1228 : // Solution and data array registration
1229 0 : M().sync(currMElem()(0,solCh0,ibln),currMElemOK()(0,solCh0,ibln));
1230 0 : V().sync(visCube(0,0,row),flagCube(0,0,row));
1231 :
1232 0 : for (Int chn=0; chn<nChanDat; chn++,V()++,dataChan++) {
1233 :
1234 0 : if (trial)
1235 0 : M().flag(V());
1236 : else
1237 : // data and solution ok, do the apply
1238 0 : M().apply(V());
1239 :
1240 : // inc soln ch axis if freq-dependent
1241 0 : if (freqDepMat())
1242 0 : M()++;
1243 :
1244 : } // chn
1245 :
1246 : // If requested update the weights
1247 : /*
1248 : if (calWt() && !trial) {
1249 : wtvec.reference(wt.array());
1250 : updateWt(wtvec,*a1,*a2);
1251 : }
1252 :
1253 : if (calWt() && !trial)
1254 : wt.next();
1255 : */
1256 : }
1257 :
1258 0 : }
1259 :
1260 0 : void VisMueller::syncCalMat(const Bool& doInv) {
1261 :
1262 0 : if (prtlev()>5) cout << " VM::syncCalMat()"
1263 0 : << " (MValid()=" << MValid() << ")" << endl;
1264 :
1265 : // If Mueller matrices now invalid, re-sync them
1266 0 : if (!MValid()) syncMueller(doInv);
1267 :
1268 0 : }
1269 :
1270 0 : void VisMueller::syncMueller(const Bool& doInv) {
1271 :
1272 : // RI
1273 : //prtlev()=10;
1274 :
1275 0 : if (prtlev()>6) cout << " VM::syncMueller()" << endl;
1276 :
1277 : // If Mueller pars are just the matrix elements:
1278 0 : if (trivialMuellerElem()) {
1279 : // Matrix Elem cache references par cache
1280 0 : currMElem().reference(currCPar());
1281 0 : currMElemOK().reference(currParOK());
1282 :
1283 : // EXCEPT, ensure uniqueness if taking the inverse
1284 : // (this makes a copy, can we avoid?)
1285 0 : if (doInv) currMElem().unique();
1286 : }
1287 : else {
1288 : // Make local matrix element cache the correct size:
1289 0 : currMElem().resize(muellerNPar(this->muellerType()),nChanMat(),nCalMat());
1290 0 : currMElem().unique(); // Ensure uniqueness!
1291 :
1292 : // OK is the shape of the M matrix itself
1293 0 : currMElemOK().resize(muellerNPar(this->muellerType()),nChanMat(),nCalMat());
1294 0 : currMElemOK().unique();
1295 0 : currMElemOK()=false;
1296 :
1297 : // The matrix state is invalid until we actually calculate them
1298 0 : invalidateM();
1299 :
1300 : // And calculate the matrix elements for all baselines
1301 0 : calcAllMueller();
1302 :
1303 : }
1304 :
1305 : // weight calibration
1306 0 : if (calWt()) syncWtScale();
1307 :
1308 : // Ensure Mueller matrix renderer is correct
1309 0 : createMueller();
1310 :
1311 : // Invert, if requested
1312 0 : if (doInv) invMueller();
1313 :
1314 : // Set matrix elements by their ok flags
1315 0 : setMatByOk();
1316 :
1317 : // Mueller matrices now valid
1318 0 : validateM();
1319 :
1320 0 : }
1321 :
1322 : // Calculate Mueller matrices for all baselines
1323 0 : void VisMueller::calcAllMueller() {
1324 :
1325 0 : if (prtlev()>6) cout << " VM::calcAllMueller" << endl;
1326 :
1327 : // Should handle OK flags in this method, and only
1328 : // do Mueller calc if OK
1329 :
1330 0 : Vector<Complex> oneMueller;
1331 0 : Vector<Bool> oneMOK;
1332 0 : Vector<Complex> onePar;
1333 0 : Vector<Bool> onePOK;
1334 :
1335 0 : ArrayIterator<Complex> Miter(currMElem(),1);
1336 0 : ArrayIterator<Bool> MOKiter(currMElemOK(),1);
1337 0 : ArrayIterator<Complex> Piter(currCPar(),1);
1338 0 : ArrayIterator<Bool> POKiter(currParOK(),1);
1339 :
1340 : // All required baselines
1341 0 : for (Int ibln=0; ibln<nCalMat(); ibln++) {
1342 :
1343 : // The following assumes that nChanPar()=1 or nChanMat()
1344 :
1345 0 : for (Int ich=0; ich<nChanMat(); ich++) {
1346 :
1347 0 : oneMueller.reference(Miter.array());
1348 0 : oneMOK.reference(MOKiter.array());
1349 0 : onePar.reference(Piter.array());
1350 0 : onePOK.reference(POKiter.array());
1351 :
1352 : // TBD What if calcOneMueller needs freq value info?
1353 :
1354 0 : calcOneMueller(oneMueller,oneMOK,onePar,onePOK);
1355 :
1356 : // Advance iterators, as required
1357 0 : Miter.next();
1358 0 : MOKiter.next();
1359 0 : if (freqDepPar()) {
1360 0 : Piter.next();
1361 0 : POKiter.next();
1362 : }
1363 :
1364 : }
1365 :
1366 : // Step to next baseline's pars if we didn't in channel loop
1367 0 : if (!freqDepPar()) {
1368 0 : Piter.next();
1369 0 : POKiter.next();
1370 : }
1371 : }
1372 :
1373 0 : }
1374 :
1375 0 : void VisMueller::calcOneMueller(Vector<Complex>& /*mat*/, Vector<Bool>& /*mOk*/,
1376 : const Vector<Complex>& /*par*/, const Vector<Bool>& /*pOk*/) {
1377 :
1378 0 : if (prtlev()>10) cout << " VM::calcOneMueller()" << endl;
1379 :
1380 : // If Mueller matrix is trivial, shouldn't get here
1381 0 : if (trivialMuellerElem())
1382 0 : throw(AipsError("Trivial Mueller Matrix logic error."));
1383 :
1384 : // Otherwise, this method apparently hasn't been specialized, as required
1385 : else
1386 0 : throw(AipsError("Unknown non-trivial Mueller-from-parameter calculation requested."));
1387 :
1388 : }
1389 :
1390 0 : void VisMueller::invMueller() {
1391 :
1392 : // This method only called in apply context!
1393 0 : AlwaysAssert((isApplied()),AipsError);
1394 :
1395 0 : if (prtlev()>6) cout << " VM::invMueller()" << endl;
1396 :
1397 0 : M().sync(currMElem()(0,0,0),currMElemOK()(0,0,0));
1398 0 : for (Int ibln=0;ibln<nCalMat();ibln++)
1399 0 : for (Int ichan=0; ichan<nChanMat(); ++ichan, M()++)
1400 : // If matrix elements look ok so far, attempt to invert
1401 : // (if invert fails, matrix is zeroed and meOk is set accordingly)
1402 0 : M().invert();
1403 :
1404 0 : }
1405 :
1406 0 : void VisMueller::setMatByOk() {
1407 :
1408 : // This method only called in apply context!
1409 0 : AlwaysAssert((isApplied()),AipsError);
1410 :
1411 0 : if (prtlev()>6)
1412 0 : cout << " VM::setMatByOk()" << endl;
1413 :
1414 0 : M().sync(currMElem()(0,0,0),currMElemOK()(0,0,0));
1415 0 : for (Int ibln=0;ibln<nCalMat();ibln++)
1416 0 : for (Int ichan=0; ichan<nChanMat(); ++ichan, M()++)
1417 : // If matrix elements look ok so far, attempt to invert
1418 : // (if invert fails, matrix is zeroed and meOk is set accordingly)
1419 0 : M().setMatByOk();
1420 :
1421 0 : }
1422 :
1423 0 : void VisMueller::createMueller() {
1424 :
1425 0 : if (prtlev()>6) cout << " VM::createMueller()" << endl;
1426 :
1427 : // Delete current Mueller if wrong type
1428 0 : Mueller::MuellerType mtype(this->muellerType());
1429 0 : if (M_[currSpw()] &&
1430 0 : M().type() != mtype)
1431 0 : delete M_[currSpw()];
1432 :
1433 : // If needed, construct the correct Mueller
1434 0 : if (!M_[currSpw()])
1435 0 : M_[currSpw()] = casa::createMueller(mtype);
1436 :
1437 :
1438 : // Nominal synchronization is with currMElem()(0,0,0);
1439 0 : M().sync(currMElem()(0,0,0),currMElemOK()(0,0,0));
1440 :
1441 0 : if (prtlev()>6) cout << " currMElem()(0,0,0) = " << currMElem()(0,0,0) << endl;
1442 0 : if (prtlev()>6) cout << " currMElem()(0,0,1) = " << currMElem()(0,0,1) << endl;
1443 :
1444 :
1445 0 : }
1446 :
1447 0 : void VisMueller::syncWtScale() {
1448 :
1449 : // Ensure proper size according to Mueller matrix type
1450 0 : switch (this->muellerType()) {
1451 0 : case Mueller::Scalar:
1452 : case Mueller::Diagonal: {
1453 0 : Int nWtCorr=muellerNPar(muellerType());
1454 0 : currWtScale().resize(nWtCorr,1,nBln());
1455 0 : break;
1456 : }
1457 0 : default: {
1458 : // Only diag and scalar versions can adjust weights
1459 : // cout<< "Turning off calWt()" << endl;
1460 0 : calWt()=false;
1461 0 : return;
1462 : break;
1463 : }
1464 : }
1465 :
1466 : // Calculate the weight scaling
1467 0 : calcWtScale();
1468 : }
1469 :
1470 :
1471 0 : void VisMueller::calcWtScale() {
1472 :
1473 : // Assumes currMElem hasn't yet been inverted
1474 :
1475 : // Insist on single channel operation
1476 0 : AlwaysAssert( (nChanMat()==1), AipsError );
1477 :
1478 0 : Cube<Float> cWS;
1479 0 : cWS.reference(currWtScale());
1480 :
1481 : // Simple pre-inversion square of currMElem
1482 0 : cWS=real(currMElem()*conj(currMElem()));
1483 0 : cWS(!currMElemOK())=1.0;
1484 :
1485 0 : }
1486 :
1487 0 : void VisMueller::updateWt(Vector<Float>& wt,const Int& a1,const Int& a2) {
1488 :
1489 : // Assumes single channel
1490 0 : Vector<Float> ws(currWtScale().xyPlane(blnidx(a1,a2)).column(0));
1491 :
1492 0 : switch (V().type()) {
1493 0 : case VisVector::One: {
1494 0 : wt(0)*=ws(0);
1495 : }
1496 0 : case VisVector::Two: {
1497 0 : for (uInt ip=0;ip<2;++ip)
1498 0 : wt(ip)*=ws(ip*3);
1499 : }
1500 : default: {
1501 0 : for (uInt ip=0;ip<wt.nelements();++ip)
1502 0 : wt(ip)*=ws(ip);
1503 : }
1504 : }
1505 :
1506 0 : }
1507 :
1508 :
1509 0 : void VisMueller::initVisMueller() {
1510 :
1511 0 : if (prtlev()>2) cout << " VM::initVisMueller()" << endl;
1512 :
1513 0 : for (Int ispw=0;ispw<nSpw(); ispw++) {
1514 0 : currMElem_[ispw] = new Cube<Complex>();
1515 0 : currMElemOK_[ispw] = new Cube<Bool>();
1516 : }
1517 0 : }
1518 :
1519 0 : void VisMueller::deleteVisMueller() {
1520 :
1521 0 : if (prtlev()>2) cout << " VM::deleteVisMueller()" << endl;
1522 :
1523 0 : for (Int ispw=0; ispw<nSpw(); ispw++) {
1524 0 : if (currMElem_[ispw]) delete currMElem_[ispw];
1525 0 : if (currMElemOK_[ispw]) delete currMElemOK_[ispw];
1526 0 : if (M_[ispw]) delete M_[ispw];
1527 : }
1528 0 : currMElem_=NULL;
1529 0 : currMElemOK_=NULL;
1530 0 : M_=NULL;
1531 0 : }
1532 :
1533 :
1534 :
1535 :
1536 : // **********************************************************
1537 : // VisJones Implementations
1538 : //
1539 :
1540 0 : VisJones::VisJones(VisSet& vs) :
1541 : VisCal(vs), VisMueller(vs),
1542 0 : J1_(vs.numberSpw(),NULL),
1543 0 : J2_(vs.numberSpw(),NULL),
1544 0 : currJElem_(vs.numberSpw(),NULL),
1545 0 : currJElemOK_(vs.numberSpw(),NULL),
1546 0 : JValid_(vs.numberSpw(),false)
1547 : {
1548 0 : if (prtlev()>2) cout << "VJ::VJ(vs)" << endl;
1549 :
1550 0 : initVisJones();
1551 0 : }
1552 :
1553 0 : VisJones::VisJones(String msname,Int MSnAnt,Int MSnSpw) :
1554 : VisCal(msname,MSnAnt,MSnSpw),
1555 : VisMueller(msname,MSnAnt,MSnSpw),
1556 : J1_(MSnSpw,NULL),
1557 : J2_(MSnSpw,NULL),
1558 : currJElem_(MSnSpw,NULL),
1559 : currJElemOK_(MSnSpw,NULL),
1560 0 : JValid_(MSnSpw,false)
1561 : {
1562 0 : if (prtlev()>2) cout << "VJ::VJ(msname,MSnAnt,MSnSpw)" << endl;
1563 :
1564 0 : initVisJones();
1565 0 : }
1566 :
1567 :
1568 0 : VisJones::VisJones(const MSMetaInfoForCal& msmc) :
1569 : VisCal(msmc),
1570 : VisMueller(msmc),
1571 0 : J1_(nSpw(),NULL),
1572 0 : J2_(nSpw(),NULL),
1573 0 : currJElem_(nSpw(),NULL),
1574 0 : currJElemOK_(nSpw(),NULL),
1575 0 : JValid_(nSpw(),False)
1576 : {
1577 0 : if (prtlev()>2) cout << "VJ::VJ(MSMetaInfoForCal)" << endl;
1578 :
1579 0 : initVisJones();
1580 0 : }
1581 :
1582 0 : VisJones::VisJones(const Int& nAnt) :
1583 : VisCal(nAnt),
1584 : VisMueller(nAnt),
1585 : J1_(1,NULL),
1586 : J2_(1,NULL),
1587 : currJElem_(1,NULL),
1588 : currJElemOK_(1,NULL),
1589 0 : JValid_(1,false)
1590 : {
1591 0 : if (prtlev()>2) cout << "VJ::VJ(i,j,k)" << endl;
1592 :
1593 0 : initVisJones();
1594 0 : }
1595 :
1596 0 : VisJones::~VisJones() {
1597 0 : if (prtlev()>2) cout << "VJ::~VJ()" << endl;
1598 :
1599 0 : deleteVisJones();
1600 0 : }
1601 :
1602 0 : Mueller::MuellerType VisJones::muellerType() {
1603 :
1604 : // Ask Mueller to give us the answer:
1605 0 : return casa::muellerType(jonesType(),V().type());
1606 :
1607 : }
1608 :
1609 0 : void VisJones::state() {
1610 :
1611 : // Call parent
1612 0 : VisMueller::state();
1613 :
1614 0 : if (applyByJones()) {
1615 0 : if (prtlev()>3) cout << "VJ::state()" << endl;
1616 0 : cout << boolalpha;
1617 : // cout << "This is an antenna-based (Jones) calibration." << endl;
1618 0 : cout << " applyByJones=" << applyByJones() << endl;
1619 0 : cout << " jonesType= " << jonesType() << endl;
1620 0 : cout << " trivialJonesElem= " << trivialJonesElem() << endl;
1621 0 : cout << " JValid() = " << JValid() << endl;
1622 :
1623 0 : cout << " currJElem().shape() = " << currJElem().shape()
1624 0 : << " (" << currJElem().data() << ")" << endl;
1625 0 : cout << " currJElemOK().shape() = " << currJElemOK().shape()
1626 0 : << " (" << currJElemOK().data() << ") "
1627 0 : << " (ntrue=" << ntrue(currJElemOK()) << ")" << endl;
1628 :
1629 0 : cout << "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -" << endl;
1630 : }
1631 0 : }
1632 :
1633 : // VisJones: PROTECTED methods
1634 :
1635 :
1636 : // Apply this calibration to VisBuffer visibilities
1637 0 : void VisJones::applyCal(VisBuffer& vb, Cube<Complex>& Vout,
1638 : Bool trial) {
1639 :
1640 0 : if (prtlev()>3) cout << " VJ::applyCal()" << endl;
1641 :
1642 : // CURRENTLY ASSUMES ONLY *ONE* TIMESTAMP IN VISBUFFER
1643 :
1644 : /*
1645 : cout << "VJ::applyCal: type= " << typeName() << " trial = "
1646 : << boolalpha << trial << " calWt = " << calWt()
1647 : << " freqDepPar() = " << freqDepPar()
1648 : << " freqDepMat() = " << freqDepMat()
1649 : << endl;
1650 : */
1651 :
1652 0 : if (applyByMueller())
1653 0 : VisMueller::applyCal(vb,Vout);
1654 : else {
1655 :
1656 : // TBD: applyByJones()=true necessarily
1657 :
1658 : // Data info/indices
1659 : Int* dataChan;
1660 0 : Bool* flagR=&vb.flagRow()(0);
1661 0 : Bool* flag=&vb.flag()(0,0);
1662 0 : Int* a1=&vb.antenna1()(0);
1663 0 : Int* a2=&vb.antenna2()(0);
1664 0 : Matrix<Float> wtmat;
1665 :
1666 : // Prepare to cal weights
1667 0 : if (!trial)
1668 0 : wtmat.reference(vb.weightMat());
1669 :
1670 0 : ArrayIterator<Float> wt(wtmat,1);
1671 0 : Vector<Float> wtvec;
1672 :
1673 : // Alert Jones matrices to whether data is scalar or not
1674 : // (this is relevant only for proper handling of flags
1675 : // in case of scalar data, for now)
1676 0 : if (V().type()==VisVector::One) {
1677 0 : J1().setScalarData(true);
1678 0 : J2().setScalarData(true);
1679 : }
1680 : else {
1681 0 : J1().setScalarData(false);
1682 0 : J2().setScalarData(false);
1683 : }
1684 :
1685 : // iterate rows
1686 0 : Int& nRow(vb.nRow());
1687 0 : Int& nChanDat(vb.nChannel());
1688 : // cout << currSpw() << " startChan() = " << startChan() << " nChanMat() = " << nChanMat() << " nChanDat="<<nChanDat <<endl;
1689 0 : for (Int row=0; row<nRow; row++,flagR++,a1++,a2++) {
1690 :
1691 : // Solution channel registration
1692 0 : Int solCh0(0);
1693 0 : dataChan=&vb.channel()(0);
1694 :
1695 : // If cal _parameters_ are not freqDep (e.g., a delay)
1696 : // the startChan() should be the same as the first data channel
1697 0 : if (freqDepMat() && !freqDepPar())
1698 0 : startChan()=(*dataChan);
1699 :
1700 : // Solution and data array registration
1701 0 : J1().sync(currJElem()(0,solCh0,*a1),currJElemOK()(0,solCh0,*a1));
1702 0 : J2().sync(currJElem()(0,solCh0,*a2),currJElemOK()(0,solCh0,*a2));
1703 0 : if (!trial)
1704 0 : V().sync(Vout(0,0,row));
1705 :
1706 :
1707 0 : for (Int chn=0; chn<nChanDat; chn++,flag++,V()++,dataChan++) {
1708 :
1709 0 : if (trial) {
1710 : // only update flag info
1711 0 : J1().applyFlag(*flag);
1712 0 : J2().applyFlag(*flag);
1713 : }
1714 : else {
1715 : // if this data channel unflagged
1716 0 : J1().applyRight(V(),*flag);
1717 0 : J2().applyLeft(V(),*flag);
1718 : }
1719 :
1720 : // inc soln ch axis if freq-dependent (and next dataChan within soln)
1721 0 : if (freqDepMat()) {
1722 0 : J1()++;
1723 0 : J2()++;
1724 : }
1725 :
1726 : } // chn
1727 :
1728 :
1729 : // If requested, update the weights
1730 0 : if (!trial && calWt()) {
1731 0 : wtvec.reference(wt.array());
1732 0 : updateWt(wtvec,*a1,*a2);
1733 : }
1734 :
1735 0 : if (!trial)
1736 0 : wt.next();
1737 :
1738 : }
1739 :
1740 : }
1741 :
1742 0 : }
1743 :
1744 : // Apply this calibration to VisBuffer visibilities
1745 0 : void VisJones::applyCal2(vi::VisBuffer2& vb,
1746 : Cube<Complex>& Vout, Cube<Float>& Wout,
1747 : Bool trial) {
1748 :
1749 0 : if (prtlev()>3) cout << " VJ::applyCal()" << endl;
1750 :
1751 : // CURRENTLY ASSUMES ONLY *ONE* TIMESTAMP IN VISBUFFER
1752 :
1753 : /*
1754 : cout << "VJ::applyCal: type= " << typeName() << " trial = "
1755 : << boolalpha << trial << " calWt = " << calWt()
1756 : << " freqDepPar() = " << freqDepPar()
1757 : << " freqDepMat() = " << freqDepMat()
1758 : << endl;
1759 : */
1760 :
1761 0 : if (applyByMueller())
1762 0 : throw(AipsError("applyByMueller call from VisJones::applyCal2 NYI."));
1763 : // VisMueller::applyCal(vb,Vout);
1764 : else {
1765 :
1766 : // TBD: applyByJones()=true necessarily
1767 :
1768 : // References to VB2's contents' _data_
1769 0 : Vector<Bool> flagRv(vb.flagRow());
1770 0 : Vector<Int> a1v(vb.antenna1());
1771 0 : Vector<Int> a2v(vb.antenna2());
1772 0 : Cube<Bool> flagCube(vb.flagCube());
1773 0 : Cube<Complex> visCube(Vout);
1774 0 : ArrayIterator<Float> wt(Wout,2);
1775 0 : Matrix<Float> wtmat;
1776 :
1777 : // Data info/indices
1778 : Int* dataChan;
1779 0 : Bool* flagR=&flagRv(0);
1780 0 : Int* a1=&a1v(0);
1781 0 : Int* a2=&a2v(0);
1782 :
1783 : // Ensure VisVector for data acces has correct form
1784 0 : Int ncorr(vb.nCorrelations());
1785 0 : if (V().type() != ncorr)
1786 0 : V().setType(visType(ncorr));
1787 :
1788 :
1789 : // Alert Jones matrices to whether data is scalar or not
1790 : // (this is relevant only for proper handling of flags
1791 : // in case of scalar data, for now)
1792 0 : if (V().type()==VisVector::One) {
1793 0 : J1().setScalarData(true);
1794 0 : J2().setScalarData(true);
1795 : }
1796 : else {
1797 0 : J1().setScalarData(false);
1798 0 : J2().setScalarData(false);
1799 : }
1800 :
1801 : // iterate rows
1802 0 : Int nRow=vb.nRows();
1803 0 : Int nChanDat=vb.nChannels();
1804 0 : Vector<Int> dataChanv(vb.getChannelNumbers(0)); // All rows have same chans
1805 : // cout << currSpw() << " startChan() = " << startChan() << " nChanMat() = " << nChanMat() << " nChanDat="<<nChanDat <<endl;
1806 0 : for (Int row=0; row<nRow; row++,flagR++,a1++,a2++) {
1807 :
1808 : // Solution channel registration
1809 0 : Int solCh0(0);
1810 0 : dataChan=&dataChanv(0);
1811 :
1812 : // If cal _parameters_ are not freqDep (e.g., a delay)
1813 : // the startChan() should be the same as the first data channel
1814 0 : if (freqDepMat() && !freqDepPar())
1815 0 : startChan()=(*dataChan);
1816 :
1817 : // Solution and data array registration
1818 0 : J1().sync(currJElem()(0,solCh0,*a1),currJElemOK()(0,solCh0,*a1));
1819 0 : J2().sync(currJElem()(0,solCh0,*a2),currJElemOK()(0,solCh0,*a2));
1820 0 : V().sync(visCube(0,0,row),flagCube(0,0,row));
1821 :
1822 0 : for (Int chn=0; chn<nChanDat; chn++,V()++,dataChan++) {
1823 :
1824 0 : if (trial) {
1825 : // only update flag info
1826 0 : J1().flagRight(V());
1827 0 : J2().flagLeft(V());
1828 : }
1829 : else {
1830 : // if this data channel unflagged
1831 0 : J1().applyRight(V());
1832 0 : J2().applyLeft(V());
1833 : }
1834 :
1835 : // inc soln ch axis if freq-dependent (and next dataChan within soln)
1836 0 : if (freqDepMat()) {
1837 0 : J1()++;
1838 0 : J2()++;
1839 : }
1840 :
1841 : } // chn
1842 :
1843 :
1844 : // If requested, update the weights
1845 0 : if (!trial && calWt()) {
1846 0 : wtmat.reference(wt.array());
1847 0 : updateWt2(wtmat,*a1,*a2);
1848 : }
1849 :
1850 0 : if (!trial)
1851 0 : wt.next();
1852 :
1853 : }
1854 :
1855 : }
1856 :
1857 0 : }
1858 :
1859 :
1860 0 : void VisJones::syncCalMat(const Bool& doInv) {
1861 :
1862 0 : if (prtlev()>5) cout << " VJ::syncCalMat()"
1863 0 : << " (JValid()=" << JValid() << ")" << endl;
1864 :
1865 : // cout << " VisCal::syncCalMat(doInv): " << currCPar().data() <<" "<< currJElem().data()
1866 : // << endl;
1867 :
1868 : // If necessary, synchronize the Jones matrices
1869 0 : if (!JValid()) syncJones(doInv);
1870 :
1871 : // cout << " VisCal::syncCalMat(doInv): " << currCPar().data() <<" "<< currJElem().data()
1872 : // << endl;
1873 :
1874 : // Be sure sync'd matrices at their origin
1875 0 : J1().origin();
1876 0 : J2().origin();
1877 :
1878 : // If requested and necessary, synchronize the Mueller matrices
1879 : // (NEVER invert Muellers, as Jones already have been)
1880 0 : if (applyByMueller() && !MValid()) syncMueller(false);
1881 :
1882 0 : }
1883 :
1884 :
1885 0 : void VisJones::syncJones(const Bool& doInv) {
1886 :
1887 0 : if (prtlev()>6) cout << " VJ::syncJones()" << endl;
1888 :
1889 : // If Jones pars are just the matrix elements:
1890 0 : if (trivialJonesElem()) {
1891 :
1892 : // Matrix Elem cache references par cache
1893 0 : currJElem().reference(currCPar());
1894 0 : currJElemOK().reference(currParOK());
1895 :
1896 : // EXCEPT, ensure uniqueness if taking the inverse
1897 : // (this makes a copy, can we avoid?)
1898 0 : if (doInv) {
1899 : //cout << " Unique: " << currJElem().data() << "-->";
1900 0 : currJElem().unique();
1901 : //cout << currJElem().data() << endl;
1902 : }
1903 : }
1904 : else {
1905 :
1906 : // Make local matrix element cache the correct size:
1907 0 : currJElem().resize(jonesNPar(jonesType()),nChanMat(),nAnt());
1908 0 : currJElem().unique(); // Ensure uniqueness!
1909 :
1910 : // OK matches size of the J matrix itself
1911 0 : currJElemOK().resize(jonesNPar(jonesType()),nChanMat(),nAnt());
1912 0 : currJElem().unique(); // Ensure uniqueness!
1913 0 : currJElem()=false;
1914 :
1915 : // The matrix state is invalid until we actually calculate them
1916 0 : invalidateJ();
1917 :
1918 : // And calculate the matrix elements for all baselines
1919 0 : calcAllJones();
1920 :
1921 : }
1922 :
1923 : /*
1924 : if (type()==VisCal::B) {
1925 : cout << typeName() << ": currJElem().shape() = " << currJElem().shape() << endl;
1926 : cout << "currJElem() = " << currJElem() << endl;
1927 : cout << "currJElemOK() = " << currJElemOK() << endl;
1928 : cout << "freqDepMat() = " << boolalpha << freqDepMat() << endl;
1929 : }
1930 : */
1931 :
1932 : // Pre-inv syncWtScale:
1933 0 : if (calWt()) syncWtScale();
1934 :
1935 : // Ensure Jones Matrix renders are ok
1936 0 : this->createJones();
1937 :
1938 : // Invert, if requested
1939 0 : if (doInv) invJones();
1940 :
1941 : // Set matrix elements according to OK flags
1942 0 : setMatByOk();
1943 :
1944 : // Jones matrices now valid
1945 0 : validateJ();
1946 0 : invalidateM(); // still invalid
1947 :
1948 0 : }
1949 :
1950 0 : void VisJones::calcAllJones() {
1951 :
1952 0 : if (prtlev()>6) cout << " VJ::calcAllJones()" << endl;
1953 :
1954 : // Should handle OK flags in this method, and only
1955 : // do Jones calc if OK
1956 :
1957 0 : Vector<Complex> oneJones;
1958 0 : Vector<Bool> oneJOK;
1959 0 : Vector<Complex> onePar;
1960 0 : Vector<Bool> onePOK;
1961 :
1962 0 : ArrayIterator<Complex> Jiter(currJElem(),1);
1963 0 : ArrayIterator<Bool> JOKiter(currJElemOK(),1);
1964 0 : ArrayIterator<Complex> Piter(currCPar(),1);
1965 0 : ArrayIterator<Bool> POKiter(currParOK(),1);
1966 :
1967 0 : for (Int iant=0; iant<nAnt(); iant++) {
1968 :
1969 : // The following assumes that nChanPar()=1 or nChanMat()
1970 :
1971 0 : for (Int ich=0; ich<nChanMat(); ich++) {
1972 :
1973 0 : oneJones.reference(Jiter.array());
1974 0 : oneJOK.reference(JOKiter.array());
1975 0 : onePar.reference(Piter.array());
1976 0 : onePOK.reference(POKiter.array());
1977 :
1978 : // Calculate the Jones matrix
1979 0 : calcOneJones(oneJones,oneJOK,onePar,onePOK);
1980 :
1981 : // Advance iterators
1982 0 : Jiter.next();
1983 0 : JOKiter.next();
1984 0 : if (freqDepPar()) {
1985 0 : Piter.next();
1986 0 : POKiter.next();
1987 : }
1988 :
1989 : }
1990 : // Step to next antenns's pars if we didn't in channel loop
1991 0 : if (!freqDepPar()) {
1992 0 : Piter.next();
1993 0 : POKiter.next();
1994 : }
1995 : }
1996 0 : }
1997 :
1998 0 : void VisJones::calcOneJones(Vector<Complex>& /*mat*/, Vector<Bool>& /*mOk*/,
1999 : const Vector<Complex>& /*par*/, const Vector<Bool>& /*pOk*/ ) {
2000 :
2001 0 : if (prtlev()>10) cout << " VJ::calcOneJones()" << endl;
2002 :
2003 : // If Jones matrices are trivial, should never reach here!
2004 0 : if (trivialJonesElem())
2005 0 : throw(AipsError("Trivial Jones Matrix logic error."));
2006 :
2007 : // Otherwise, this method apparently hasn't been specialized, as required
2008 : else
2009 0 : throw(AipsError("Unknown non-trivial Jones-from-parameter calculation requested."));
2010 :
2011 : }
2012 :
2013 0 : void VisJones::invJones() {
2014 :
2015 0 : if (prtlev()>6) cout << " VJ::invJones()" << endl;
2016 :
2017 0 : J1().sync(currJElem()(0,0,0),currJElemOK()(0,0,0));
2018 0 : for (Int iant=0;iant<nAnt();iant++)
2019 0 : for (Int ichan=0; ichan<nChanMat();++ichan,J1()++)
2020 : // If matrix elements look ok so far, attempt to invert
2021 : // (if invert fails, currJElemOK will be set accordingly)
2022 0 : J1().invert();
2023 :
2024 0 : }
2025 :
2026 0 : void VisJones::setMatByOk() {
2027 :
2028 0 : if (prtlev()>6)
2029 0 : cout << " VJ::setMatByOk" << endl;
2030 :
2031 0 : J1().sync(currJElem()(0,0,0),currJElemOK()(0,0,0));
2032 0 : for (Int iant=0;iant<nAnt();iant++)
2033 0 : for (Int ichan=0; ichan<nChanMat();++ichan,J1()++)
2034 : // If matrix elements look ok so far, attempt to invert
2035 : // (if invert fails, currJElemOK will be set accordingly)
2036 0 : J1().setMatByOk();
2037 :
2038 0 : }
2039 :
2040 0 : void VisJones::calcAllMueller() {
2041 :
2042 0 : if (prtlev()>6) cout << " VJ::calcAllMueller()" << endl;
2043 :
2044 0 : M().sync(currMElem()(0,0,0));
2045 0 : for (Int a1=0;a1<nAnt();++a1) {
2046 0 : J1().sync(currJElem()(0,0,a1),currJElemOK()(0,0,a1));
2047 0 : for (Int a2=a1;a2<nAnt();++a2) {
2048 0 : J2().sync(currJElem()(0,0,a2),currJElemOK()(0,0,a2));
2049 0 : for (Int ich=0;ich<nChanMat();ich++,J1()++,J2()++,M()++)
2050 0 : M().fromJones(J1(),J2());
2051 : }
2052 : }
2053 :
2054 0 : }
2055 :
2056 0 : void VisJones::createJones() {
2057 :
2058 0 : if (prtlev()>6) cout << " VJ::createJones()" << endl;
2059 :
2060 : // Delete current Jones if wrong type
2061 0 : Jones::JonesType jtype(jonesType());
2062 :
2063 0 : if (J1_[currSpw()] &&
2064 0 : J1().type() != jtype)
2065 0 : delete J1_[currSpw()];
2066 :
2067 0 : if (J2_[currSpw()] &&
2068 0 : J2().type() != jtype)
2069 0 : delete J2_[currSpw()];
2070 :
2071 : // If needed, construct the correct Jones
2072 0 : if (!J1_[currSpw()])
2073 0 : J1_[currSpw()] = casa::createJones(jtype);
2074 :
2075 0 : if (!J2_[currSpw()])
2076 0 : J2_[currSpw()] = casa::createJones(jtype);
2077 :
2078 :
2079 : // Nominal synchronization is with currJElem()(0,0,0):
2080 0 : J1().sync(currJElem()(0,0,0),currJElemOK()(0,0,0));
2081 0 : J2().sync(currJElem()(0,0,0),currJElemOK()(0,0,0));
2082 :
2083 0 : }
2084 :
2085 0 : void VisJones::syncWtScale() {
2086 :
2087 : // cout << "VJ::syncWtScale (" << typeName() << ")" << endl;
2088 :
2089 :
2090 : // Ensure proper size according to Jones matrix type
2091 0 : switch (this->jonesType()) {
2092 0 : case Jones::Scalar:
2093 : case Jones::Diagonal: {
2094 0 : Int nWtScale=jonesNPar(jonesType());
2095 0 : currWtScale().resize(nWtScale,1,nAnt());
2096 0 : break;
2097 : }
2098 0 : default: {
2099 : // Only diag and scalar versions can adjust weights
2100 : // cout<< "Turning off calWt()" << endl;
2101 0 : calWt()=false;
2102 0 : return;
2103 : break;
2104 : }
2105 : }
2106 :
2107 : // Calculate the weight scale factors (specializable)
2108 0 : calcWtScale();
2109 :
2110 :
2111 : // cout << "VJ::syncWtScale: currWtScale() = " << currWtScale() << endl;
2112 :
2113 : }
2114 :
2115 0 : void VisJones::calcWtScale() {
2116 :
2117 : // Assumes currJElem has not yet been inverted!
2118 :
2119 : // Insist on single channel operation
2120 0 : AlwaysAssert( (nChanMat()==1), AipsError );
2121 :
2122 : // Just a reference
2123 0 : Cube<Float> cWS(currWtScale());
2124 :
2125 : // We use simple (pre-inversion) square of currJElem
2126 0 : cWS=real(currJElem()*conj(currJElem()));
2127 0 : cWS(!currJElemOK())=1.0;
2128 :
2129 0 : }
2130 :
2131 0 : void VisJones::updateWt(Vector<Float>& wt,const Int& a1,const Int& a2) {
2132 :
2133 0 : Vector<Float> ws1(currWtScale().xyPlane(a1).column(0));
2134 0 : Vector<Float> ws2(currWtScale().xyPlane(a2).column(0));
2135 :
2136 : /*
2137 : if (a1==0 && a2==1) {
2138 : cout << "jonestype() = " << jonesType()
2139 : << " jonesNPar(jonesType()) = " << jonesNPar(jonesType())
2140 : << " nPar() = " << nPar()
2141 : << endl;
2142 : cout << currSpw() << " "
2143 : << a1 << " "
2144 : << a2 << " "
2145 : << wt << " "
2146 : << ws1 << " "
2147 : << ws2 << " "
2148 : << ws1(0/nPar()) << " "
2149 : << ws2(0%nPar()) << " ";
2150 : }
2151 : */
2152 :
2153 0 : switch(jonesType()) {
2154 0 : case Jones::Scalar: {
2155 : // pol-indep corrections very simple; all correlations
2156 : // corrected by same value
2157 0 : Float ws=(ws1(0)*ws2(0));
2158 0 : wt*=ws;
2159 0 : break;
2160 : }
2161 0 : case Jones::Diagonal: {
2162 0 : switch (V().type()) {
2163 0 : case VisVector::Two: {
2164 0 : for (uInt ip=0;ip<2;++ip)
2165 0 : wt(ip)*=(ws1(ip)*ws2(ip));
2166 0 : break;
2167 : }
2168 0 : default: {
2169 : // NB: This always will apply the first weight scale info the a single corr
2170 0 : for (uInt ip=0;ip<wt.nelements();++ip) {
2171 0 : wt(ip)*=ws1(ip/2);
2172 0 : wt(ip)*=ws2(ip%2);
2173 : }
2174 0 : break;
2175 : }
2176 : }
2177 0 : break;
2178 : }
2179 0 : default:
2180 : // We don't calibrate weights for General Jones matrices (yet)
2181 0 : break;
2182 : }
2183 :
2184 : /*
2185 : if (a1==0 && a2==1)
2186 : cout << " ---> " << wt << endl;
2187 : */
2188 0 : }
2189 :
2190 : // WEIGHT_SPECTRUM-capable version
2191 0 : void VisJones::updateWt2(Matrix<Float>& wt,const Int& a1,const Int& a2) {
2192 :
2193 0 : Int nVco=wt.shape()[0]; // ==nChanDat()?
2194 0 : Int nVch=wt.shape()[1]; // ==nCorrDat()?
2195 :
2196 0 : Matrix<Float> ws1(currWtScale().xyPlane(a1));
2197 0 : Matrix<Float> ws2(currWtScale().xyPlane(a2));
2198 :
2199 0 : Int nCch=ws1.shape()[1]; // ==nChanMat()?
2200 :
2201 : // Channel counts either idential, or single-chan calibration
2202 : // alwaysAssert( nVch==nCch || nCch==1, AipsError);
2203 :
2204 0 : switch(jonesType()) {
2205 0 : case Jones::Scalar: {
2206 : // pol-indep corrections very simple; all correlations
2207 : // corrected by same value, per channel
2208 0 : for (Int ich=0;ich<nVch;++ich) {
2209 0 : Int iCch=ich%nCch; // 0 or ich
2210 0 : Float ws=(ws1(0,iCch)*ws2(0,iCch));
2211 0 : Matrix<Float> wtm(wt(Slice(),Slice(ich,1,1)));
2212 0 : wtm*=ws;
2213 : }
2214 0 : break;
2215 : }
2216 0 : case Jones::Diagonal: {
2217 :
2218 0 : Int nCorrPerPol=max(nVco/2,1); // >=1
2219 0 : for (Int ich=0;ich<nVch;++ich) {
2220 0 : Int iCch=ich%nCch; // 0 or ich
2221 0 : for (Int ico=0;ico<nVco;++ico)
2222 0 : wt(ico,ich) *= ( ws1(ico/nCorrPerPol,iCch) *
2223 0 : ws2(ico%2,iCch) );
2224 : }
2225 0 : break;
2226 : }
2227 0 : default:
2228 : // We don't calibrate weights for General Jones matrices (yet)
2229 0 : break;
2230 : }
2231 :
2232 0 : }
2233 :
2234 0 : void VisJones::initVisJones() {
2235 :
2236 0 : if (prtlev()>2) cout << " VJ::initVisJones()" << endl;
2237 :
2238 0 : for (Int ispw=0;ispw<nSpw(); ispw++) {
2239 0 : currJElem_[ispw] = new Cube<Complex>();
2240 0 : currJElemOK_[ispw] = new Cube<Bool>();
2241 : }
2242 0 : }
2243 :
2244 0 : void VisJones::deleteVisJones() {
2245 :
2246 0 : if (prtlev()>2) cout << " VJ::deleteVisJones()" << endl;
2247 :
2248 0 : for (Int ispw=0; ispw<nSpw(); ispw++) {
2249 0 : if (currJElem_[ispw]) delete currJElem_[ispw];
2250 0 : if (currJElemOK_[ispw]) delete currJElemOK_[ispw];
2251 0 : if (J1_[ispw]) delete J1_[ispw];
2252 0 : if (J2_[ispw]) delete J2_[ispw];
2253 :
2254 : }
2255 0 : currJElem_=NULL;
2256 0 : currJElemOK_=NULL;
2257 0 : J1_=NULL;
2258 0 : J2_=NULL;
2259 0 : }
2260 :
2261 : } //# NAMESPACE CASA - END
2262 :
|
