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