Line data Source code
1 : //# StandardVisCal.cc: Implementation of Standard VisCal types
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003,2011
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 :
27 : #include <synthesis/MeasurementComponents/StandardVisCal.h>
28 : #include <synthesis/MeasurementComponents/CalCorruptor.h>
29 :
30 : #include <msvis/MSVis/VisBuffer.h>
31 : #include <msvis/MSVis/VisBuffAccumulator.h>
32 : #include <synthesis/CalTables/CTIter.h>
33 : #include <casacore/ms/MeasurementSets/MSColumns.h>
34 : #include <synthesis/MeasurementEquations/VisEquation.h>
35 : #include <synthesis/MeasurementComponents/SolveDataBuffer.h>
36 : #include <casacore/scimath/Fitting/LSQFit.h>
37 : #include <casacore/scimath/Fitting/LinearFit.h>
38 : #include <casacore/scimath/Functionals/CompiledFunction.h>
39 : #include <casacore/scimath/Functionals/Polynomial.h>
40 : #include <casacore/scimath/Mathematics/AutoDiff.h>
41 : #include <casacore/casa/BasicMath/Math.h>
42 : #include <casacore/tables/TaQL/ExprNode.h>
43 :
44 : #include <casacore/casa/Arrays/ArrayMath.h>
45 : #include <casacore/casa/Arrays/MatrixMath.h>
46 : #include <casacore/casa/BasicSL/String.h>
47 : #include <casacore/casa/Utilities/Assert.h>
48 : #include <casacore/casa/Utilities/GenSort.h>
49 : #include <casacore/casa/Exceptions/Error.h>
50 : #include <casacore/casa/OS/Memory.h>
51 : #include <casacore/casa/System/Aipsrc.h>
52 :
53 : #include <sstream>
54 :
55 : #include <casacore/measures/Measures/MCBaseline.h>
56 : #include <casacore/measures/Measures/MDirection.h>
57 : #include <casacore/measures/Measures/MEpoch.h>
58 : #include <casacore/measures/Measures/MeasTable.h>
59 :
60 : #include <casacore/casa/Logging/LogMessage.h>
61 : #include <casacore/casa/Logging/LogSink.h>
62 : // math.h ?
63 :
64 : using namespace casacore;
65 : namespace casa { //# NAMESPACE CASA - BEGIN
66 :
67 :
68 : // **********************************************************
69 : // PJones
70 : //
71 :
72 0 : PJones::PJones(VisSet& vs) :
73 : VisCal(vs),
74 : VisMueller(vs),
75 : VisJones(vs),
76 : pjonestype_(Jones::Diagonal),
77 0 : pa_()
78 : {
79 0 : if (prtlev()>2) cout << "P::P(vs)" << endl;
80 0 : }
81 :
82 0 : PJones::PJones(String msname,Int MSnAnt,Int MSnSpw) :
83 : VisCal(msname,MSnAnt,MSnSpw),
84 : VisMueller(msname,MSnAnt,MSnSpw),
85 : VisJones(msname,MSnAnt,MSnSpw),
86 : pjonestype_(Jones::Diagonal),
87 0 : pa_()
88 : {
89 0 : if (prtlev()>2) cout << "P::P(msname,MSnAnt,MSnSpw)" << endl;
90 0 : }
91 :
92 0 : PJones::PJones(const MSMetaInfoForCal& msmc) :
93 : VisCal(msmc),
94 : VisMueller(msmc),
95 : VisJones(msmc),
96 : pjonestype_(Jones::Diagonal),
97 0 : pa_()
98 : {
99 0 : if (prtlev()>2) cout << "P::P(msmc)" << endl;
100 0 : }
101 :
102 :
103 0 : PJones::~PJones() {
104 0 : if (prtlev()>2) cout << "P::~P()" << endl;
105 0 : }
106 :
107 : // PJones needs to know pol basis and feed pa
108 0 : void PJones::syncMeta(const VisBuffer& vb) {
109 :
110 : // Call parent (sets currTime())
111 0 : VisJones::syncMeta(vb);
112 :
113 : // Basis
114 0 : if (vb.corrType()(0)==5) // Circulars
115 0 : pjonestype_=Jones::Diagonal;
116 0 : else if (vb.corrType()(0)==9) // Linears
117 0 : pjonestype_=Jones::General;
118 :
119 : // Get parallactic angle from the vb:
120 0 : pa().resize(nAnt());
121 0 : pa() = vb.feed_pa(currTime());
122 :
123 0 : }
124 :
125 : // PJones needs to know pol basis and feed pa
126 0 : void PJones::syncMeta2(const vi::VisBuffer2& vb) {
127 :
128 : // Call parent (sets currTime())
129 0 : VisJones::syncMeta2(vb);
130 :
131 : // Basis
132 0 : if (vb.correlationTypes()(0)==5) // Circulars
133 0 : pjonestype_=Jones::Diagonal;
134 0 : else if (vb.correlationTypes()(0)==9) // Linears
135 0 : pjonestype_=Jones::General;
136 :
137 : // Get parallactic angle from the vb:
138 0 : pa().resize(nAnt());
139 0 : pa() = vb.feedPa(currTime());
140 :
141 0 : }
142 :
143 0 : void PJones::calcPar() {
144 :
145 0 : if (prtlev()>6) cout << " VC::calcPar()" << endl;
146 :
147 : // Initialize parameter arrays
148 0 : currCPar().resize(1,1,nAnt());
149 0 : currParOK().resize(1,1,nAnt());
150 0 : currParOK()=true;
151 :
152 : // Fill currCPar() with exp(i*pa)
153 0 : Float* a=pa().data();
154 0 : Complex* cp=currCPar().data();
155 0 : Double ang(0.0);
156 0 : for (Int iant=0;iant<nAnt();++iant,++a,++cp) {
157 0 : ang=Double(*a);
158 0 : (*cp) = Complex(cos(ang),sin(ang)); // as a complex number
159 : }
160 : // cout << "ang = " << ang << endl;
161 :
162 : // Pars now valid, matrices not
163 0 : validateP();
164 0 : invalidateJ();
165 :
166 0 : }
167 :
168 : // Calculate a single Jones matrix by some means from parameters
169 0 : void PJones::calcOneJones(Vector<Complex>& mat, Vector<Bool>& mOk,
170 : const Vector<Complex>& par, const Vector<Bool>& pOk ) {
171 :
172 0 : if (prtlev()>10) cout << " P::calcOneJones()" << endl;
173 :
174 :
175 0 : if (pOk(0)) {
176 :
177 0 : switch (jonesType()) {
178 : // Circular version:
179 0 : case Jones::Diagonal: {
180 0 : mat(0)=conj(par(0)); // exp(-ia)
181 0 : mat(1)=par(0); // exp(ia)
182 0 : mOk=true;
183 0 : break;
184 : }
185 : // Linear version:
186 0 : case Jones::General: {
187 0 : Float a=arg(par(0));
188 0 : mat(0)=mat(3)=cos(a);
189 0 : mat(1)=sin(a);
190 0 : mat(2)=-mat(1);
191 0 : mOk=true;
192 0 : break;
193 : }
194 0 : default:
195 0 : throw(AipsError("PJones doesn't know if it is Diag (Circ) or General (Lin)"));
196 : break;
197 :
198 : }
199 :
200 : }
201 0 : }
202 :
203 :
204 :
205 : // **********************************************************
206 : // TJones Implementations
207 : //
208 :
209 0 : TJones::TJones(VisSet& vs) :
210 : VisCal(vs), // virtual base
211 : VisMueller(vs), // virtual base
212 : SolvableVisJones(vs), // immediate parent
213 0 : tcorruptor_p(NULL)
214 : {
215 0 : if (prtlev()>2) cout << "T::T(vs)" << endl;
216 0 : }
217 :
218 0 : TJones::TJones(String msname,Int MSnAnt,Int MSnSpw) :
219 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
220 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
221 : SolvableVisJones(msname,MSnAnt,MSnSpw), // immediate parent
222 0 : tcorruptor_p(NULL)
223 : {
224 0 : if (prtlev()>2) cout << "T::T(msname,MSnAnt,MSnSpw)" << endl;
225 0 : }
226 :
227 0 : TJones::TJones(const MSMetaInfoForCal& msmc) :
228 : VisCal(msmc),
229 : VisMueller(msmc),
230 : SolvableVisJones(msmc),
231 0 : tcorruptor_p(NULL)
232 : {
233 0 : if (prtlev()>2) cout << "T::T(msmc)" << endl;
234 0 : }
235 :
236 :
237 0 : TJones::TJones(const Int& nAnt) :
238 : VisCal(nAnt),
239 : VisMueller(nAnt),
240 : SolvableVisJones(nAnt),
241 0 : tcorruptor_p(NULL)
242 : {
243 0 : if (prtlev()>2) cout << "T::T(nAnt)" << endl;
244 0 : }
245 :
246 0 : TJones::~TJones() {
247 0 : if (prtlev()>2) cout << "T::~T()" << endl;
248 0 : }
249 :
250 0 : void TJones::guessPar(VisBuffer& vb) {
251 :
252 0 : if (prtlev()>4) cout << " T::guessPar(vb)" << endl;
253 :
254 : // Assumes: 1. corrs in canonical order
255 : // 2. vb has 1 channel (has been freq-averaged)
256 :
257 :
258 : // Make an antenna-based guess at T
259 : // Correlation membership-dependence
260 : // nCorr = 1: use icorr=0
261 : // nCorr = 2: use icorr=[0,1]
262 : // nCorr = 4: use icorr=[0,3]
263 :
264 0 : Int nCorr(1);
265 0 : Int nDataCorr(vb.visCube().shape()(0));
266 0 : Vector<Int> corridx(1,0);
267 0 : if (nDataCorr==2) {
268 0 : nCorr=2;
269 0 : corridx.resize(nCorr);
270 0 : corridx(0)=0;
271 0 : corridx(1)=1;
272 : }
273 0 : else if (nDataCorr==4) {
274 0 : nCorr=2;
275 0 : corridx.resize(nCorr);
276 0 : corridx(0)=0;
277 0 : corridx(1)=3;
278 : }
279 :
280 : // Find out which ants are available
281 : // TBD: count nominal guessant rows, insist not much less than nAnt
282 0 : Vector<Bool> antok(nAnt(),false);
283 0 : Vector<Bool> rowok(vb.nRow(),false);
284 0 : for (Int irow=0;irow<vb.nRow();++irow) {
285 : // Is this row ok?
286 0 : rowok(irow)= (!vb.flagRow()(irow) &&
287 0 : (vb.antenna1()(irow)!=vb.antenna2()(irow)) &&
288 0 : nfalse(vb.flag().column(irow))> 0 );
289 0 : if (rowok(irow)) {
290 0 : antok(vb.antenna1()(irow))=true;
291 0 : antok(vb.antenna2()(irow))=true;
292 : }
293 : }
294 :
295 : // Assume refant is the target ant, for starters
296 0 : Int guessant(refant());
297 :
298 : // If no refant specified, or no data for refant
299 : // base first guess on first good ant
300 0 : if (guessant<0 || !antok(guessant)) {
301 0 : guessant=0;
302 0 : while (antok(guessant)<1) guessant++;
303 : }
304 :
305 0 : AlwaysAssert(guessant>-1,AipsError);
306 :
307 0 : Cube<Complex>& V(vb.visCube());
308 0 : Float amp(0.0),ampave(0.0);
309 0 : Int namp(0);
310 0 : solveCPar()=Complex(0.0);
311 0 : for (Int irow=1;irow<vb.nRow();++irow) { // why not row zero? copied from Calibrator
312 :
313 0 : if (rowok(irow)) {
314 0 : Int a1=vb.antenna1()(irow);
315 0 : Int a2=vb.antenna2()(irow);
316 :
317 : // If this row contains the guessant
318 0 : if (a1 == guessant || a2==guessant) {
319 :
320 0 : for (Int icorr=0;icorr<nCorr;icorr++) {
321 0 : Complex& Vi(V(corridx(icorr),0,irow));
322 0 : amp=abs(Vi);
323 0 : if (amp>0.0f) {
324 0 : if (a1 == guessant)
325 0 : solveCPar()(0,0,a2)+=(conj(Vi)/amp/Float(nCorr));
326 : // solveCPar()(0,0,a2)+=(conj(Vi)/Float(nCorr));
327 : else
328 0 : solveCPar()(0,0,a1)+=((Vi)/amp/Float(nCorr));
329 : // solveCPar()(0,0,a1)+=((Vi)/Float(nCorr));
330 :
331 0 : ampave+=amp;
332 0 : namp++;
333 : // cout << " " << abs(Vi) << " " << arg(Vi)*180.0/C::pi << endl;
334 : }
335 : }
336 : } // guessant in row
337 : } // rowok
338 : } // irow
339 :
340 : // cout << "Guess:" << endl
341 : // << "amp = " << amplitude(solveCPar())
342 : // << endl;
343 :
344 :
345 : // Scale them by the mean amplitude
346 :
347 0 : if (namp>0) {
348 0 : ampave/=Float(namp);
349 0 : ampave=sqrt(ampave);
350 : // solveCPar()*=Complex(ampave);
351 0 : solveCPar()/=Complex(ampave);
352 0 : solveCPar()(0,0,guessant)=Complex(ampave);
353 0 : solveCPar()(LogicalArray(amplitude(solveCPar())==0.0f)) = Complex(ampave);
354 : }
355 : else
356 0 : solveCPar()=Complex(0.3);
357 :
358 0 : solveParOK()=true;
359 :
360 : /*
361 : cout << "Guess:" << endl
362 : << "amp = " << amplitude(solveCPar())
363 : << "phase = " << phase(solveCPar())
364 : << endl;
365 : */
366 0 : }
367 :
368 :
369 0 : void TJones::guessPar(SDBList& sdbs,const Bool&) {
370 :
371 0 : if (prtlev()>4) cout << " T::guessPar(sdbs)" << endl;
372 :
373 : // Assumes: 1. corrs in canonical order
374 :
375 : // Use just the first SDB in the SDBList, for now
376 0 : SolveDataBuffer& sdb(sdbs(0));
377 :
378 0 : AlwaysAssert(sdb.nChannels()==1,AipsError);
379 :
380 : // Make an antenna-based guess at T
381 : // Correlation membership-dependence
382 : // nCorr = 1: use icorr=0
383 : // nCorr = 2: use icorr=[0,1]
384 : // nCorr = 4: use icorr=[0,3]
385 :
386 0 : Int nCorr(1);
387 0 : Int nDataCorr(sdb.nCorrelations());
388 0 : Vector<Int> corridx(1,0);
389 0 : if (nDataCorr==2) {
390 0 : nCorr=2;
391 0 : corridx.resize(nCorr);
392 0 : corridx(0)=0;
393 0 : corridx(1)=1;
394 : }
395 0 : else if (nDataCorr==4) {
396 0 : nCorr=2;
397 0 : corridx.resize(nCorr);
398 0 : corridx(0)=0;
399 0 : corridx(1)=3;
400 : }
401 :
402 : // Find out which ants are available
403 : // TBD: count nominal guessant rows, insist not much less than nAnt
404 0 : Int nRow=sdb.nRows();
405 0 : Vector<Bool> antok(nAnt(),False);
406 0 : Vector<Bool> rowok(nRow,False);
407 0 : for (Int irow=0;irow<nRow;++irow) {
408 0 : const Int& a1(sdb.antenna1()(irow));
409 0 : const Int& a2(sdb.antenna2()(irow));
410 : // Is this row ok?
411 0 : rowok(irow)= (!sdb.flagRow()(irow) &&
412 0 : (a1!=a2) &&
413 0 : nfalse(sdb.flagCube().xyPlane(irow))> 0 );
414 0 : if (rowok(irow)) {
415 0 : antok(a1)=True;
416 0 : antok(a2)=True;
417 : }
418 : }
419 :
420 : // Assume refant is the target ant, for starters
421 0 : Int guessant(refant());
422 :
423 : // If no refant specified, or no data for refant
424 : // base first guess on first good ant
425 0 : if (guessant<0 || !antok(guessant)) {
426 0 : guessant=0;
427 0 : while (antok(guessant)<1) guessant++;
428 : }
429 :
430 0 : AlwaysAssert(guessant>-1,AipsError);
431 :
432 0 : const Cube<Complex>& V(sdb.visCubeCorrected());
433 0 : Float amp(0.0),ampave(0.0);
434 0 : Int namp(0);
435 0 : solveCPar()=Complex(0.0);
436 0 : for (Int irow=1;irow<nRow;++irow) { // why not row zero? copied from Calibrator
437 :
438 0 : if (rowok(irow)) {
439 0 : const Int& a1(sdb.antenna1()(irow));
440 0 : const Int& a2(sdb.antenna2()(irow));
441 :
442 : // If this row contains the guessant
443 0 : if (a1 == guessant || a2==guessant) {
444 :
445 0 : for (Int icorr=0;icorr<nCorr;icorr++) {
446 0 : const Complex& Vi(V(corridx(icorr),0,irow));
447 0 : amp=abs(Vi);
448 0 : if (amp>0.0f) {
449 0 : if (a1 == guessant)
450 0 : solveCPar()(0,0,a2)+=(conj(Vi)/amp/Float(nCorr));
451 : // solveCPar()(0,0,a2)+=(conj(Vi)/Float(nCorr));
452 : else
453 0 : solveCPar()(0,0,a1)+=((Vi)/amp/Float(nCorr));
454 : // solveCPar()(0,0,a1)+=((Vi)/Float(nCorr));
455 :
456 0 : ampave+=amp;
457 0 : namp++;
458 : // cout << " " << abs(Vi) << " " << arg(Vi)*180.0/C::pi << endl;
459 : }
460 : }
461 : } // guessant in row
462 : } // rowok
463 : } // irow
464 :
465 : // cout << "Guess:" << endl
466 : // << "amp = " << amplitude(solveCPar())
467 : // << endl;
468 :
469 :
470 : // Scale them by the mean amplitude
471 :
472 0 : if (namp>0) {
473 0 : ampave/=Float(namp);
474 0 : ampave=sqrt(ampave);
475 : // solveCPar()*=Complex(ampave);
476 0 : solveCPar()/=Complex(ampave);
477 0 : solveCPar()(0,0,guessant)=Complex(ampave);
478 0 : solveCPar()(LogicalArray(amplitude(solveCPar())==0.0f)) = Complex(ampave);
479 : }
480 : else
481 0 : solveCPar()=Complex(0.3);
482 :
483 0 : solveParOK()=True;
484 :
485 : /*
486 : cout << "Guess:" << endl
487 : << "amp = " << amplitude(solveCPar())
488 : << "phase = " << phase(solveCPar())
489 : << endl;
490 : */
491 0 : }
492 :
493 : // Fill the trivial DJ matrix elements
494 0 : void TJones::initTrivDJ() {
495 :
496 0 : if (prtlev()>4) cout << " T::initTrivDJ" << endl;
497 :
498 : // Must be trivial
499 0 : AlwaysAssert((trivialDJ()),AipsError);
500 :
501 : // This is the unit matrix
502 : // TBD: could we use a Jones::Unit type instead?
503 0 : diffJElem().resize(IPosition(4,1,1,1,1));
504 0 : diffJElem()=Complex(1.0);
505 :
506 0 : }
507 :
508 :
509 :
510 :
511 :
512 :
513 :
514 : //============================================================
515 :
516 :
517 0 : void TJones::createCorruptor(const VisIter& vi, const Record& simpar, const Int nSim) {
518 0 : LogIO os(LogOrigin("T", "createCorruptor()", WHERE));
519 :
520 0 : tcorruptor_p = new AtmosCorruptor(nSim);
521 0 : corruptor_p=tcorruptor_p;
522 :
523 : // call generic parent to set corr,spw,etc info
524 0 : SolvableVisCal::createCorruptor(vi,simpar,nSim);
525 :
526 0 : Int Seed(1234);
527 0 : if (simpar.isDefined("seed")) {
528 0 : Seed=simpar.asInt("seed");
529 : }
530 :
531 0 : Int rxType(0); // 0=2SB, 1=DSB
532 0 : if (simpar.isDefined("rxType")) {
533 0 : rxType=simpar.asInt("rxType");
534 : }
535 :
536 0 : Float Beta(1.1); // exponent for generalized 1/f noise
537 0 : if (simpar.isDefined("beta")) {
538 0 : Beta=simpar.asFloat("beta");
539 : }
540 :
541 0 : if (simpar.isDefined("mean_pwv"))
542 0 : tcorruptor_p->mean_pwv() = simpar.asFloat("mean_pwv");
543 :
544 0 : if (tcorruptor_p->mean_pwv()<=0)
545 0 : throw(AipsError("AtmCorruptor attempted initialization with undefined PWV"));
546 :
547 0 : if (simpar.isDefined("mode")) {
548 0 : if (prtlev()>2)
549 0 : cout << "initializing T:Corruptor with mode " << simpar.asString("mode") << endl;
550 0 : String simMode=simpar.asString("mode");
551 :
552 0 : if (simMode == "test")
553 0 : tcorruptor_p->initialize(rxType);
554 0 : else if (simMode == "individual" or simMode == "screen") {
555 :
556 0 : Float Scale(1.); // RELATIVE scale of fluctuations (to mean_pwv)
557 0 : if (simpar.isDefined("delta_pwv")) {
558 0 : if (simpar.asFloat("delta_pwv")>1.) {
559 0 : Scale=1.;
560 0 : os << LogIO::WARN << " decreasing PWV fluctuation magnitude to 100% of the mean PWV " << LogIO::POST;
561 : } else {
562 0 : Scale=simpar.asFloat("delta_pwv");
563 : }
564 : } else {
565 0 : os << LogIO::WARN << " setting PWV fluctuation magnitude to 15% of the mean PWV " << LogIO::POST;
566 0 : Scale=0.15;
567 : }
568 :
569 0 : os << " PWV fluctuations = " << Scale << " of mean PWV which is " << simpar.asFloat("mean_pwv") << "mm " << LogIO::POST;
570 :
571 :
572 : // slot_times for a fBM-based corruption need to be even even if solTimes are not
573 : // so will define startTime and stopTime and reset nsim() here.
574 :
575 0 : if (simpar.isDefined("startTime")) {
576 0 : corruptor_p->startTime() = simpar.asDouble("startTime");
577 : } else {
578 0 : throw(AipsError("start/stop time not defined"));
579 : }
580 0 : if (simpar.isDefined("stopTime")) {
581 0 : corruptor_p->stopTime() = simpar.asDouble("stopTime");
582 : } else {
583 0 : throw(AipsError("start/stop time not defined"));
584 : }
585 :
586 : // RI todo T::createCorr make min screen granularity a user parameter
587 0 : Float minFBM_interval = 0.1; // generate screens on >0.1s intervals
588 0 : Float defaultFBM_interval = 10.;
589 0 : Float fBM_interval=max(interval(), minFBM_interval);
590 0 : if (interval() <= 0.){ // use default value
591 0 : fBM_interval = defaultFBM_interval;
592 0 : } else if ( (fBM_interval - interval()) > 1E-5 ){
593 0 : os << LogIO::WARN << " Requested phase screen time granularity (" << interval()
594 0 : << " s) is too small! Will use the minimum permitted value: " << minFBM_interval << " s." << LogIO::POST;
595 : }
596 :
597 0 : corruptor_p->setEvenSlots(fBM_interval);
598 :
599 0 : if (simpar.asString("mode") == "individual")
600 0 : tcorruptor_p->initialize(Seed,Beta,Scale,rxType);
601 0 : else if (simpar.asString("mode") == "screen") {
602 0 : const MSAntennaColumns& antcols(vi.msColumns().antenna());
603 0 : if (simpar.isDefined("windspeed")) {
604 0 : tcorruptor_p->windspeed()=simpar.asFloat("windspeed");
605 0 : tcorruptor_p->initialize(Seed,Beta,Scale,rxType,antcols);
606 : } else
607 0 : throw(AipsError("Unknown wind speed for T:Corruptor"));
608 : }
609 :
610 0 : } else if (simMode == "tsys-atm" or simMode == "tsys-manual") {
611 : // NEW 20100818 change from Mf to Tf
612 : // M corruptor initialization didn't matter M or Mf here - it checks mode in
613 : // the Atmoscorruptor init.
614 0 : tcorruptor_p->initialize(vi,simpar,VisCal::T,rxType);
615 0 : extraTag()="NoiseScale"; // collapseForSim catches this
616 :
617 : } else
618 0 : throw(AipsError("Unknown mode for T:Corruptor"));
619 : } else {
620 0 : throw(AipsError("No Mode specified for T:Corruptor."));
621 : }
622 0 : }
623 :
624 :
625 :
626 :
627 :
628 : // **********************************************************
629 : // TfJones Implementations
630 : //
631 :
632 0 : TfJones::TfJones(VisSet& vs) :
633 : VisCal(vs), // virtual base
634 : VisMueller(vs), // virtual base
635 0 : TJones(vs) // immediate parent
636 : {
637 0 : if (prtlev()>2) cout << "Tf::Tf(vs)" << endl;
638 0 : }
639 :
640 0 : TfJones::TfJones(String msname,Int MSnAnt,Int MSnSpw) :
641 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
642 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
643 0 : TJones(msname,MSnAnt,MSnSpw) // immediate parent
644 : {
645 0 : if (prtlev()>2) cout << "Tf::Tf(msname,MSnAnt,MSnSpw)" << endl;
646 0 : }
647 :
648 0 : TfJones::TfJones(const MSMetaInfoForCal& msmc) :
649 : VisCal(msmc), // virtual base
650 : VisMueller(msmc), // virtual base
651 0 : TJones(msmc) // immediate parent
652 : {
653 0 : if (prtlev()>2) cout << "Tf::Tf(msmc)" << endl;
654 0 : }
655 :
656 :
657 0 : TfJones::TfJones(const Int& nAnt) :
658 : VisCal(nAnt),
659 : VisMueller(nAnt),
660 0 : TJones(nAnt)
661 : {
662 0 : if (prtlev()>2) cout << "Tf::Tf(nAnt)" << endl;
663 0 : }
664 :
665 0 : TfJones::~TfJones() {
666 0 : if (prtlev()>2) cout << "Tf::~Tf()" << endl;
667 0 : }
668 :
669 :
670 : // **********************************************************
671 : // GJones Implementations
672 : //
673 :
674 0 : GJones::GJones(VisSet& vs) :
675 : VisCal(vs), // virtual base
676 : VisMueller(vs), // virtual base
677 : SolvableVisJones(vs), // immediate parent
678 0 : gcorruptor_p(NULL)
679 : {
680 0 : if (prtlev()>2) cout << "G::G(vs)" << endl;
681 0 : }
682 :
683 0 : GJones::GJones(String msname,Int MSnAnt,Int MSnSpw) :
684 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
685 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
686 : SolvableVisJones(msname,MSnAnt,MSnSpw), // immediate parent
687 0 : gcorruptor_p(NULL)
688 : {
689 0 : if (prtlev()>2) cout << "G::G(msname,MSnAnt,MSnSpw)" << endl;
690 0 : }
691 :
692 0 : GJones::GJones(const MSMetaInfoForCal& msmc) :
693 : VisCal(msmc), // virtual base
694 : VisMueller(msmc), // virtual base
695 : SolvableVisJones(msmc), // immediate parent
696 0 : gcorruptor_p(NULL)
697 : {
698 0 : if (prtlev()>2) cout << "G::G(msmc)" << endl;
699 0 : }
700 :
701 0 : GJones::GJones(const Int& nAnt) :
702 : VisCal(nAnt),
703 : VisMueller(nAnt),
704 : SolvableVisJones(nAnt),
705 0 : gcorruptor_p(NULL)
706 : {
707 0 : if (prtlev()>2) cout << "G::G(nAnt)" << endl;
708 0 : }
709 :
710 0 : GJones::~GJones() {
711 0 : if (prtlev()>2) cout << "G::~G()" << endl;
712 0 : }
713 :
714 0 : void GJones::setSolve(const Record& solve) {
715 :
716 : // call parent to get general stuff
717 0 : SolvableVisJones::setSolve(solve);
718 :
719 : // parse solmode and rmsthresh
720 0 : solmode()=String("");
721 0 : if (solve.isDefined("solmode"))
722 0 : solmode()=solve.asString("solmode");
723 0 : solmode().upcase();
724 0 : rmsthresh().resize(0);
725 0 : if (solve.isDefined("rmsthresh")) {
726 0 : Vector<Double> rmsth;
727 0 : rmsth=solve.asArrayDouble("rmsthresh"); // parse from record as Vector<Double>
728 0 : Int nrms=rmsth.nelements();
729 0 : rmsthresh().resize(0);
730 0 : if (nrms>0) {
731 0 : rmsthresh().resize(nrms);
732 0 : convertArray(rmsthresh(), rmsth); // convert into Float array
733 : }
734 : }
735 0 : Int nrms=rmsthresh().nelements();
736 0 : if (nrms==0 && solmode().contains("R")) {
737 0 : rmsthresh().assign(Vector<Float>(std::vector<Float>({7.0,5.0,4.0,3.5,3.0,2.8,2.6,2.4,2.2,2.5})));
738 : }
739 :
740 0 : if (typeName()=="G Jones" || typeName()=="T Jones") {
741 0 : if (solmode().contains("L1"))
742 0 : logSink() << "Doing L1 solution." << LogIO::POST;
743 0 : if (solmode().contains("R"))
744 0 : logSink() << "Doing iterative outlier rejection using thresholds=" << rmsthresh() << LogIO::POST;
745 : }
746 :
747 :
748 0 : }
749 :
750 :
751 0 : void GJones::guessPar(VisBuffer& vb) {
752 :
753 0 : if (prtlev()>4) cout << " G::guessPar(vb)" << endl;
754 :
755 : // Make a guess at antenna-based G
756 : // Correlation membership-dependencexm
757 : // assumes corrs in canonical order
758 : // nCorr = 1: use icorr=0
759 : // nCorr = 2: use icorr=[0,1]
760 : // nCorr = 4: use icorr=[0,3]
761 :
762 0 : Int nCorr(2);
763 0 : Int nDataCorr(vb.visCube().shape()(0));
764 0 : Vector<Int> corridx(nCorr,0);
765 0 : if (nDataCorr==2) {
766 0 : corridx(0)=0;
767 0 : corridx(1)=1;
768 : }
769 0 : else if (nDataCorr==4) {
770 0 : corridx(0)=0;
771 0 : corridx(1)=3;
772 : }
773 :
774 0 : Int guesschan(vb.visCube().shape()(1)-1);
775 :
776 : // cout << "guesschan = " << guesschan << endl;
777 : // cout << "nCorr = " << nCorr << endl;
778 : // cout << "corridx = " << corridx << endl;
779 :
780 :
781 : // Find out which ants are available
782 0 : Vector<Int> antok(nAnt(),0);
783 0 : Vector<Bool> rowok(vb.nRow(),false);
784 0 : for (Int irow=0;irow<vb.nRow();++irow) {
785 : // Is this row ok
786 0 : rowok(irow)= (!vb.flagRow()(irow) &&
787 0 : vb.antenna1()(irow)!=vb.antenna2()(irow) &&
788 0 : nfalse(vb.flag().column(irow))> 0 );
789 0 : if (rowok(irow)) {
790 0 : antok(vb.antenna1()(irow))++;
791 0 : antok(vb.antenna2()(irow))++;
792 : }
793 : }
794 :
795 : // Assume refant is the target ant, for starters
796 0 : Int guessant(refant());
797 : // Int guessant(-1);
798 :
799 : // If no refant specified, or no data for refant
800 : // base first guess on first good ant
801 0 : if (guessant<0 || antok(guessant)<1) {
802 0 : guessant=0;
803 0 : while (antok(guessant)<1) guessant++;
804 : }
805 :
806 : // cout << "antok = " << antok << endl;
807 :
808 : // cout << "guessant = " << guessant << " (" << currSpw() << ")" << endl;
809 :
810 0 : AlwaysAssert(guessant>-1,AipsError);
811 :
812 0 : Cube<Complex>& V(vb.visCube());
813 0 : Float amp(0.0),ampave(0.0);
814 0 : Int namp(0);
815 0 : solveCPar()=Complex(0.0);
816 0 : for (Int irow=1;irow<vb.nRow();++irow) {
817 :
818 0 : if (rowok(irow) && !vb.flag()(guesschan,irow)) {
819 0 : Int a1=vb.antenna1()(irow);
820 0 : Int a2=vb.antenna2()(irow);
821 :
822 : // If this row contains the guessant
823 0 : if (a1 == guessant || a2==guessant) {
824 :
825 : // cout << a1 << " " << a2 << " "
826 : // << vb.visCube().xyPlane(irow).column(guesschan) << " "
827 : // << amplitude(vb.visCube().xyPlane(irow).column(guesschan)) << " "
828 : // << endl;
829 :
830 0 : for (Int icorr=0;icorr<nCorr;icorr++) {
831 0 : Complex& Vi(V(corridx(icorr),guesschan,irow));
832 0 : amp=abs(Vi);
833 0 : if (amp>0.0f) {
834 0 : if (a1 == guessant)
835 0 : solveCPar()(icorr,0,a2)=conj(Vi);
836 : else
837 0 : solveCPar()(icorr,0,a1)=(Vi);
838 :
839 0 : ampave+=amp;
840 0 : namp++;
841 : }
842 : }
843 : } // guessant
844 : } // rowok
845 : } // irow
846 :
847 : // Scale them by the mean amplitude
848 :
849 0 : if (namp>0) {
850 0 : ampave/=Float(namp);
851 0 : ampave=sqrt(ampave);
852 : // solveCPar()*=Complex(ampave);
853 0 : solveCPar()/=Complex(ampave);
854 0 : solveCPar()(0,0,guessant)=solveCPar()(1,0,guessant)=Complex(ampave);
855 0 : solveCPar()(LogicalArray(amplitude(solveCPar())==0.0f)) = Complex(ampave);
856 : }
857 : else
858 0 : solveCPar()=Complex(0.3);
859 :
860 0 : solveParOK()=true;
861 :
862 : //For scalar data, Set "other" pol soln to zero
863 0 : if (nDataCorr == 1)
864 0 : solveCPar()(IPosition(3,1,0,0),IPosition(3,1,0,nAnt()-1))=Complex(0.0);
865 :
866 : /*
867 : cout << "Guess:" << endl;
868 : cout << "amplitude(solveCPar()) = " << amplitude(solveCPar()) << endl;
869 : cout << "phases = " << phase(solveCPar())*180.0/C::pi << endl;
870 : cout << "solveParOK() = " << solveParOK() << endl;
871 : */
872 0 : }
873 :
874 0 : void GJones::guessPar(SDBList& sdbs, const Bool& corrDepFlags) {
875 :
876 0 : if (prtlev()>4) cout << " G::guessPar(sdbs)" << endl;
877 :
878 : // Use just the first SDB in the SDBList, for now
879 0 : SolveDataBuffer& sdb(sdbs(0));
880 :
881 : // Make a guess at antenna-based G
882 : // Correlation membership-dependencexm
883 : // assumes corrs in canonical order
884 : // nCorr = 1: use icorr=0
885 : // nCorr = 2: use icorr=[0,1]
886 : // nCorr = 4: use icorr=[0,3]
887 :
888 0 : Int nCorr(2);
889 0 : Int nDataCorr(sdb.nCorrelations());
890 0 : Vector<Int> corridx(nCorr,0);
891 0 : if (nDataCorr==2) {
892 0 : corridx(0)=0;
893 0 : corridx(1)=1;
894 : }
895 0 : else if (nDataCorr==4) {
896 0 : corridx(0)=0;
897 0 : corridx(1)=3;
898 : }
899 :
900 0 : Int guesschan(sdb.nChannels()-1);
901 :
902 : /*
903 : cout << "guesschan = " << guesschan << endl;
904 : cout << "nCorr = " << nCorr << endl;
905 : cout << "corridx = " << corridx << endl;
906 : */
907 :
908 : // Find out which ants are available
909 0 : Int nRow=sdb.nRows();
910 0 : Vector<Int> antok(nAnt(),0);
911 0 : Vector<Bool> rowok(nRow,False);
912 0 : for (Int irow=0;irow<nRow;++irow) {
913 :
914 0 : const Int& a1(sdb.antenna1()(irow));
915 0 : const Int& a2(sdb.antenna2()(irow));
916 :
917 : // Is this row ok?
918 0 : rowok(irow)= (!sdb.flagRow()(irow) &&
919 0 : a1!=a2);
920 :
921 0 : if (!corrDepFlags) {
922 : // All relevant correlations must be good
923 0 : for (Int icorr=0;icorr<nCorr;++icorr)
924 0 : rowok(irow)=(rowok(irow) && !sdb.flagCube()(corridx[icorr],guesschan,irow));
925 : }
926 :
927 0 : if (rowok(irow)) {
928 0 : antok(a1)++;
929 0 : antok(a2)++;
930 : }
931 : }
932 :
933 : // Assume refant is the target ant, for starters
934 0 : Int guessant(refant());
935 : // Int guessant(-1);
936 :
937 : // If no refant specified, or no data for refant
938 : // base first guess on first good ant
939 0 : if (guessant<0 || antok(guessant)<1) {
940 0 : guessant=0;
941 0 : while (antok(guessant)<1) guessant++;
942 : }
943 :
944 : // cout << "antok = " << antok << endl;
945 :
946 : // cout << "guessant = " << guessant << " (" << currSpw() << ")" << endl;
947 :
948 0 : AlwaysAssert(guessant>-1,AipsError);
949 :
950 0 : const Cube<Complex>& V(sdb.visCubeCorrected());
951 0 : Float amp(0.0),ampave(0.0);
952 0 : Int namp(0);
953 0 : solveCPar()=Complex(0.0);
954 0 : for (Int irow=0;irow<nRow;++irow) {
955 :
956 0 : if (rowok(irow)) {
957 0 : const Int& a1=sdb.antenna1()(irow);
958 0 : const Int& a2=sdb.antenna2()(irow);
959 :
960 : // If this row contains the guessant
961 0 : if (a1 == guessant || a2==guessant) {
962 :
963 : // cout << a1 << " " << a2 << " "
964 : // << sdb.visCube().xyPlane(irow).column(guesschan) << " "
965 : // << amplitude(sdb.visCube().xyPlane(irow).column(guesschan)) << " "
966 : // << endl;
967 :
968 0 : for (Int icorr=0;icorr<nCorr;icorr++) {
969 0 : const Complex& Vi(V(corridx(icorr),guesschan,irow));
970 0 : amp=abs(Vi);
971 0 : if (amp>0.0f) {
972 0 : if (a1 == guessant)
973 0 : solveCPar()(icorr,0,a2)=conj(Vi);
974 : else
975 0 : solveCPar()(icorr,0,a1)=(Vi);
976 :
977 0 : ampave+=amp;
978 0 : namp++;
979 : }
980 : }
981 : } // guessant
982 : } // rowok
983 : } // irow
984 :
985 : // Scale them by the mean amplitude
986 :
987 0 : if (namp>0) {
988 0 : ampave/=Float(namp);
989 0 : ampave=sqrt(ampave);
990 : // solveCPar()*=Complex(ampave);
991 0 : solveCPar()/=Complex(ampave);
992 0 : solveCPar()(0,0,guessant)=solveCPar()(1,0,guessant)=Complex(ampave);
993 0 : solveCPar()(LogicalArray(amplitude(solveCPar())==0.0f)) = Complex(ampave);
994 : }
995 : else
996 0 : solveCPar()=Complex(0.3);
997 :
998 0 : solveParOK()=True;
999 :
1000 : //For scalar data, Set "other" pol soln to zero
1001 0 : if (nDataCorr == 1)
1002 0 : solveCPar()(IPosition(3,1,0,0),IPosition(3,1,0,nAnt()-1))=Complex(0.0);
1003 :
1004 : /*
1005 : cout << "Guess:" << endl;
1006 : cout << "amplitude(solveCPar()) = " << amplitude(solveCPar()) << endl;
1007 : cout << "phases = " << phase(solveCPar())*180.0/C::pi << endl;
1008 : cout << "solveParOK() = " << solveParOK() << endl;
1009 : */
1010 0 : }
1011 :
1012 : // Fill the trivial DJ matrix elements
1013 0 : void GJones::initTrivDJ() {
1014 :
1015 0 : if (prtlev()>4) cout << " G::initTrivDJ" << endl;
1016 :
1017 : // Must be trivial
1018 0 : AlwaysAssert((trivialDJ()),AipsError);
1019 :
1020 : // 1 0 0 0
1021 : // 0 0 & 0 1
1022 : //
1023 0 : if (diffJElem().nelements()==0) {
1024 0 : diffJElem().resize(IPosition(4,2,2,1,1));
1025 0 : diffJElem()=0.0;
1026 0 : diffJElem()(IPosition(4,0,0,0,0))=Complex(1.0);
1027 0 : diffJElem()(IPosition(4,1,1,0,0))=Complex(1.0);
1028 : }
1029 :
1030 0 : }
1031 :
1032 :
1033 0 : void GJones::createCorruptor(const VisIter& vi, const Record& simpar, const Int nSim) {
1034 : {
1035 :
1036 0 : LogIO os(LogOrigin("G", "createCorruptor()", WHERE));
1037 0 : if (prtlev()>2) cout << " G::createCorruptor()" << endl;
1038 :
1039 0 : gcorruptor_p = new GJonesCorruptor(nSim);
1040 0 : corruptor_p = gcorruptor_p;
1041 :
1042 : // call generic parent to set corr,spw,etc info
1043 0 : SolvableVisCal::createCorruptor(vi,simpar,nSim);
1044 :
1045 0 : Int Seed(1234);
1046 0 : if (simpar.isDefined("seed")) {
1047 0 : Seed=simpar.asInt("seed");
1048 : }
1049 :
1050 0 : if (simpar.isDefined("tsys")) {
1051 0 : gcorruptor_p->tsys() = simpar.asFloat("tsys");
1052 : }
1053 :
1054 0 : if (simpar.isDefined("mode")) {
1055 0 : if (prtlev()>2)
1056 0 : cout << "initializing GCorruptor with mode " << simpar.asString("mode") << endl;
1057 :
1058 : // slot_times for a fBM-based corruption need to be even even if solTimes are not
1059 : // so will define startTime and stopTime and reset nsim() here.
1060 :
1061 0 : if (simpar.isDefined("startTime")) {
1062 0 : corruptor_p->startTime() = simpar.asDouble("startTime");
1063 : } else {
1064 0 : throw(AipsError("start/stop time not defined"));
1065 : }
1066 0 : if (simpar.isDefined("stopTime")) {
1067 0 : corruptor_p->stopTime() = simpar.asDouble("stopTime");
1068 : } else {
1069 0 : throw(AipsError("start/stop time not defined"));
1070 : }
1071 :
1072 0 : if (simpar.asString("mode")=="fbm") {
1073 :
1074 0 : Float Beta(1.1); // exponent for generalized 1/f noise
1075 0 : if (simpar.isDefined("beta")) {
1076 0 : Beta=simpar.asFloat("beta");
1077 : }
1078 :
1079 0 : Float Scale(.15); // scale of fluctuations
1080 0 : if (simpar.isDefined("amplitude")) {
1081 0 : Scale=simpar.asFloat("amplitude");
1082 0 : if (Scale>=.9) {
1083 0 : os << LogIO::WARN << " decreasing gain fluctuations from " << Scale << " to 0.9 " << LogIO::POST;
1084 0 : Scale=.9;
1085 : }
1086 : }
1087 :
1088 0 : Float fBM_interval=max(interval(),5.); // generate screens on 5s intervals or longer
1089 0 : corruptor_p->setEvenSlots(fBM_interval);
1090 0 : gcorruptor_p->initialize(Seed,Beta,Scale);
1091 :
1092 0 : } else if (simpar.asString("mode")=="random") {
1093 :
1094 0 : Complex Scale(0.1,0.1); // scale of fluctuations
1095 0 : if (simpar.isDefined("camp")) {
1096 0 : Scale=simpar.asComplex("camp");
1097 : }
1098 0 : gcorruptor_p->initialize(Seed,Scale);
1099 :
1100 0 : } else throw AipsError("incompatible mode "+simpar.asString("mode"));
1101 :
1102 :
1103 : } else
1104 0 : throw(AipsError("Unknown mode for GJonesCorruptor"));
1105 : }
1106 0 : }
1107 :
1108 :
1109 :
1110 :
1111 :
1112 :
1113 :
1114 : // **********************************************************
1115 : // BJones Implementations
1116 : //
1117 :
1118 0 : BJones::BJones(VisSet& vs) :
1119 : VisCal(vs), // virtual base
1120 : VisMueller(vs), // virtual base
1121 : GJones(vs), // immediate parent
1122 0 : maxchangap_p(0)
1123 : {
1124 0 : if (prtlev()>2) cout << "B::B(vs)" << endl;
1125 0 : }
1126 :
1127 0 : BJones::BJones(String msname,Int MSnAnt,Int MSnSpw) :
1128 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
1129 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
1130 : GJones(msname,MSnAnt,MSnSpw), // immediate parent
1131 0 : maxchangap_p(0)
1132 : {
1133 0 : if (prtlev()>2) cout << "B::B(msname,MSnAnt,MSnSpw)" << endl;
1134 0 : }
1135 :
1136 0 : BJones::BJones(const MSMetaInfoForCal& msmc) :
1137 : VisCal(msmc), // virtual base
1138 : VisMueller(msmc), // virtual base
1139 : GJones(msmc), // immediate parent
1140 0 : maxchangap_p(0)
1141 : {
1142 0 : if (prtlev()>2) cout << "B::B(msmc)" << endl;
1143 0 : }
1144 :
1145 0 : BJones::BJones(const Int& nAnt) :
1146 : VisCal(nAnt),
1147 : VisMueller(nAnt),
1148 : GJones(nAnt),
1149 0 : maxchangap_p(0)
1150 : {
1151 0 : if (prtlev()>2) cout << "B::B(nAnt)" << endl;
1152 0 : }
1153 :
1154 0 : BJones::~BJones() {
1155 0 : if (prtlev()>2) cout << "B::~B()" << endl;
1156 0 : }
1157 :
1158 0 : void BJones::setSolve(const Record& solve) {
1159 :
1160 : // call parent to get general stuff
1161 0 : GJones::setSolve(solve);
1162 :
1163 0 : if (solmode()!="") {
1164 0 : solmode()="";
1165 0 : cout << "solmode options not yet supported for B solutions; ignoring." << endl;
1166 : }
1167 :
1168 : // get max chan gap from user
1169 0 : maxchangap_p=0;
1170 0 : if (solve.isDefined("maxgap"))
1171 0 : maxchangap_p=solve.asInt("maxgap");
1172 :
1173 0 : }
1174 :
1175 0 : void BJones::normalize() {
1176 :
1177 : // Only if we have a CalTable, and it is not empty
1178 0 : if (ct_ && ct_->nrow()>0) {
1179 :
1180 : // TBD: trap attempts to normalize a caltable containing FPARAM (non-Complex)?
1181 :
1182 0 : logSink() << "Normalizing solutions per spw, time, ant, pol with " << solNorm().normtypeString()
1183 : << " in amplitude (and center channel in phase)."
1184 0 : << LogIO::POST;
1185 :
1186 : // Bandpass is normalized per spw, time, antenna, and pol
1187 0 : Block<String> col(3);
1188 0 : col[0]="SPECTRAL_WINDOW_ID";
1189 0 : col[1]="TIME";
1190 0 : col[2]="ANTENNA1";
1191 0 : CTIter ctiter(*ct_,col);
1192 :
1193 : // Cube iteration axes are pol and ant (ant is degenerate)
1194 0 : IPosition itax(2,0,2);
1195 :
1196 0 : while (!ctiter.pastEnd()) {
1197 0 : Cube<Bool> fl(ctiter.flag());
1198 :
1199 0 : if (nfalse(fl)>0) {
1200 0 : Cube<Complex> p(ctiter.cparam());
1201 0 : ArrayIterator<Complex> soliter(p,itax,false);
1202 0 : ArrayIterator<Bool> fliter(fl,itax,false);
1203 0 : Int ipol(0);
1204 0 : while (!soliter.pastEnd()) {
1205 0 : Complex normfactor=normSolnArray(soliter.array(),!fliter.array(),true); // Do phase
1206 0 : logSink() << " Normalization factor (" << solNorm().normtypeString() << ") for"
1207 : << " spw=" << ctiter.thisSpw()
1208 0 : << " time=" << MVTime(ctiter.thisTime()/C::day).string(MVTime::YMD,7)
1209 : << " ant=" << ctiter.thisAntenna1()
1210 : << " pol " << ipol%2
1211 0 : << " = " << abs(normfactor) << ", " << arg(normfactor)*180.0/C::pi << "deg"
1212 0 : << LogIO::POST;
1213 0 : soliter.next();
1214 0 : fliter.next();
1215 0 : ++ipol;
1216 : }
1217 :
1218 : // record result...
1219 0 : ctiter.setcparam(p);
1220 : }
1221 0 : ctiter.next();
1222 : }
1223 : }
1224 :
1225 0 : }
1226 :
1227 0 : void BJones::globalPostSolveTinker() {
1228 :
1229 : // Call parent to do more general things
1230 0 : SolvableVisJones::globalPostSolveTinker();
1231 :
1232 : // Fill gaps in channels, if necessary
1233 0 : if (maxchangap_p>0)
1234 0 : fillChanGaps();
1235 :
1236 0 : }
1237 :
1238 0 : void BJones::fillChanGaps() {
1239 :
1240 : // TBD: Can this be consolidated with normalization (should be done before norm!)
1241 :
1242 : logSink() << "Filling in flagged solution channels by interpolation."
1243 0 : << LogIO::POST;
1244 :
1245 : // Iteration axes (norm per spw, pol, ant, timestamp)
1246 0 : IPosition itax(3,0,2,3);
1247 :
1248 :
1249 : // Only if we have a CalTable, and it is not empty
1250 0 : if (ct_ && ct_->nrow()>0) {
1251 :
1252 : // TBD: trap attempts to normalize a caltable containing FPARAM (non-Complex)?
1253 :
1254 : logSink() << "Normalizing solutions per spw, pol, ant, time."
1255 0 : << LogIO::POST;
1256 :
1257 : // In this generic version, one normalization factor per spw
1258 0 : Block<String> col(3);
1259 0 : col[0]="SPECTRAL_WINDOW_ID";
1260 0 : col[1]="TIME";
1261 0 : col[2]="ANTENNA1";
1262 0 : CTIter ctiter(*ct_,col);
1263 :
1264 : // Cube iteration axes are pol and ant
1265 0 : IPosition itax(2,0,2);
1266 :
1267 0 : while (!ctiter.pastEnd()) {
1268 0 : Cube<Bool> fl(ctiter.flag());
1269 :
1270 0 : if (nfalse(fl)>0) {
1271 0 : Cube<Complex> p(ctiter.cparam());
1272 0 : ArrayIterator<Complex> soliter(p,itax,false);
1273 0 : ArrayIterator<Bool> fliter(fl,itax,false);
1274 0 : Array<Bool> sok;
1275 0 : while (!soliter.pastEnd()) {
1276 0 : Array<Bool> thfl(fliter.array());
1277 0 : sok.assign(!thfl);
1278 0 : fillChanGapArray(soliter.array(),sok);
1279 0 : thfl=!sok; // sok is revised (shapes necessarily match)
1280 0 : soliter.next();
1281 0 : fliter.next();
1282 : }
1283 :
1284 : // record result...
1285 0 : ctiter.setcparam(p);
1286 0 : ctiter.setflag(fl);
1287 : }
1288 0 : ctiter.next();
1289 : }
1290 : }
1291 0 : }
1292 :
1293 0 : void BJones::fillChanGapArray(Array<Complex>& sol,
1294 : Array<Bool>& solOK) {
1295 :
1296 : // TBD: Do this with InterpolateArray1D a la CTPatchedInterp::resampleInFreq
1297 :
1298 : // Make the arrays 1D
1299 0 : Vector<Complex> solv(sol.reform(IPosition(1,sol.nelements())));
1300 0 : Vector<Bool> solOKv(solOK.reform(IPosition(1,solOK.nelements())));
1301 :
1302 0 : Int nChan(solv.nelements());
1303 0 : Bool done(false);
1304 0 : Int ich(0), ch1(-1), ch2(-1);
1305 0 : Int dch(1);
1306 : Float a0, da, a, p0, dp, p, fch;
1307 :
1308 : // Advance to first unflagged channel
1309 0 : while(!solOKv(ich) && ich<nChan) ++ich;
1310 :
1311 : // Found no unflagged channels, so signal escape
1312 0 : if (ich==nChan) done=true;
1313 :
1314 : // done turns true if we reach nChan, and nothing more to do
1315 0 : while (!done) {
1316 :
1317 : // Advance to next flagged channel
1318 0 : while(solOKv(ich) && ich<nChan) ++ich;
1319 :
1320 0 : if (ich<nChan) {
1321 :
1322 : // Left boundary of gap
1323 0 : ch1=ich-1; // (NB: above logic prevents ch1 < 0)
1324 :
1325 : // Find right boundary:
1326 0 : ch2=ich+1;
1327 0 : while (!solOKv(ch2) && ch2<nChan) ++ch2;
1328 :
1329 0 : if (ch2<nChan) {
1330 :
1331 : // The span of the interpolation (in channels)
1332 0 : dch=ch2-ch1;
1333 :
1334 : //cout << ch1 << " " << ch2 << " " << dch << endl;
1335 :
1336 : // Interpolate only if gap is narrower than maxchangap_p
1337 0 : if (dch<=maxchangap_p+1) {
1338 :
1339 : // calculate interp params
1340 0 : a0=abs(solv(ch1));
1341 0 : da=abs(solv(ch2))-a0;
1342 0 : p0=arg(solv(ch1));
1343 0 : dp=arg(solv(ch2))-p0;
1344 0 : if (dp>C::pi) dp-=(2*C::pi);
1345 0 : if (dp<-C::pi) dp+=(2*C::pi);
1346 :
1347 : //cout << a0 << " " << da << " " << p0 << " " << dp << endl;
1348 :
1349 :
1350 : // interpolate the intervening channels
1351 0 : while (ich<ch2) {
1352 0 : fch=Float(ich-ch1)/Float(dch);
1353 0 : a=a0 + da*fch;
1354 0 : p=p0 + dp*fch;
1355 :
1356 : // cout << " " << ich << " " << a << " " << p << endl;
1357 :
1358 0 : solv(ich)=a*Complex(cos(p),sin(p));
1359 0 : solOKv(ich)=true;
1360 0 : ++ich;
1361 : }
1362 :
1363 : // Begin looking for new gap on next round
1364 0 : ++ich;
1365 :
1366 : }
1367 : else
1368 : // we skipped a gap, look beyond it on next round
1369 0 : ich=ch2+1;
1370 : }
1371 : else
1372 : // Reached nChan looking for ch2
1373 0 : done=true;
1374 : }
1375 : else
1376 : // Reach nChan looking for gaps
1377 0 : done=true;
1378 : } // done
1379 :
1380 0 : }
1381 :
1382 0 : void BJones::calcWtScale() {
1383 :
1384 : // cout << "BJones::calcWtScale()" << endl;
1385 :
1386 :
1387 : // Access pre-chan-interp'd bandpass amplitude/flags
1388 0 : Cube<Float> amps;
1389 0 : Cube<Bool> ampfl;
1390 :
1391 :
1392 0 : if (cpp_) {
1393 0 : Cube<Float> tr;
1394 0 : cpp_->getTresult(tr,ampfl,currObs(),currScan(),currField(),-1,currSpw());
1395 0 : amps.reference(tr(Slice(0,2,2),Slice(),Slice()));
1396 : }
1397 0 : else if (ci_) {
1398 0 : amps=Cube<Float>(ci_->tresultF(currObs(),currScan(),currField(),currSpw()))(Slice(0,2,2),Slice(),Slice());
1399 0 : ampfl=ci_->tresultFlag(currObs(),currScan(),currField(),currSpw());
1400 : }
1401 : else
1402 : // BPOLY?
1403 0 : throw(AipsError("Can't BJones::calcWtScale because there is no solution interpolation...."));
1404 :
1405 :
1406 : // Initialize
1407 0 : currWtScale()=1.0;
1408 :
1409 :
1410 0 : IPosition ash(amps.shape());
1411 0 : ash(1)=1; // only one channel in the weights
1412 0 : Cube<Float> cWS(currWtScale());
1413 :
1414 0 : IPosition it3(2,0,2);
1415 0 : ArrayIterator<Float> A(amps,it3,false);
1416 0 : ArrayIterator<Bool> Aok(ampfl,it3,false);
1417 0 : ArrayIterator<Float> cWSi(cWS,it3,false);
1418 :
1419 0 : while (!A.pastEnd()) {
1420 :
1421 : // the weight scale factor is just the square of the
1422 : // freq-dep normalization
1423 0 : cWSi.array()*=(float)pow(calcPowerNorm(A.array(),!Aok.array()),2);
1424 :
1425 0 : A.next();
1426 0 : Aok.next();
1427 0 : cWSi.next();
1428 :
1429 : }
1430 :
1431 0 : }
1432 :
1433 :
1434 :
1435 :
1436 :
1437 : // **********************************************************
1438 : // JJones Implementations
1439 : //
1440 :
1441 0 : JJones::JJones(VisSet& vs) :
1442 : VisCal(vs), // virtual base
1443 : VisMueller(vs), // virtual base
1444 0 : SolvableVisJones(vs) // immediate parent
1445 : {
1446 0 : if (prtlev()>2) cout << "J::J(vs)" << endl;
1447 0 : }
1448 :
1449 0 : JJones::JJones(String msname,Int MSnAnt,Int MSnSpw) :
1450 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
1451 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
1452 0 : SolvableVisJones(msname,MSnAnt,MSnSpw) // immediate parent
1453 : {
1454 0 : if (prtlev()>2) cout << "J::J(msname,MSnAnt,MSnSpw)" << endl;
1455 0 : }
1456 :
1457 0 : JJones::JJones(const MSMetaInfoForCal& msmc) :
1458 : VisCal(msmc), // virtual base
1459 : VisMueller(msmc), // virtual base
1460 0 : SolvableVisJones(msmc) // immediate parent
1461 : {
1462 0 : if (prtlev()>2) cout << "J::J(msmc)" << endl;
1463 0 : }
1464 :
1465 0 : JJones::JJones(const Int& nAnt) :
1466 : VisCal(nAnt),
1467 : VisMueller(nAnt),
1468 0 : SolvableVisJones(nAnt)
1469 : {
1470 0 : if (prtlev()>2) cout << "J::J(nAnt)" << endl;
1471 0 : }
1472 :
1473 0 : JJones::~JJones() {
1474 0 : if (prtlev()>2) cout << "J::~J()" << endl;
1475 0 : }
1476 :
1477 0 : void JJones::setSolve(const Record& solvepar) {
1478 :
1479 : // Call parent
1480 0 : SolvableVisJones::setSolve(solvepar);
1481 :
1482 : // For J insist preavg is meaningful (5 minutes or user-supplied)
1483 0 : if (preavg()<0.0)
1484 0 : preavg()=300.0;
1485 :
1486 0 : }
1487 :
1488 :
1489 0 : void JJones::guessPar(VisBuffer& vb) {
1490 :
1491 0 : if (prtlev()>4) cout << " ::guessPar(vb)" << endl;
1492 :
1493 : // Make a guess at antenna-based J
1494 : // Correlation membership-dependence
1495 : // assumes corrs in canonical order
1496 : // nCorr = 1: use icorr=0
1497 : // nCorr = 2: use icorr=[0,1]
1498 : // nCorr = 4: use icorr=[0,3]
1499 :
1500 : // This method sets the off-diag = 0.0,
1501 : // and the on-diag as if this were G
1502 :
1503 0 : Int nCorr(2);
1504 0 : Int nDataCorr(vb.visCube().shape()(0));
1505 0 : Vector<Int> corridx(nCorr,0);
1506 0 : if (nDataCorr==2) {
1507 0 : corridx(0)=0;
1508 0 : corridx(1)=1;
1509 : }
1510 0 : else if (nDataCorr==4) {
1511 0 : corridx(0)=0;
1512 0 : corridx(1)=3;
1513 : }
1514 :
1515 0 : Cube<Complex>& V(vb.visCube());
1516 0 : Float amp(0.0),ampave(0.0);
1517 0 : Int namp(0);
1518 0 : solveCPar()=Complex(0.0);
1519 0 : for (Int irow=1;irow<nAnt();++irow) {
1520 :
1521 0 : for (Int icorr=0;icorr<nCorr;icorr++) {
1522 0 : Complex& Vi(V(corridx(icorr),0,irow));
1523 0 : amp=abs(Vi);
1524 0 : if (amp>0.0f) {
1525 0 : solveCPar()(3*icorr,0,irow)=(conj(Vi)/amp);
1526 0 : ampave+=amp;
1527 0 : namp++;
1528 : // cout << " " << abs(Vi) << " " << arg(Vi)*180.0/C::pi << endl;
1529 : }
1530 : }
1531 :
1532 : }
1533 :
1534 : // Scale them by the mean amplitude
1535 0 : ampave/=Float(namp);
1536 0 : ampave=sqrt(ampave);
1537 0 : solveCPar()*=Complex(ampave);
1538 0 : solveParOK()=true;
1539 :
1540 : // cout << "post-guess:" << endl;
1541 : // cout << "solveCPar() = " << solveCPar() << endl;
1542 : // cout << "phases = " << phase(solveCPar())*180.0/C::pi << endl;
1543 : // cout << "solveParOK() = " << solveParOK() << endl;
1544 :
1545 0 : }
1546 :
1547 : // Fill the trivial DJ matrix elements
1548 0 : void JJones::initTrivDJ() {
1549 :
1550 0 : if (prtlev()>4) cout << " J::initTrivDJ" << endl;
1551 :
1552 : // Must be trivial
1553 0 : AlwaysAssert((trivialDJ()),AipsError);
1554 :
1555 : // 1 0 0 1 0 0 0 0
1556 : // 0 0 0 0 1 0 0 1
1557 :
1558 0 : diffJElem().resize(IPosition(4,4,4,1,1));
1559 0 : diffJElem()=0.0;
1560 0 : diffJElem()(IPosition(4,0,0,0,0))=Complex(1.0);
1561 0 : diffJElem()(IPosition(4,1,1,0,0))=Complex(1.0);
1562 0 : diffJElem()(IPosition(4,2,2,0,0))=Complex(1.0);
1563 0 : diffJElem()(IPosition(4,3,3,0,0))=Complex(1.0);
1564 :
1565 0 : }
1566 :
1567 :
1568 : // **********************************************************
1569 : // JfJones Implementations
1570 : //
1571 :
1572 0 : JfJones::JfJones(VisSet& vs) :
1573 : VisCal(vs), // virtual base
1574 : VisMueller(vs), // virtual base
1575 0 : JJones(vs) // immediate parent
1576 : {
1577 0 : if (prtlev()>2) cout << "Jf::Jf(vs)" << endl;
1578 0 : }
1579 :
1580 0 : JfJones::JfJones(String msname,Int MSnAnt,Int MSnSpw) :
1581 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
1582 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
1583 0 : JJones(msname,MSnAnt,MSnSpw) // immediate parent
1584 : {
1585 0 : if (prtlev()>2) cout << "Jf::Jf(msname,MSnAnt,MSnSpw)" << endl;
1586 0 : }
1587 :
1588 0 : JfJones::JfJones(const MSMetaInfoForCal& msmc) :
1589 : VisCal(msmc), // virtual base
1590 : VisMueller(msmc), // virtual base
1591 0 : JJones(msmc) // immediate parent
1592 : {
1593 0 : if (prtlev()>2) cout << "Jf::Jf(msmc)" << endl;
1594 0 : }
1595 :
1596 0 : JfJones::JfJones(const Int& nAnt) :
1597 : VisCal(nAnt),
1598 : VisMueller(nAnt),
1599 0 : JJones(nAnt)
1600 : {
1601 0 : if (prtlev()>2) cout << "Jf::Jf(nAnt)" << endl;
1602 0 : }
1603 :
1604 0 : JfJones::~JfJones() {
1605 0 : if (prtlev()>2) cout << "Jf::~Jf()" << endl;
1606 0 : }
1607 :
1608 :
1609 : /*
1610 :
1611 : // **********************************************************
1612 : // MMueller: baseline-based (closure) solution
1613 : //
1614 :
1615 : MMueller::MMueller(VisSet& vs) :
1616 : VisCal(vs), // virtual base
1617 : VisMueller(vs), // virtual base
1618 : SolvableVisMueller(vs) // immediate parent
1619 : {
1620 : if (prtlev()>2) cout << "M::M(vs)" << endl;
1621 : }
1622 :
1623 : MMueller::MMueller(String msname,Int MSnAnt,Int MSnSpw) :
1624 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
1625 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
1626 : SolvableVisMueller(msname,MSnAnt,MSnSpw) // immediate parent
1627 : {
1628 : if (prtlev()>2) cout << "M::M(msname,MSnAnt,MSnSpw)" << endl;
1629 : }
1630 :
1631 : MMueller::MMueller(const MSMetaInfoForCal& msmc) :
1632 : VisCal(msmc), // virtual base
1633 : VisMueller(msmc), // virtual base
1634 : SolvableVisMueller(msmc) // immediate parent
1635 : {
1636 : if (prtlev()>2) cout << "M::M(msmc)" << endl;
1637 : }
1638 :
1639 : MMueller::MMueller(const Int& nAnt) :
1640 : VisCal(nAnt),
1641 : VisMueller(nAnt),
1642 : SolvableVisMueller(nAnt)
1643 : {
1644 : if (prtlev()>2) cout << "M::M(nAnt)" << endl;
1645 : }
1646 :
1647 : MMueller::~MMueller() {
1648 : if (prtlev()>2) cout << "M::~M()" << endl;
1649 : }
1650 :
1651 : void MMueller::setApply(const Record& apply) {
1652 :
1653 : SolvableVisCal::setApply(apply);
1654 :
1655 : // Force calwt to false for now
1656 : calWt()=false;
1657 :
1658 : }
1659 :
1660 : void MMueller::newselfSolve(VisSet& vs, VisEquation& ve) {
1661 :
1662 : if (prtlev()>4) cout << " M::selfSolve(ve)" << endl;
1663 :
1664 : // Inform logger/history
1665 : logSink() << "Solving for " << typeName()
1666 : << LogIO::POST;
1667 :
1668 : // Initialize the svc according to current VisSet
1669 : // (this counts intervals, sizes CalSet)
1670 : Vector<Int> nChunkPerSol;
1671 : Int nSol = sizeUpSolve(vs,nChunkPerSol);
1672 :
1673 : // Create the Caltable
1674 : createMemCalTable();
1675 :
1676 : // The iterator, VisBuffer
1677 : VisIter& vi(vs.iter());
1678 : VisBuffer vb(vi);
1679 :
1680 : // cout << "nSol = " << nSol << endl;
1681 : // cout << "nChunkPerSol = " << nChunkPerSol << endl;
1682 :
1683 : Vector<Int> slotidx(vs.numberSpw(),-1);
1684 :
1685 : Int nGood(0);
1686 : vi.originChunks();
1687 : for (Int isol=0;isol<nSol && vi.moreChunks();++isol) {
1688 :
1689 : // Arrange to accumulate
1690 : VisBuffAccumulator vba(nAnt(),preavg(),false);
1691 :
1692 : for (Int ichunk=0;ichunk<nChunkPerSol(isol);++ichunk) {
1693 :
1694 : // Current _chunk_'s spw
1695 : Int spw(vi.spectralWindow());
1696 :
1697 : // Abort if we encounter a spw for which a priori cal not available
1698 : if (!ve.spwOK(spw))
1699 : throw(AipsError("Pre-applied calibration not available for at least 1 spw. Check spw selection carefully."));
1700 :
1701 :
1702 : // Collapse each timestamp in this chunk according to VisEq
1703 : // with calibration and averaging
1704 : for (vi.origin(); vi.more(); vi++) {
1705 :
1706 : // Force read of the field Id
1707 : vb.fieldId();
1708 :
1709 : // Apply the channel mask
1710 : this->applyChanMask(vb);
1711 :
1712 : // This forces the data/model/wt I/O, and applies
1713 : // any prior calibrations
1714 : ve.collapse(vb);
1715 :
1716 : // If permitted/required by solvable component, normalize
1717 : if (normalizable())
1718 : vb.normalize();
1719 :
1720 : // If this solve not freqdep, and channels not averaged yet, do so
1721 : if (!freqDepMat() && vb.nChannel()>1)
1722 : vb.freqAveCubes();
1723 :
1724 : // Accumulate collapsed vb in a time average
1725 : vba.accumulate(vb);
1726 : }
1727 : // Advance the VisIter, if possible
1728 : if (vi.moreChunks()) vi.nextChunk();
1729 :
1730 : }
1731 :
1732 : // Finalize the averged VisBuffer
1733 : vba.finalizeAverage();
1734 :
1735 : // The VisBuffer to solve with
1736 : VisBuffer& svb(vba.aveVisBuff());
1737 :
1738 : // Make data amp- or phase-only
1739 : enforceAPonData(svb);
1740 :
1741 : // Establish meta-data for this interval
1742 : // (some of this may be used _during_ solve)
1743 : // (this sets currSpw() in the SVC)
1744 : Bool vbOk=syncSolveMeta(svb,-1);
1745 :
1746 : Int thisSpw=spwMap()(svb.spectralWindow());
1747 : slotidx(thisSpw)++;
1748 :
1749 : // We are actually solving for all channels simultaneously
1750 : solveCPar().reference(solveAllCPar());
1751 : solveParOK().reference(solveAllParOK());
1752 : solveParErr().reference(solveAllParErr());
1753 : solveParSNR().reference(solveAllParSNR());
1754 :
1755 : // Fill solveCPar() with 1, nominally, and flagged
1756 : solveCPar()=Complex(1.0);
1757 : solveParOK()=false;
1758 :
1759 : if (vbOk && svb.nRow()>0) {
1760 :
1761 : // Insist that channel,row shapes match
1762 : IPosition visshape(svb.visCube().shape());
1763 : AlwaysAssert(solveCPar().shape().getLast(2)==visshape.getLast(2),AipsError);
1764 :
1765 : // Zero flagged data
1766 : IPosition vblc(3,0,0,0);
1767 : IPosition vtrc(visshape); vtrc-=1;
1768 : Int nCorr(visshape(0));
1769 : for (Int i=0;i<nCorr;++i) {
1770 : vblc(0)=vtrc(0)=i;
1771 : svb.visCube()(vblc,vtrc).reform(visshape.getLast(2))(svb.flag())=Complex(1.0);
1772 : }
1773 :
1774 : // Form correct slice of visCube to copy to solveCPar
1775 : IPosition vcblc(3,0,0,0);
1776 : IPosition vctrc(svb.visCube().shape()); vctrc-=1;
1777 : IPosition vcstr(3,1,1,1);
1778 :
1779 : IPosition spblc(3,0,0,0);
1780 : IPosition sptrc(solveCPar().shape()); sptrc-=1;
1781 :
1782 : IPosition flshape(svb.flag().shape());
1783 :
1784 : switch (nCorr) {
1785 : case 1: {
1786 : // fill 1st par only
1787 : spblc(0)=sptrc(0)=0;
1788 : solveCPar()(spblc,sptrc)=svb.visCube();
1789 : // first pol flags
1790 : solveParOK()(spblc,sptrc).reform(flshape)=!svb.flag();
1791 : break;
1792 : }
1793 : case 2: {
1794 : // shapes match
1795 : solveCPar()=svb.visCube();
1796 : spblc(0)=sptrc(0)=0;
1797 : solveParOK()(spblc,sptrc).reform(flshape)=!svb.flag();
1798 : spblc(0)=sptrc(0)=1;
1799 : solveParOK()(spblc,sptrc).reform(flshape)=!svb.flag();
1800 :
1801 : break;
1802 : }
1803 : case 4: {
1804 : // Slice visCube with stride
1805 : vcstr(0)=3;
1806 : solveCPar()=svb.visCube()(vcblc,vctrc,vcstr);
1807 : spblc(0)=sptrc(0)=0;
1808 : solveParOK()(spblc,sptrc).reform(flshape)=!svb.flag();
1809 : spblc(0)=sptrc(0)=1;
1810 : solveParOK()(spblc,sptrc).reform(flshape)=!svb.flag();
1811 :
1812 : break;
1813 : }
1814 : default:
1815 : throw(AipsError("Problem in MMueller::selfSolve."));
1816 : break;
1817 : }
1818 :
1819 : nGood++;
1820 :
1821 : // record in in-memory caltable
1822 : keepNCT();
1823 : }
1824 : }
1825 :
1826 : logSink() << " Found good "
1827 : << typeName() << " solutions in "
1828 : << nGood << " intervals."
1829 : << LogIO::POST;
1830 :
1831 : // Store whole of result in a caltable
1832 : if (nGood==0)
1833 : logSink() << "No output calibration table written."
1834 : << LogIO::POST;
1835 : else {
1836 :
1837 : // Do global post-solve tinkering (e.g., normalization)
1838 : globalPostSolveTinker();
1839 :
1840 : // write the table
1841 : storeNCT();
1842 : }
1843 :
1844 : }
1845 :
1846 : void MMueller::globalPostSolveTinker() {
1847 :
1848 : // normalize, if requested
1849 : if (solnorm()) normalize();
1850 :
1851 : }
1852 :
1853 : void MMueller::createCorruptor(const VisIter& vi, const Record& simpar, const Int nSim)
1854 : {
1855 : LogIO os(LogOrigin("MM", "createCorruptor()", WHERE));
1856 :
1857 : if (prtlev()>2) cout << " MM::createCorruptor()" << endl;
1858 : os << LogIO::DEBUG1 << " MM::createCorruptor()"
1859 : << LogIO::POST;
1860 :
1861 : atmcorruptor_p = new AtmosCorruptor();
1862 : corruptor_p = atmcorruptor_p;
1863 :
1864 : // call generic parent to set corr,spw,etc info
1865 : SolvableVisCal::createCorruptor(vi,simpar,nSim);
1866 :
1867 : Int rxType(0); // 0=2SB, 1=DSB
1868 : if (simpar.isDefined("rxType")) {
1869 : rxType=simpar.asInt("rxType");
1870 : }
1871 :
1872 : // this is the M type corruptor - maybe we should make the corruptor
1873 : // take the VC as an argument
1874 : atmcorruptor_p->initialize(vi,simpar,VisCal::M,rxType);
1875 : }
1876 :
1877 :
1878 :
1879 :
1880 :
1881 :
1882 : // **********************************************************
1883 : // MfMueller: freq-dep MMueller
1884 : //
1885 :
1886 : MfMueller::MfMueller(VisSet& vs) :
1887 : VisCal(vs), // virtual base
1888 : VisMueller(vs), // virtual base
1889 : MMueller(vs) // immediate parent
1890 : {
1891 : if (prtlev()>2) cout << "Mf::Mf(vs)" << endl;
1892 : }
1893 :
1894 : MfMueller::MfMueller(String msname,Int MSnAnt,Int MSnSpw) :
1895 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
1896 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
1897 : MMueller(msname,MSnAnt,MSnSpw) // immediate parent
1898 : {
1899 : if (prtlev()>2) cout << "Mf::Mf(msname,MSnAnt,MSnSpw)" << endl;
1900 : }
1901 :
1902 : MfMueller::MfMueller(const MSMetaInfoForCal& msmc) :
1903 : VisCal(msmc), // virtual base
1904 : VisMueller(msmc), // virtual base
1905 : MMueller(msmc) // immediate parent
1906 : {
1907 : if (prtlev()>2) cout << "Mf::Mf(msmc)" << endl;
1908 : }
1909 :
1910 : MfMueller::MfMueller(const Int& nAnt) :
1911 : VisCal(nAnt),
1912 : VisMueller(nAnt),
1913 : MMueller(nAnt)
1914 : {
1915 : if (prtlev()>2) cout << "Mf::Mf(nAnt)" << endl;
1916 : }
1917 :
1918 : MfMueller::~MfMueller() {
1919 : if (prtlev()>2) cout << "Mf::~Mf()" << endl;
1920 : }
1921 :
1922 : void MfMueller::normalize() {
1923 :
1924 : // This is just like BJones
1925 : // TBD: consolidate (via generalized SVC::normalize(Block<String> cols) )
1926 :
1927 : // Only if we have a CalTable, and it is not empty
1928 : if (ct_ && ct_->nrow()>0) {
1929 :
1930 : // TBD: trap attempts to normalize a caltable containing FPARAM (non-Complex)?
1931 :
1932 : logSink() << "Normalizing solutions per spw, pol, baseline, time"
1933 : << LogIO::POST;
1934 :
1935 : Block<String> col(4);
1936 : col[0]="SPECTRAL_WINDOW_ID";
1937 : col[1]="TIME";
1938 : col[2]="ANTENNA1";
1939 : col[3]="ANTENNA2";
1940 : CTIter ctiter(*ct_,col);
1941 :
1942 : // Cube iteration axes are pol and ant
1943 : IPosition itax(2,0,2);
1944 :
1945 : while (!ctiter.pastEnd()) {
1946 : Cube<Bool> fl(ctiter.flag());
1947 :
1948 : if (nfalse(fl)>0) {
1949 : Cube<Complex> p(ctiter.cparam());
1950 : ArrayIterator<Complex> soliter(p,itax,false);
1951 : ArrayIterator<Bool> fliter(fl,itax,false);
1952 : while (!soliter.pastEnd()) {
1953 : normSolnArray(soliter.array(),!fliter.array(),true); // Do phase
1954 : soliter.next();
1955 : fliter.next();
1956 : }
1957 :
1958 : // record result...
1959 : ctiter.setcparam(p);
1960 : }
1961 : ctiter.next();
1962 : }
1963 : }
1964 : cout << "End of MfMueller::normalize()" << endl;
1965 : }
1966 :
1967 : */
1968 :
1969 : // **********************************************************
1970 : // TOpac
1971 : //
1972 :
1973 0 : TOpac::TOpac(VisSet& vs) :
1974 : VisCal(vs),
1975 : VisMueller(vs),
1976 : TJones(vs),
1977 0 : za_()
1978 : {
1979 0 : if (prtlev()>2) cout << "TOpac::TOpac(vs)" << endl;
1980 0 : }
1981 :
1982 0 : TOpac::TOpac(String msname,Int MSnAnt,Int MSnSpw) :
1983 : VisCal(msname,MSnAnt,MSnSpw),
1984 : VisMueller(msname,MSnAnt,MSnSpw),
1985 : TJones(msname,MSnAnt,MSnSpw),
1986 0 : za_()
1987 : {
1988 0 : if (prtlev()>2) cout << "TOpac::TOpac(msname,MSnAnt,MSnSpw)" << endl;
1989 0 : }
1990 :
1991 0 : TOpac::TOpac(const MSMetaInfoForCal& msmc) :
1992 : VisCal(msmc),
1993 : VisMueller(msmc),
1994 : TJones(msmc),
1995 0 : za_()
1996 : {
1997 0 : if (prtlev()>2) cout << "TOpac::TOpac(msmc)" << endl;
1998 0 : }
1999 :
2000 0 : TOpac::~TOpac() {
2001 0 : if (prtlev()>2) cout << "TOpac::~TOpac()" << endl;
2002 0 : }
2003 :
2004 0 : void TOpac::setApply(const Record& applypar) {
2005 :
2006 : // TBD: Handle spwmap properly?
2007 :
2008 : // Prepare zenith angle storage
2009 0 : za().resize(nAnt());
2010 0 : za().set(0.0);
2011 :
2012 0 : String table("");
2013 0 : if (applypar.isDefined("table"))
2014 0 : table=applypar.asString("table");
2015 :
2016 0 : if (table!="")
2017 : // Attempt new-fashioned opacity table
2018 0 : SolvableVisCal::setApply(applypar);
2019 : else {
2020 :
2021 : // We are applying
2022 0 : setSolved(false);
2023 0 : setApplied(true);
2024 :
2025 0 : LogMessage message;
2026 0 : { ostringstream o;
2027 0 : o<< "Invoking opacity application without a caltable (e.g., " << endl
2028 0 : << " via opacity=[...] in calibration tasks) will soon be disabled." << endl
2029 0 : << " Please begin using gencal to generate an opacity caltable, " << endl
2030 0 : << " and then supply that table in the standard manner." << endl;
2031 0 : message.message(o);
2032 0 : message.priority(LogMessage::WARN);
2033 0 : logSink().post(message);
2034 : }
2035 :
2036 : // Detect and extract opacities from applypar record (old way)
2037 0 : if (applypar.isDefined("opacity")) {
2038 0 : opacity_.resize();
2039 0 : opacity_=applypar.asArrayDouble("opacity");
2040 : }
2041 0 : Int nopac=opacity_.nelements();
2042 :
2043 0 : if (nopac>0 && sum(opacity_)>0.0) {
2044 :
2045 : // Old-fashioned opacity_ is non-trivial, so adopt them
2046 :
2047 0 : if (nopac<nSpw()) {
2048 : // Resize (with copy) to match nSpw,
2049 : // duplicating last specified entry
2050 0 : opacity_.resize(nSpw(),true);
2051 0 : opacity_(IPosition(1,nopac),IPosition(1,nSpw()-1))=opacity_(nopac-1);
2052 : }
2053 :
2054 0 : Int oldspw; oldspw=currSpw();
2055 0 : for (Int ispw=0;ispw<nSpw();++ispw) {
2056 0 : currSpw()=ispw;
2057 0 : currRPar().resize(1,1,nAnt());
2058 0 : currRPar()=Float(opacity_(ispw));
2059 0 : currParOK().resize(1,1,nAnt());
2060 0 : currParOK()=true;
2061 : }
2062 0 : currSpw()=oldspw;
2063 :
2064 : }
2065 : else
2066 0 : throw(AipsError("No opacity info specified."));
2067 : }
2068 :
2069 0 : }
2070 :
2071 0 : String TOpac::applyinfo() {
2072 :
2073 0 : if (opacity_.nelements()==0)
2074 0 : return TJones::applyinfo();
2075 : else {
2076 0 : ostringstream o;
2077 0 : o << typeName();
2078 0 : o << ": opacity=" << opacity_;
2079 0 : o << boolalpha
2080 0 : << " calWt=" << calWt();
2081 0 : return String(o);
2082 : }
2083 : }
2084 :
2085 : // TOpac needs zenith angle (old VB version)
2086 0 : void TOpac::syncMeta(const VisBuffer& vb) {
2087 :
2088 : // Call parent (sets currTime())
2089 0 : TJones::syncMeta(vb);
2090 :
2091 : // Current time's zenith angle...
2092 0 : za().resize(nAnt());
2093 0 : Vector<MDirection> antazel(vb.azel(currTime()));
2094 0 : Double* a=za().data();
2095 0 : for (Int iant=0;iant<nAnt();++iant,++a)
2096 0 : (*a)=C::pi_2 - antazel(iant).getAngle().getValue()(1);
2097 :
2098 0 : }
2099 :
2100 : // TOpac needs zenith angle (VB2 version)
2101 0 : void TOpac::syncMeta2(const vi::VisBuffer2& vb) {
2102 :
2103 : // Call parent (sets currTime())
2104 0 : TJones::syncMeta2(vb);
2105 :
2106 : // Current time's zenith angle...
2107 0 : za().resize(nAnt());
2108 0 : Vector<MDirection> antazel(vb.azel(currTime()));
2109 0 : Double* a=za().data();
2110 0 : for (Int iant=0;iant<nAnt();++iant,++a)
2111 0 : (*a)=C::pi_2 - antazel(iant).getAngle().getValue()(1);
2112 :
2113 0 : }
2114 :
2115 :
2116 0 : void TOpac::calcPar() {
2117 :
2118 0 : if (prtlev()>6) cout << " TOpac::calcPar()" << endl;
2119 :
2120 : // If we are interpolating from a table, get opacity(time)
2121 0 : if (ci_ || cpp_)
2122 0 : SolvableVisCal::calcPar();
2123 : else
2124 0 : throw(AipsError("Error in TOpac::calcPar()"));
2125 :
2126 : // Pars now valid, matrices not yet
2127 0 : validateP();
2128 0 : invalidateJ(); // Force new calculation of za-dep matrix elements
2129 :
2130 0 : }
2131 :
2132 :
2133 0 : void TOpac::calcAllJones() {
2134 :
2135 0 : if (prtlev()>6) cout << " TOpac::calcAllJones()" << endl;
2136 :
2137 : // Nominally no opacity
2138 0 : currJElem()=Complex(1.0);
2139 0 : currJElemOK()=currParOK();
2140 :
2141 0 : Complex* J=currJElem().data();
2142 0 : Float* op=currRPar().data();
2143 0 : Bool* opok=currParOK().data();
2144 0 : Double* a=za().data();
2145 0 : for (Int iant=0; iant<nAnt(); ++iant,++J,++op,++opok,++a) {
2146 0 : if ((*opok) && (*a)<C::pi_2)
2147 0 : (*J) = Complex(sqrt(exp(-1.0 * Double(*op)/cos(*a))));
2148 : }
2149 :
2150 0 : }
2151 :
2152 :
2153 : // **********************************************************
2154 : // TfOpac
2155 : //
2156 :
2157 0 : TfOpac::TfOpac(VisSet& vs) :
2158 : VisCal(vs), // virtual base
2159 : VisMueller(vs), // virtual base
2160 0 : TOpac(vs) // immediate parent
2161 : {
2162 0 : if (prtlev()>2) cout << "TfOpac::TfOpac(vs)" << endl;
2163 0 : }
2164 :
2165 0 : TfOpac::TfOpac(String msname,Int MSnAnt,Int MSnSpw) :
2166 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
2167 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
2168 0 : TOpac(msname,MSnAnt,MSnSpw) // immediate parent
2169 : {
2170 0 : if (prtlev()>2) cout << "TfOpac::TfOpac(msname,MSnAnt,MSnSpw)" << endl;
2171 0 : }
2172 :
2173 0 : TfOpac::TfOpac(const MSMetaInfoForCal& msmc) :
2174 : VisCal(msmc), // virtual base
2175 : VisMueller(msmc), // virtual base
2176 0 : TOpac(msmc) // immediate parent
2177 : {
2178 0 : if (prtlev()>2) cout << "TfOpac::TfOpac(msmc)" << endl;
2179 0 : }
2180 :
2181 :
2182 : //TfOpac::TfOpac(const Int& nAnt) :
2183 : // VisCal(nAnt),
2184 : // VisMueller(nAnt),
2185 : // TOpac(nAnt)
2186 : //{
2187 : // if (prtlev()>2) cout << "TfOpac::TfOpac(nAnt)" << endl;
2188 : //}
2189 :
2190 0 : TfOpac::~TfOpac() {
2191 0 : if (prtlev()>2) cout << "TfOpac::~TfOpac()" << endl;
2192 0 : }
2193 :
2194 0 : void TfOpac::calcAllJones() {
2195 :
2196 0 : if (prtlev()>6) cout << " TfOpac::calcAllJones()" << endl;
2197 :
2198 : // Nominally no opacity
2199 0 : currJElem()=Complex(1.0);
2200 0 : currJElemOK()=currParOK();
2201 :
2202 0 : Complex* J=currJElem().data();
2203 0 : Float* op=currRPar().data();
2204 0 : Bool* opok=currParOK().data();
2205 0 : Double* a=za().data();
2206 0 : for (Int iant=0; iant<nAnt(); ++iant,++a) {
2207 0 : for (Int ich=0; ich<nChanMat(); ich++, J++, op++, opok++) {
2208 0 : if ((*opok) && (*a)<C::pi_2)
2209 0 : (*J) = Complex(sqrt(exp(-1.0 * Double(*op)/cos(*a))));
2210 : }
2211 : }
2212 0 : }
2213 :
2214 :
2215 0 : void TfOpac::calcWtScale() {
2216 :
2217 : // Initialize - not used for TfOpac, but we need to overwrite the
2218 : // single channel version or it will throw an exception
2219 0 : currWtScale()=1.0;
2220 :
2221 0 : }
2222 :
2223 :
2224 : // **********************************************************
2225 : // MMueller: baseline-based (closure) solution
2226 : //
2227 :
2228 0 : MMueller::MMueller(VisSet& vs) :
2229 : VisCal(vs), // virtual base
2230 : VisMueller(vs), // virtual base
2231 : SolvableVisMueller(vs), // immediate parent
2232 0 : useGenGathSolve_p(false) // VisSet-driven ctor
2233 : {
2234 0 : if (prtlev()>2) cout << "M::M(vs)" << endl;
2235 0 : }
2236 :
2237 0 : MMueller::MMueller(String msname,Int MSnAnt,Int MSnSpw) :
2238 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
2239 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
2240 : SolvableVisMueller(msname,MSnAnt,MSnSpw), // immediate parent
2241 0 : useGenGathSolve_p(true) // modern ctor
2242 : {
2243 0 : if (prtlev()>2) cout << "M::M(msname,MSnAnt,MSnSpw)" << endl;
2244 0 : }
2245 :
2246 0 : MMueller::MMueller(const MSMetaInfoForCal& msmc) :
2247 : VisCal(msmc), // virtual base
2248 : VisMueller(msmc), // virtual base
2249 : SolvableVisMueller(msmc), // immediate parent
2250 0 : useGenGathSolve_p(true) // modern ctor
2251 : {
2252 0 : if (prtlev()>2) cout << "M::M(msmc)" << endl;
2253 0 : }
2254 :
2255 0 : MMueller::MMueller(const Int& nAnt) :
2256 : VisCal(nAnt),
2257 : VisMueller(nAnt),
2258 : SolvableVisMueller(nAnt),
2259 0 : useGenGathSolve_p(true) // modern ctor
2260 : {
2261 0 : if (prtlev()>2) cout << "M::M(nAnt)" << endl;
2262 0 : }
2263 :
2264 0 : MMueller::~MMueller() {
2265 0 : if (prtlev()>2) cout << "M::~M()" << endl;
2266 0 : }
2267 :
2268 0 : void MMueller::setApply(const Record& apply) {
2269 :
2270 0 : SolvableVisCal::setApply(apply);
2271 :
2272 : // Force calwt to false for now
2273 0 : calWt()=false;
2274 :
2275 0 : }
2276 :
2277 0 : void MMueller::newselfSolve(VisSet& vs, VisEquation& ve) {
2278 :
2279 0 : if (prtlev()>4) cout << " M::selfSolve(ve)" << endl;
2280 :
2281 0 : AlwaysAssert(!useGenGathSolve_p,AipsError);
2282 :
2283 : // Inform logger/history
2284 0 : logSink() << "Solving for " << typeName()
2285 0 : << LogIO::POST;
2286 :
2287 : // Initialize the svc according to current VisSet
2288 : // (this counts intervals, sizes CalSet)
2289 0 : Vector<Int> nChunkPerSol;
2290 0 : Int nSol = sizeUpSolve(vs,nChunkPerSol);
2291 :
2292 : // Create the Caltable
2293 0 : createMemCalTable();
2294 :
2295 : // The iterator, VisBuffer
2296 0 : VisIter& vi(vs.iter());
2297 0 : VisBuffer vb(vi);
2298 :
2299 : // cout << "nSol = " << nSol << endl;
2300 : // cout << "nChunkPerSol = " << nChunkPerSol << endl;
2301 :
2302 0 : Vector<Int> slotidx(vs.numberSpw(),-1);
2303 :
2304 0 : Int nGood(0);
2305 0 : vi.originChunks();
2306 0 : for (Int isol=0;isol<nSol && vi.moreChunks();++isol) {
2307 :
2308 : // Arrange to accumulate
2309 0 : VisBuffAccumulator vba(nAnt(),preavg(),false);
2310 :
2311 0 : for (Int ichunk=0;ichunk<nChunkPerSol(isol);++ichunk) {
2312 :
2313 : // Current _chunk_'s spw
2314 0 : Int spw(vi.spectralWindow());
2315 :
2316 : // Abort if we encounter a spw for which a priori cal not available
2317 0 : if (!ve.spwOK(spw))
2318 0 : throw(AipsError("Pre-applied calibration not available for at least 1 spw. Check spw selection carefully."));
2319 :
2320 :
2321 : // Collapse each timestamp in this chunk according to VisEq
2322 : // with calibration and averaging
2323 0 : for (vi.origin(); vi.more(); vi++) {
2324 :
2325 : // Force read of the field Id
2326 0 : vb.fieldId();
2327 :
2328 : // Apply the channel mask
2329 0 : this->applyChanMask(vb);
2330 :
2331 : // This forces the data/model/wt I/O, and applies
2332 : // any prior calibrations
2333 0 : ve.collapse(vb);
2334 :
2335 : // If permitted/required by solvable component, normalize
2336 0 : if (normalizable())
2337 0 : vb.normalize();
2338 :
2339 : // If this solve not freqdep, and channels not averaged yet, do so
2340 0 : if (!freqDepMat() && vb.nChannel()>1)
2341 0 : vb.freqAveCubes();
2342 :
2343 : // Accumulate collapsed vb in a time average
2344 0 : vba.accumulate(vb);
2345 : }
2346 : // Advance the VisIter, if possible
2347 0 : if (vi.moreChunks()) vi.nextChunk();
2348 :
2349 : }
2350 :
2351 : // Finalize the averged VisBuffer
2352 0 : vba.finalizeAverage();
2353 :
2354 : // The VisBuffer to solve with
2355 0 : VisBuffer& svb(vba.aveVisBuff());
2356 :
2357 : // Make data amp- or phase-only
2358 0 : enforceAPonData(svb);
2359 :
2360 : // Establish meta-data for this interval
2361 : // (some of this may be used _during_ solve)
2362 : // (this sets currSpw() in the SVC)
2363 0 : Bool vbOk=syncSolveMeta(svb,-1);
2364 :
2365 0 : Int thisSpw=spwMap()(svb.spectralWindow());
2366 0 : slotidx(thisSpw)++;
2367 :
2368 : // We are actually solving for all channels simultaneously
2369 0 : solveCPar().reference(solveAllCPar());
2370 0 : solveParOK().reference(solveAllParOK());
2371 0 : solveParErr().reference(solveAllParErr());
2372 0 : solveParSNR().reference(solveAllParSNR());
2373 :
2374 : // Fill solveCPar() with 1, nominally, and flagged
2375 0 : solveCPar()=Complex(1.0);
2376 0 : solveParOK()=false;
2377 :
2378 0 : if (vbOk && svb.nRow()>0) {
2379 :
2380 : // Insist that channel,row shapes match
2381 0 : IPosition visshape(svb.visCube().shape());
2382 0 : AlwaysAssert(solveCPar().shape().getLast(2)==visshape.getLast(2),AipsError);
2383 :
2384 : // Zero flagged data
2385 0 : IPosition vblc(3,0,0,0);
2386 0 : IPosition vtrc(visshape); vtrc-=1;
2387 0 : Int nCorr(visshape(0));
2388 0 : for (Int i=0;i<nCorr;++i) {
2389 0 : vblc(0)=vtrc(0)=i;
2390 0 : svb.visCube()(vblc,vtrc).reform(visshape.getLast(2))(svb.flag())=Complex(1.0);
2391 : }
2392 :
2393 : // Form correct slice of visCube to copy to solveCPar
2394 0 : IPosition vcblc(3,0,0,0);
2395 0 : IPosition vctrc(svb.visCube().shape()); vctrc-=1;
2396 0 : IPosition vcstr(3,1,1,1);
2397 :
2398 0 : IPosition spblc(3,0,0,0);
2399 0 : IPosition sptrc(solveCPar().shape()); sptrc-=1;
2400 :
2401 0 : IPosition flshape(svb.flag().shape());
2402 :
2403 0 : switch (nCorr) {
2404 0 : case 1: {
2405 : // fill 1st par only
2406 0 : spblc(0)=sptrc(0)=0;
2407 0 : solveCPar()(spblc,sptrc)=svb.visCube();
2408 : // first pol flags
2409 0 : solveParOK()(spblc,sptrc).reform(flshape)=!svb.flag();
2410 0 : break;
2411 : }
2412 0 : case 2: {
2413 : // shapes match
2414 0 : solveCPar()=svb.visCube();
2415 0 : spblc(0)=sptrc(0)=0;
2416 0 : solveParOK()(spblc,sptrc).reform(flshape)=!svb.flag();
2417 0 : spblc(0)=sptrc(0)=1;
2418 0 : solveParOK()(spblc,sptrc).reform(flshape)=!svb.flag();
2419 :
2420 0 : break;
2421 : }
2422 0 : case 4: {
2423 : // Slice visCube with stride
2424 0 : vcstr(0)=3;
2425 0 : solveCPar()=svb.visCube()(vcblc,vctrc,vcstr);
2426 0 : spblc(0)=sptrc(0)=0;
2427 0 : solveParOK()(spblc,sptrc).reform(flshape)=!svb.flag();
2428 0 : spblc(0)=sptrc(0)=1;
2429 0 : solveParOK()(spblc,sptrc).reform(flshape)=!svb.flag();
2430 :
2431 0 : break;
2432 : }
2433 0 : default:
2434 0 : throw(AipsError("Problem in MMueller::selfSolve."));
2435 : break;
2436 : }
2437 :
2438 0 : nGood++;
2439 :
2440 : // record in in-memory caltable
2441 0 : keepNCT();
2442 : }
2443 : }
2444 :
2445 : logSink() << " Found good "
2446 0 : << typeName() << " solutions in "
2447 : << nGood << " intervals."
2448 0 : << LogIO::POST;
2449 :
2450 : // Store whole of result in a caltable
2451 0 : if (nGood==0)
2452 : logSink() << "No output calibration table written."
2453 0 : << LogIO::POST;
2454 : else {
2455 :
2456 : // Do global post-solve tinkering (e.g., normalization)
2457 0 : globalPostSolveTinker();
2458 :
2459 : // write the table
2460 0 : storeNCT();
2461 : }
2462 :
2463 0 : }
2464 :
2465 0 : void MMueller::globalPostSolveTinker() {
2466 :
2467 : // normalize, if requested
2468 0 : if (solnorm()) normalize();
2469 :
2470 0 : }
2471 :
2472 0 : void MMueller::selfSolveOne(SDBList& sdbs) {
2473 :
2474 0 : AlwaysAssert(useGenGathSolve_p,AipsError);
2475 :
2476 : // Call the implementation...
2477 0 : this->solveOne(sdbs);
2478 :
2479 0 : }
2480 :
2481 0 : void MMueller::createCorruptor(const VisIter& vi, const Record& simpar, const Int nSim)
2482 : {
2483 0 : LogIO os(LogOrigin("MM", "createCorruptor()", WHERE));
2484 :
2485 0 : if (prtlev()>2) cout << " MM::createCorruptor()" << endl;
2486 : os << LogIO::DEBUG1 << " MM::createCorruptor()"
2487 0 : << LogIO::POST;
2488 :
2489 0 : atmcorruptor_p = new AtmosCorruptor();
2490 0 : corruptor_p = atmcorruptor_p;
2491 :
2492 : // call generic parent to set corr,spw,etc info
2493 0 : SolvableVisCal::createCorruptor(vi,simpar,nSim);
2494 :
2495 0 : Int rxType(0); // 0=2SB, 1=DSB
2496 0 : if (simpar.isDefined("rxType")) {
2497 0 : rxType=simpar.asInt("rxType");
2498 : }
2499 :
2500 : // this is the M type corruptor - maybe we should make the corruptor
2501 : // take the VC as an argument
2502 0 : atmcorruptor_p->initialize(vi,simpar,VisCal::M,rxType);
2503 0 : }
2504 :
2505 0 : void MMueller::solveOne(SDBList& sdbs) {
2506 :
2507 0 : AlwaysAssert(useGenGathSolve_p,AipsError);
2508 :
2509 : // This just a simple average of the parallel-hand
2510 :
2511 0 : Int nSDB=sdbs.nSDB();
2512 :
2513 : //cout << "nSDB=" << nSDB << endl;
2514 :
2515 : // We are actually solving for all channels simultaneously
2516 0 : solveCPar().reference(solveAllCPar());
2517 0 : solveParOK().reference(solveAllParOK());
2518 0 : solveParErr().reference(solveAllParErr());
2519 0 : solveParSNR().reference(solveAllParSNR());
2520 :
2521 : // Fill solveCPar() with 0, nominally, and flagged
2522 0 : solveCPar()=Complex(0.0);
2523 0 : solveParOK()=false;
2524 :
2525 0 : Cube<Float> sumwt(solveCPar().shape(),0.0f);
2526 0 : Int nChan=sdbs.nChannels();
2527 0 : AlwaysAssert(nChan==nChanPar(),AipsError); // channelization should be consistent
2528 0 : Int nCorr=sdbs.nCorrelations();
2529 0 : Int dCorr=( nCorr==4 ? 3 : 1 );
2530 :
2531 0 : for (Int isdb=0;isdb<nSDB;++isdb) {
2532 0 : SolveDataBuffer& sdb(sdbs(isdb));
2533 :
2534 0 : Int nRow=sdb.nRows();
2535 :
2536 : // Zero flagged data
2537 0 : Cube<Complex> vc(sdb.visCubeCorrected()); // non-const access
2538 0 : vc(sdb.flagCube())=Complex(0.0);
2539 :
2540 0 : const Vector<Int>& a1(sdb.antenna1()), a2(sdb.antenna2());
2541 :
2542 0 : Cube<Complex> vis,sol;
2543 0 : Cube<Float> wt,swt;
2544 0 : Cube<Bool> fl,solok;
2545 0 : Slice chsl; // whole slice on chan axis
2546 0 : for (Int irow=0;irow<nRow;++irow) {
2547 0 : Slice vrsl(irow,1,1), srsl(blnidx(a1(irow),a2(irow)),1,1);
2548 0 : for (Int icor=0;icor<nCorr;icor+=dCorr) {
2549 0 : Slice vcsl(icor,1,1), scsl(icor%2,1,1);
2550 :
2551 0 : vis.reference(sdb.visCubeCorrected()(vcsl,chsl,vrsl));
2552 0 : wt.reference(sdb.weightSpectrum()(vcsl,chsl,vrsl));
2553 0 : fl.reference(sdb.flagCube()(vcsl,chsl,vrsl));
2554 0 : sol.reference(solveCPar()(scsl,chsl,srsl));
2555 0 : solok.reference(solveParOK()(scsl,chsl,srsl));
2556 0 : swt.reference(sumwt(scsl,chsl,srsl));
2557 :
2558 : // Accumulate
2559 0 : sol+=(vis*wt);
2560 0 : swt+=wt;
2561 0 : solok|=!fl;
2562 : }
2563 : }
2564 : }
2565 :
2566 : // Normalize
2567 0 : sumwt(!solveParOK())=(1.0);
2568 0 : solveCPar()/=sumwt;
2569 :
2570 : // TBD.... (move to MMueller::postSolveTinker, since this is wrong for AMueller!)
2571 : //solveCPar()(!solveParOK())=Complex(1.0); // set flagged ones to 1.0
2572 :
2573 0 : }
2574 :
2575 :
2576 : // **********************************************************
2577 : // MfMueller: freq-dep MMueller
2578 : //
2579 :
2580 0 : MfMueller::MfMueller(VisSet& vs) :
2581 : VisCal(vs), // virtual base
2582 : VisMueller(vs), // virtual base
2583 0 : MMueller(vs) // immediate parent
2584 : {
2585 0 : if (prtlev()>2) cout << "Mf::Mf(vs)" << endl;
2586 0 : }
2587 :
2588 0 : MfMueller::MfMueller(String msname,Int MSnAnt,Int MSnSpw) :
2589 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
2590 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
2591 0 : MMueller(msname,MSnAnt,MSnSpw) // immediate parent
2592 : {
2593 0 : if (prtlev()>2) cout << "Mf::Mf(msname,MSnAnt,MSnSpw)" << endl;
2594 0 : }
2595 :
2596 0 : MfMueller::MfMueller(const MSMetaInfoForCal& msmc) :
2597 : VisCal(msmc), // virtual base
2598 : VisMueller(msmc), // virtual base
2599 0 : MMueller(msmc) // immediate parent
2600 : {
2601 0 : if (prtlev()>2) cout << "Mf::Mf(msmc)" << endl;
2602 0 : }
2603 :
2604 0 : MfMueller::MfMueller(const Int& nAnt) :
2605 : VisCal(nAnt),
2606 : VisMueller(nAnt),
2607 0 : MMueller(nAnt)
2608 : {
2609 0 : if (prtlev()>2) cout << "Mf::Mf(nAnt)" << endl;
2610 0 : }
2611 :
2612 0 : MfMueller::~MfMueller() {
2613 0 : if (prtlev()>2) cout << "Mf::~Mf()" << endl;
2614 0 : }
2615 :
2616 0 : void MfMueller::normalize() {
2617 :
2618 : // This is just like BJones
2619 : // TBD: consolidate (via generalized SVC::normalize(Block<String> cols) )
2620 :
2621 : // Only if we have a CalTable, and it is not empty
2622 0 : if (ct_ && ct_->nrow()>0) {
2623 :
2624 : // TBD: trap attempts to normalize a caltable containing FPARAM (non-Complex)?
2625 :
2626 : logSink() << "Normalizing solutions per spw, pol, baseline, time"
2627 0 : << LogIO::POST;
2628 :
2629 0 : Block<String> col(4);
2630 0 : col[0]="SPECTRAL_WINDOW_ID";
2631 0 : col[1]="TIME";
2632 0 : col[2]="ANTENNA1";
2633 0 : col[3]="ANTENNA2";
2634 0 : CTIter ctiter(*ct_,col);
2635 :
2636 : // Cube iteration axes are pol and ant
2637 0 : IPosition itax(2,0,2);
2638 :
2639 0 : while (!ctiter.pastEnd()) {
2640 0 : Cube<Bool> fl(ctiter.flag());
2641 :
2642 0 : if (nfalse(fl)>0) {
2643 0 : Cube<Complex> p(ctiter.cparam());
2644 0 : ArrayIterator<Complex> soliter(p,itax,false);
2645 0 : ArrayIterator<Bool> fliter(fl,itax,false);
2646 0 : while (!soliter.pastEnd()) {
2647 0 : normSolnArray(soliter.array(),!fliter.array(),true); // Do phase
2648 0 : soliter.next();
2649 0 : fliter.next();
2650 : }
2651 :
2652 : // record result...
2653 0 : ctiter.setcparam(p);
2654 : }
2655 0 : ctiter.next();
2656 : }
2657 : }
2658 0 : cout << "End of MfMueller::normalize()" << endl;
2659 0 : }
2660 :
2661 : } //# NAMESPACE CASA - END
|
