Line data Source code
1 : //# DJones.cc: Implementation of polarization-related calibration 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/DJones.h>
28 : #include <synthesis/MeasurementComponents/CalCorruptor.h>
29 : #include <synthesis/MeasurementComponents/MSMetaInfoForCal.h>
30 : #include <synthesis/MeasurementComponents/SolveDataBuffer.h>
31 : #include <msvis/MSVis/VisBuffer.h>
32 : #include <msvis/MSVis/VisBuffAccumulator.h>
33 : #include <synthesis/CalTables/CTIter.h>
34 : #include <synthesis/MeasurementEquations/VisEquation.h>
35 : #include <casacore/scimath/Fitting/LSQFit.h>
36 : #include <casacore/scimath/Fitting/LinearFit.h>
37 : #include <casacore/scimath/Functionals/CompiledFunction.h>
38 : #include <casacore/scimath/Functionals/Polynomial.h>
39 : #include <casacore/scimath/Mathematics/AutoDiff.h>
40 : #include <casacore/casa/BasicMath/Math.h>
41 : #include <casacore/tables/TaQL/ExprNode.h>
42 :
43 : #include <casacore/casa/Arrays/ArrayMath.h>
44 : #include <casacore/casa/Arrays/MaskArrMath.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 : // DJones Implementations
70 : //
71 :
72 0 : DJones::DJones(VisSet& vs) :
73 : VisCal(vs), // virtual base
74 : VisMueller(vs), // virtual base
75 : SolvableVisJones(vs), // immediate parent
76 0 : solvePol_(0)
77 : {
78 0 : if (prtlev()>2) cout << "D::D(vs)" << endl;
79 :
80 0 : }
81 :
82 :
83 0 : DJones::DJones(String msname,Int MSnAnt,Int MSnSpw) :
84 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
85 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
86 : SolvableVisJones(msname,MSnAnt,MSnSpw), // immediate parent
87 0 : solvePol_(0)
88 : {
89 0 : if (prtlev()>2) cout << "D::D(msname,MSnAnt,MSnSpw)" << endl;
90 0 : }
91 :
92 0 : DJones::DJones(const MSMetaInfoForCal& msmc) :
93 : VisCal(msmc), // virtual base
94 : VisMueller(msmc), // virtual base
95 : SolvableVisJones(msmc), // immediate parent
96 0 : solvePol_(0)
97 : {
98 0 : if (prtlev()>2) cout << "D::D(msmc)" << endl;
99 0 : }
100 :
101 :
102 0 : DJones::DJones(const Int& nAnt) :
103 : VisCal(nAnt),
104 : VisMueller(nAnt),
105 : SolvableVisJones(nAnt),
106 0 : solvePol_(0)
107 : {
108 0 : if (prtlev()>2) cout << "D::D(nAnt)" << endl;
109 :
110 0 : }
111 :
112 0 : DJones::~DJones() {
113 0 : if (prtlev()>2) cout << "D::~D()" << endl;
114 0 : }
115 :
116 0 : void DJones::setApply(const Record& apply) {
117 :
118 0 : SolvableVisJones::setApply(apply);
119 :
120 : // Force calwt to false for now
121 0 : calWt()=false;
122 :
123 0 : }
124 :
125 :
126 0 : void DJones::setSolve(const Record& solvepar) {
127 :
128 : // Call parent
129 0 : SolvableVisJones::setSolve(solvepar);
130 :
131 : // Determine if we are solving for source pol or not
132 0 : if (solvepar.isDefined("type")) {
133 0 : String type = solvepar.asString("type");
134 0 : if (type.contains("QU")) {
135 0 : solvePol_=2;
136 0 : logSink() << "Will solve for source polarization (Q,U)" << LogIO::POST;
137 : }
138 0 : else if (type.contains("X")) {
139 0 : solvePol_=1;
140 0 : logSink() << "Will solve for source polarization position angle correction" << LogIO::POST;
141 : }
142 : else
143 0 : solvePol_=0;
144 0 : }
145 :
146 0 : logSink() << "Using only cross-hand data for instrumental polarization solution." << LogIO::POST;
147 :
148 : // For D insist preavg is meaningful (5 minutes or user-supplied)
149 0 : if (preavg()<0.0)
150 0 : preavg()=300.0;
151 :
152 : /*
153 : cout << boolalpha;
154 : cout << endl << endl;
155 : cout << "useGenericGatherForSolve = " << useGenericGatherForSolve() << endl;
156 : cout << "useGenericSolveOne = " << useGenericSolveOne() << endl;
157 : */
158 :
159 0 : }
160 :
161 :
162 0 : void DJones::calcOneJones(Vector<Complex>& mat, Vector<Bool>& mOk,
163 : const Vector<Complex>& par, const Vector<Bool>& pOk) {
164 :
165 0 : if (prtlev()>10) cout << " D::calcOneJones(vb)" << endl;
166 :
167 : // Only if both pars are good, form matrix
168 0 : if (allEQ(pOk,true)) {
169 :
170 : // On-diag = 1
171 0 : mat(0)=mat(3)=Complex(1.0);
172 : // Off-diag = par
173 0 : mat(1)=par(0);
174 0 : mat(2)=par(1);
175 :
176 0 : mOk=true;
177 :
178 : }
179 : else {
180 0 : mat=Complex(0.0);
181 0 : mOk=false;
182 : }
183 0 : }
184 :
185 0 : void DJones::guessPar(VisBuffer& vb) {
186 :
187 0 : if (prtlev()>4) cout << " D::guessPar(vb)" << endl;
188 :
189 : // TBD: Should we use a common wt for the X-hands??
190 :
191 : // First guess is zero D-terms
192 0 : solveCPar()=0.0;
193 0 : solveParOK()=true;
194 :
195 0 : if (jonesType()==Jones::GenLinear) {
196 0 : vb.weightMat().row(0)=0.0;
197 0 : vb.weightMat().row(3)=0.0;
198 : }
199 :
200 0 : if (solvePol()) {
201 :
202 : // solvePol() tells us how many source pol parameters
203 0 : srcPolPar().resize(solvePol());
204 :
205 : // The following assumes the MODEL_DATA has been
206 : // corrupted by P
207 :
208 0 : LSQFit fit(2,LSQComplex());
209 0 : Vector<Complex> ce(2);
210 0 : ce(0)=Complex(1.0);
211 0 : Complex d,md;
212 0 : Float wt,a(0.0);
213 0 : for (Int irow=0;irow<vb.nRow();++irow) {
214 0 : if (!vb.flagRow()(irow) &&
215 0 : vb.antenna1()(irow)!=vb.antenna2()(irow)) {
216 0 : for (Int ich=0;ich<vb.nChannel();++ich) {
217 0 : if (!vb.flag()(ich,irow)) {
218 0 : for (Int icorr=1;icorr<3;++icorr) {
219 0 : md=vb.modelVisCube()(icorr,ich,irow);
220 0 : if (icorr==2) md=conj(md);
221 0 : a=abs(md);
222 0 : if (a>0.0) {
223 0 : wt=Double(vb.weightMat()(icorr,irow));
224 0 : if (wt>0.0) {
225 0 : d=vb.visCube()(icorr,ich,irow);
226 0 : if (icorr==2) d=conj(d);
227 0 : if (abs(d)>0.0) {
228 :
229 0 : ce(1)=md;
230 0 : fit.makeNorm(ce.data(),wt,d,LSQFit::COMPLEX);
231 :
232 : } // abs(d)>0
233 : } // wt>0
234 : } // a>0
235 : } // icorr
236 : } // !flag
237 : } // ich
238 : } // !flagRow
239 : } // row
240 :
241 : uInt rank;
242 0 : Bool ok = fit.invert(rank);
243 :
244 0 : Complex sol[2];
245 0 : if (ok)
246 0 : fit.solve(sol);
247 : else
248 0 : throw(AipsError("Source polarization solution is singular; try solving for D-terms only."));
249 :
250 0 : if (solvePol()==1 && a>0.0)
251 0 : srcPolPar()(0)=Complex(arg(sol[1]));
252 0 : else if (solvePol()==2) {
253 0 : srcPolPar()(0)=Complex(real(sol[1]));
254 0 : srcPolPar()(1)=Complex(imag(sol[1]));
255 : }
256 :
257 : // Log source polarization solution
258 0 : reportSolvedQU();
259 :
260 :
261 0 : } // solvePol()?
262 :
263 0 : }
264 :
265 :
266 0 : void DJones::setUpForPolSolve(vi::VisBuffer2& vb) {
267 :
268 : // Divide model and data by (scalar) stokes I model (which may be resolved!),
269 : // and set model cross-hands to (1,0) so we can solve for fractional
270 : // pol factors.
271 :
272 : // NB: this is circ-basis-specific!!
273 :
274 : // NB: this method assumes vCC.set(vC) (corrected data already set)
275 :
276 : // Only if solving for Q an U
277 : // (leave cross-hands alone if just solving for X)
278 0 : if (solvePol()>1) {
279 :
280 0 : Int nCorr(vb.nCorrelations());
281 0 : Vector<Float> ampCorr(nCorr);
282 0 : Vector<Int> n(nCorr,0);
283 0 : Complex sI(0.0);
284 0 : for (rownr_t irow=0;irow<vb.nRows();++irow) {
285 0 : if (!vb.flagRow()(irow)) {
286 0 : ampCorr=0.0f;
287 0 : n=0;
288 0 : for (Int ich=0;ich<vb.nChannels();++ich) {
289 :
290 0 : Vector<Bool> fl(vb.flagCube()(Slice(),Slice(ich,1,1),Slice(irow,1,1)));
291 0 : Vector<Float> wt(vb.weightSpectrum()(Slice(),Slice(ich,1,1),Slice(irow,1,1)));
292 0 : Vector<Complex> vCM(vb.visCubeModel()(Slice(),Slice(ich,1,1),Slice(irow,1,1)));
293 0 : Vector<Complex> vCC(vb.visCubeCorrected()(Slice(),Slice(ich,1,1),Slice(irow,1,1)));
294 :
295 0 : if (sum(fl)>0) {
296 : // one or more correlations is flagged; cannot use
297 0 : fl.set(true);
298 0 : wt.set(0.0f);
299 0 : vCC.set(0.0f);
300 : }
301 : else {
302 : // all corrs ok, so do the calculation
303 0 : sI=(vCM(0)+vCM(3))/Complex(2.0);
304 0 : if (abs(sI)>0.0) {
305 0 : vCC/=sI;
306 0 : wt*=square(abs(sI));
307 : }
308 : else {
309 : // model stokes I is zero, cannot use
310 0 : fl.set(true);
311 0 : wt.set(0.0f);
312 0 : vCC.set(0.0f);
313 : }
314 : }
315 : // Make model unity (so parang corruption works)
316 : // NB: this is circ-basis specific!
317 : // Q: Is this correct for p-hands? (E.g., lin basis, or V!=0 in circ basis?)
318 0 : vCM.set(Complex(1.0));
319 :
320 0 : } // ich
321 : } // !flagRow
322 : } // irow
323 0 : }
324 :
325 0 : }
326 :
327 0 : void DJones::guessPar(SDBList& sdbs,const casacore::Bool& corrDepFlags) {
328 :
329 0 : if (prtlev()>4) cout << " D::guessPar(sdbs,corrDepFlags)" << endl;
330 :
331 0 : if (corrDepFlags)
332 0 : throw(AipsError("Hmmm, corrDepFlags=True in DJones::guessPar(sdbs,corrDepFlags"));
333 :
334 : // Call (old) corrDepFlags-agnostic version
335 0 : DJones::guessPar(sdbs);
336 :
337 0 : }
338 :
339 :
340 0 : void DJones::guessPar(SDBList& sdbs) {
341 :
342 : // TBD: Should we use a common wt for the X-hands??
343 :
344 : // First guess is zero D-terms
345 0 : solveCPar()=0.0;
346 0 : solveParOK()=true;
347 :
348 0 : if (jonesType()==Jones::GenLinear) {
349 : // Zero p-hand weights
350 0 : for (int isdb=0;isdb<sdbs.nSDB();++isdb) {
351 0 : SolveDataBuffer& sdb(sdbs(isdb));
352 0 : Cube<Float> wtS(sdb.weightSpectrum());
353 0 : wtS(Slice(0,2,3),Slice(),Slice())=0.0;
354 0 : }
355 : }
356 :
357 0 : if (solvePol()) {
358 :
359 0 : Int nSDB=sdbs.nSDB();
360 :
361 0 : Int nChan=sdbs.nChannels(); // insists on uniformity over SDBs
362 :
363 : // solvePol() tells us how many source pol parameters
364 0 : srcPolPar().resize(solvePol());
365 :
366 : // The following assumes the MODEL_DATA has been
367 : // corrupted by P
368 :
369 0 : LSQFit fit(2,LSQComplex());
370 0 : Vector<Complex> ce(2);
371 0 : ce(0)=Complex(1.0);
372 0 : Complex d,md;
373 0 : Float wt,a(0.0);
374 0 : for (int isdb=0;isdb<nSDB;++isdb) {
375 0 : SolveDataBuffer& sdb(sdbs(isdb));
376 0 : for (Int irow=0;irow<sdb.nRows();++irow) {
377 0 : if (!sdb.flagRow()(irow) &&
378 0 : sdb.antenna1()(irow)!=sdb.antenna2()(irow)) {
379 0 : for (Int ich=0;ich<nChan;++ich) {
380 0 : for (Int icorr=1;icorr<3;++icorr) {
381 0 : if (!sdb.flagCube()(icorr,ich,irow)) {
382 0 : md=sdb.visCubeModel()(icorr,ich,irow);
383 0 : if (icorr==2) md=conj(md);
384 0 : a=abs(md);
385 0 : if (a>0.0) {
386 0 : wt=Double(sdb.weightSpectrum()(icorr,ich,irow));
387 0 : if (wt>0.0) {
388 0 : d=sdb.visCubeCorrected()(icorr,ich,irow);
389 0 : if (icorr==2) d=conj(d);
390 0 : if (abs(d)>0.0) {
391 :
392 0 : ce(1)=md;
393 0 : fit.makeNorm(ce.data(),wt,d,LSQFit::COMPLEX);
394 :
395 : } // abs(d)>0
396 : } // wt>0
397 : } // a>0
398 : } // icorr
399 : } // !flag
400 : } // ich
401 : } // !flagRow
402 : } // row
403 : } // isdb
404 :
405 : uInt rank;
406 0 : Bool ok = fit.invert(rank);
407 :
408 0 : Complex sol[2];
409 0 : if (ok)
410 0 : fit.solve(sol);
411 : else
412 0 : throw(AipsError("Source polarization solution is singular; try solving for D-terms only."));
413 :
414 0 : Complex Pol(1.0); // for updating the model
415 0 : if (solvePol()==1 && a>0.0) {
416 0 : srcPolPar()(0)=Complex(arg(sol[1]));
417 0 : Pol=exp(Complex(0.0,real(srcPolPar()(0)))); // same as sol[1]?
418 : }
419 0 : else if (solvePol()==2) {
420 0 : srcPolPar()(0)=Complex(real(sol[1])); // Q
421 0 : srcPolPar()(1)=Complex(imag(sol[1])); // U
422 0 : Pol=Complex(real(srcPolPar()(0)),real(srcPolPar()(1))); // same as real(sol[1]),imag(sol[1])?
423 : }
424 :
425 : // Log source polarization solution
426 0 : reportSolvedQU();
427 :
428 : // update model data with this source pol estimate
429 0 : Array<Complex> RL,LR;
430 0 : Complex cPol=conj(Pol); // for LR
431 0 : for (int isdb=0;isdb<nSDB;++isdb) {
432 0 : SolveDataBuffer& sdb(sdbs(isdb));
433 0 : RL.reference(sdb.visCubeModel()(Slice(1,1,1),Slice(),Slice())); // icorr=1
434 0 : RL*=Pol;
435 0 : LR.reference(sdb.visCubeModel()(Slice(2,1,1),Slice(),Slice())); // icorr=2
436 0 : LR*=cPol;
437 : }
438 0 : } // solvePol()?
439 :
440 0 : }
441 :
442 :
443 0 : void DJones::updatePar(const Vector<Complex> dCalPar,
444 : const Vector<Complex> dSrcPar) {
445 :
446 : // Enforce no change in source parameters
447 : // before calling generic version
448 0 : Vector<Complex> dsrcpar(dSrcPar.shape());
449 0 : dsrcpar=Complex(0.0);
450 0 : SolvableVisJones::updatePar(dCalPar,dsrcpar);
451 :
452 0 : }
453 :
454 0 : void DJones::formSolveSNR() {
455 :
456 0 : solveParSNR()=0.0;
457 :
458 0 : for (Int iant=0;iant<nAnt();++iant)
459 0 : for (Int ipar=0;ipar<nPar();++ipar) {
460 0 : if (solveParOK()(ipar,0,iant) &&
461 0 : solveParErr()(ipar,0,iant)>0.0f) {
462 0 : solveParSNR()(ipar,0,iant)=1.0f/solveParErr()(ipar,0,iant);
463 : }
464 : else
465 : // Ensure F if Err<=0 (OK?)
466 0 : solveParOK()(ipar,0,iant)=false;
467 : }
468 0 : }
469 :
470 0 : void DJones::globalPostSolveTinker() {
471 :
472 : // call parent
473 0 : SolvableVisJones::globalPostSolveTinker();
474 :
475 : // if not freqdep, report results to the logger
476 0 : logResults();
477 :
478 0 : }
479 :
480 :
481 0 : void DJones::applyRefAnt() {
482 :
483 0 : SolvableVisJones::applyRefAnt();
484 0 : return;
485 :
486 : }
487 :
488 0 : void DJones::logResults() {
489 :
490 : // Don't bother, if the Ds are channelized
491 0 : if (freqDepPar()) return;
492 :
493 0 : Vector<String> rl(2);
494 0 : rl(0)="R: ";
495 0 : rl(1)="L: ";
496 :
497 0 : MSAntennaColumns antcol(ct_->antenna());
498 0 : Vector<String> antNames(antcol.name().getColumn());
499 :
500 0 : logSink() << "The instrumental polarization solutions are: " << LogIO::POST;
501 :
502 0 : logSink().output().precision(4);
503 :
504 0 : Block<String> cols(3);
505 0 : cols[0]="SPECTRAL_WINDOW_ID";
506 0 : cols[1]="TIME";
507 0 : cols[2]="ANTENNA1";
508 0 : ROCTIter ctiter(*ct_,cols);
509 :
510 0 : Int lastspw(-1);
511 0 : Double lasttime(-1.0);
512 0 : while (!ctiter.pastEnd()) {
513 0 : Int ispw=ctiter.thisSpw();
514 0 : Double time=ctiter.thisTime();
515 0 : Int a1=ctiter.thisAntenna1();
516 0 : Vector<Complex> sol;
517 0 : sol=ctiter.cparam().xyPlane(0).column(0);
518 0 : Vector<Bool> fl;
519 0 : fl=ctiter.flag().xyPlane(0).column(0);
520 :
521 0 : if (ispw!=lastspw)
522 0 : logSink() << " Spw " << ispw << ":" << endl;
523 0 : if (time !=lasttime)
524 0 : logSink() << " Time " << MVTime(time/C::day).string(MVTime::YMD,7) << ":" << endl;
525 :
526 0 : logSink().output().setf(ios::left, ios::adjustfield);
527 :
528 0 : logSink() << " Ant=" << antNames(a1) << ": ";
529 0 : for (Int ipar=0;ipar<2;++ipar) {
530 0 : logSink() << rl(ipar);
531 0 : if (!fl(ipar)) {
532 0 : logSink() << "A=";
533 0 : logSink().output().width(10);
534 0 : logSink() << abs(sol(ipar));
535 0 : logSink() << " P=";
536 0 : logSink().output().width(8);
537 0 : logSink() << arg(sol(ipar))*180.0/C::pi;
538 0 : if (ipar==0) logSink() << " ; ";
539 : }
540 : else {
541 0 : logSink().output().width(26);
542 0 : logSink() << "(flagged)" << " ";
543 : }
544 : } // ipol
545 0 : logSink() << endl;
546 0 : logSink() << LogIO::POST;
547 :
548 0 : lastspw=ispw;
549 0 : lasttime=time;
550 0 : ctiter.next();
551 0 : } // ctiter
552 0 : }
553 :
554 :
555 :
556 : // Fill the trivial DJ matrix elements
557 0 : void DJones::initTrivDJ() {
558 :
559 0 : if (prtlev()>4) cout << " D::initTrivDJ" << endl;
560 :
561 : // Must be trivial
562 0 : AlwaysAssert((trivialDJ()),AipsError);
563 :
564 : // 0 1 0 0
565 : // 0 0 & 1 0
566 :
567 0 : if (jonesType()==Jones::General) {
568 0 : diffJElem().resize(IPosition(4,4,2,1,1));
569 0 : diffJElem()=0.0;
570 0 : diffJElem()(IPosition(4,1,0,0,0))=Complex(1.0);
571 0 : diffJElem()(IPosition(4,2,1,0,0))=Complex(1.0);
572 : }
573 : else {
574 0 : diffJElem().resize(IPosition(4,2,2,1,1));
575 0 : diffJElem()=0.0;
576 0 : diffJElem()(IPosition(4,0,0,0,0))=Complex(1.0);
577 0 : diffJElem()(IPosition(4,1,1,0,0))=Complex(1.0);
578 : }
579 :
580 0 : }
581 :
582 :
583 :
584 0 : void DJones::createCorruptor(const VisIter& vi, const Record& simpar, const Int nSim)
585 : {
586 :
587 0 : LogIO os(LogOrigin("D", "createCorruptor()", WHERE));
588 0 : if (prtlev()>2) cout << " D::createCorruptor()" << endl;
589 :
590 : // this may not be the best place for this:
591 0 : solvePol_=2;
592 :
593 : // no nSim since not time dependent (yet)
594 0 : dcorruptor_p = new DJonesCorruptor();
595 0 : corruptor_p = dcorruptor_p;
596 :
597 : // call generic parent to set corr,spw,etc info
598 0 : SolvableVisCal::createCorruptor(vi,simpar,nSim);
599 :
600 0 : Int Seed(1234);
601 0 : if (simpar.isDefined("seed")) {
602 0 : Seed=simpar.asInt("seed");
603 : }
604 :
605 0 : Complex Scale(0.1,0.1); // scale of fluctuations
606 0 : if (simpar.isDefined("camp")) {
607 0 : Scale=simpar.asComplex("camp");
608 : }
609 :
610 0 : Complex Offset(0.,0.);
611 0 : if (simpar.isDefined("offset")) {
612 0 : Offset=simpar.asComplex("offset");
613 : }
614 :
615 0 : dcorruptor_p->initialize(Seed,Scale,Offset);
616 :
617 0 : }
618 :
619 :
620 :
621 : // **********************************************************
622 : // DfJones Implementations
623 : //
624 :
625 0 : DfJones::DfJones(VisSet& vs) :
626 : VisCal(vs), // virtual base
627 : VisMueller(vs), // virtual base
628 0 : DJones(vs) // immediate parent
629 : {
630 0 : if (prtlev()>2) cout << "Df::Df(vs)" << endl;
631 0 : }
632 :
633 0 : DfJones::DfJones(String msname,Int MSnAnt,Int MSnSpw) :
634 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
635 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
636 0 : DJones(msname,MSnAnt,MSnSpw) // immediate parent
637 : {
638 0 : if (prtlev()>2) cout << "Df::Df(msname,MSnAnt,MSnSpw)" << endl;
639 0 : }
640 :
641 0 : DfJones::DfJones(const MSMetaInfoForCal& msmc) :
642 : VisCal(msmc), // virtual base
643 : VisMueller(msmc), // virtual base
644 0 : DJones(msmc) // immediate parent
645 : {
646 0 : if (prtlev()>2) cout << "Df::Df(msmc)" << endl;
647 0 : }
648 :
649 0 : DfJones::DfJones(const Int& nAnt) :
650 : VisCal(nAnt),
651 : VisMueller(nAnt),
652 0 : DJones(nAnt)
653 : {
654 0 : if (prtlev()>2) cout << "Df::Df(nAnt)" << endl;
655 0 : }
656 :
657 0 : DfJones::~DfJones() {
658 0 : if (prtlev()>2) cout << "Df::~Df()" << endl;
659 0 : }
660 :
661 :
662 :
663 : // **********************************************************
664 : // DlinJones Implementations
665 : //
666 :
667 : // Constructor
668 0 : DlinJones::DlinJones(VisSet& vs) :
669 : VisCal(vs), // virtual base
670 : VisMueller(vs), // virtual base
671 0 : DJones(vs) // immediate parent
672 : {
673 0 : if (prtlev()>2) cout << "Dlin::Dlin(vs)" << endl;
674 0 : }
675 :
676 0 : DlinJones::DlinJones(String msname,Int MSnAnt,Int MSnSpw) :
677 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
678 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
679 0 : DJones(msname,MSnAnt,MSnSpw) // immediate parent
680 : {
681 0 : if (prtlev()>2) cout << "Dlin::Dlin(msname,MSnAnt,MSnSpw)" << endl;
682 0 : }
683 :
684 0 : DlinJones::DlinJones(const MSMetaInfoForCal& msmc) :
685 : VisCal(msmc), // virtual base
686 : VisMueller(msmc), // virtual base
687 0 : DJones(msmc) // immediate parent
688 : {
689 0 : if (prtlev()>2) cout << "Dlin::Dlin(msmc)" << endl;
690 0 : }
691 :
692 0 : DlinJones::DlinJones(const Int& nAnt) :
693 : VisCal(nAnt),
694 : VisMueller(nAnt),
695 0 : DJones(nAnt)
696 : {
697 0 : if (prtlev()>2) cout << "Dlin::Dlin(nAnt)" << endl;
698 0 : }
699 :
700 0 : DlinJones::~DlinJones() {
701 0 : if (prtlev()>2) cout << "Dlin::~Dlin()" << endl;
702 0 : }
703 :
704 :
705 : // **********************************************************
706 : // DflinJones Implementations
707 : //
708 :
709 : // Constructor
710 0 : DflinJones::DflinJones(VisSet& vs) :
711 : VisCal(vs), // virtual base
712 : VisMueller(vs), // virtual base
713 0 : DlinJones(vs) // immediate parent
714 : {
715 0 : if (prtlev()>2) cout << "Dflin::Dflin(vs)" << endl;
716 0 : }
717 :
718 0 : DflinJones::DflinJones(String msname,Int MSnAnt,Int MSnSpw) :
719 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
720 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
721 0 : DlinJones(msname,MSnAnt,MSnSpw) // immediate parent
722 : {
723 0 : if (prtlev()>2) cout << "Dflin::Dflin(msname,MSnAnt,MSnSpw)" << endl;
724 0 : }
725 :
726 0 : DflinJones::DflinJones(const MSMetaInfoForCal& msmc) :
727 : VisCal(msmc), // virtual base
728 : VisMueller(msmc), // virtual base
729 0 : DlinJones(msmc) // immediate parent
730 : {
731 0 : if (prtlev()>2) cout << "Dflin::Dflin(msmc)" << endl;
732 0 : }
733 :
734 0 : DflinJones::DflinJones(const Int& nAnt) :
735 : VisCal(nAnt),
736 : VisMueller(nAnt),
737 0 : DlinJones(nAnt)
738 : {
739 0 : if (prtlev()>2) cout << "Dflin::Dflin(nAnt)" << endl;
740 0 : }
741 :
742 0 : DflinJones::~DflinJones() {
743 0 : if (prtlev()>2) cout << "Dflin::~Dflin()" << endl;
744 0 : }
745 :
746 :
747 :
748 :
749 : // **********************************************************
750 : // DllsJones Implementations
751 : //
752 :
753 : // Constructor
754 0 : DllsJones::DllsJones(VisSet& vs) :
755 : VisCal(vs), // virtual base
756 : VisMueller(vs), // virtual base
757 0 : DJones(vs) // immediate parent
758 : {
759 0 : if (prtlev()>2) cout << "Dlls::Dlls(vs)" << endl;
760 0 : }
761 :
762 0 : DllsJones::DllsJones(String msname,Int MSnAnt,Int MSnSpw) :
763 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
764 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
765 0 : DJones(msname,MSnAnt,MSnSpw) // immediate parent
766 : {
767 0 : if (prtlev()>2) cout << "Dlls::Dlls(msname,MSnAnt,MSnSpw)" << endl;
768 0 : }
769 :
770 0 : DllsJones::DllsJones(const MSMetaInfoForCal& msmc) :
771 : VisCal(msmc), // virtual base
772 : VisMueller(msmc), // virtual base
773 0 : DJones(msmc) // immediate parent
774 : {
775 0 : if (prtlev()>2) cout << "Dlls::Dlls(msmc)" << endl;
776 0 : }
777 :
778 0 : DllsJones::DllsJones(const Int& nAnt) :
779 : VisCal(nAnt),
780 : VisMueller(nAnt),
781 0 : DJones(nAnt)
782 : {
783 0 : if (prtlev()>2) cout << "Dlls::Dlls(nAnt)" << endl;
784 0 : }
785 :
786 0 : DllsJones::~DllsJones() {
787 0 : if (prtlev()>2) cout << "Dlls::~Dlls()" << endl;
788 0 : }
789 :
790 0 : void DllsJones::solveOneVB(const VisBuffer& vb) {
791 :
792 0 : Int nChan=vb.nChannel();
793 :
794 0 : solveAllCPar()=Complex(0.0);
795 0 : solveAllParOK()=false;
796 :
797 : // Solve per chan
798 0 : Vector<Complex> cEq(2);
799 0 : Vector<uInt> cEqIdx(2);
800 0 : Vector<Complex> obs(2);
801 0 : Vector<Float> wt(2);
802 0 : for (Int ich=0;ich<nChan;++ich) {
803 : //Int currCh=vb.channel()(ich);
804 : //cout << "Ch=" << currCh;
805 :
806 0 : solveCPar().reference(solveAllCPar()(Slice(),Slice(ich,1,1),Slice()));
807 0 : solveParOK().reference(solveAllParOK()(Slice(),Slice(ich,1,1),Slice()));
808 0 : solveParOK().set(false);
809 :
810 0 : LSQFit lsq(2*nAnt(),LSQComplex(),0);
811 :
812 0 : Int ng(0);
813 0 : for (Int irow=0;irow<vb.nRow();++irow) {
814 0 : if (!vb.flagRow()(irow)) {
815 :
816 : // Discern this row's antennas, avoid ACs
817 0 : Int a1=vb.antenna1()(irow);
818 0 : Int a2=vb.antenna2()(irow);
819 0 : if (a1!=a2) {
820 :
821 : // If not flagged...
822 0 : if (!vb.flag()(ich,irow)) {
823 0 : ++ng;
824 :
825 : // simple-minded OK-setting, for now
826 0 : solveParOK().xyPlane(a1).set(true);
827 0 : solveParOK().xyPlane(a2).set(true);
828 :
829 : // Discern D indices for this baseline
830 0 : Int xi=2*a1;
831 0 : Int yi=2*a1+1;
832 0 : Int xj=2*a2;
833 0 : Int yj=2*a2+1;
834 : /*
835 : cout << "currCh="<<currCh
836 : << " irow="<<irow
837 : << " a1="<<a1 << " a2="<<a2
838 : << " XY: xi/yj*="<< xi<<"/"<<yj
839 : << " YX: yi/xj*="<< yi<<"/"<<xj
840 : << endl;
841 : */
842 :
843 0 : wt(0)=vb.weightMat()(1,irow);
844 0 : wt(1)=vb.weightMat()(2,irow);
845 :
846 0 : cEq(0)=vb.modelVisCube()(0,ich,irow);
847 0 : cEq(1)=vb.modelVisCube()(3,ich,irow);
848 0 : obs(0)=vb.visCube()(1,ich,irow)-vb.modelVisCube()(1,ich,irow);
849 0 : obs(1)=vb.visCube()(2,ich,irow)-vb.modelVisCube()(2,ich,irow);
850 :
851 0 : cEqIdx(0)=yj; cEqIdx(1)=xi;
852 0 : lsq.makeNorm(2,cEqIdx.data(),cEq.data(),wt(0),obs(0),LSQFit::Complex());
853 0 : cEqIdx(0)=yi; cEqIdx(1)=xj;
854 0 : lsq.makeNorm(2,cEqIdx.data(),cEq.data(),wt(1),conj(obs(1)),LSQFit::Complex());
855 :
856 : }
857 :
858 : } // !AC
859 : } // !flagRow
860 : } // irow
861 :
862 0 : if (ng>0) {
863 :
864 : // Turn the crank on the solver
865 : uInt rank;
866 0 : lsq.invert(rank,false); // true); // SVD?
867 0 : Cube<Complex> Dest(2,1,nAnt());
868 :
869 0 : lsq.solve(Dest.data());
870 :
871 : //if (Int(rank)<nAnt())
872 : //cout << " Need to reference!!";
873 :
874 : // Conjugate Dy's because we solved for their conjugates
875 0 : Dest(Slice(1,1,1),Slice(),Slice())=conj(Dest(Slice(1,1,1),Slice(),Slice()));
876 :
877 : // Store the result
878 0 : solveCPar()=Dest;
879 : //solveParOK().set(true); // we set this more 'carefully' above
880 0 : }
881 : //cout << endl;
882 :
883 0 : } // ich
884 :
885 0 : }
886 :
887 0 : void DllsJones::solveOneSDB(SolveDataBuffer& sdb) {
888 :
889 0 : Int nChan=sdb.nChannels();
890 :
891 0 : solveAllCPar()=Complex(0.0);
892 0 : solveAllParOK()=false;
893 :
894 : // Solve per chan
895 0 : Vector<Complex> cEq(2);
896 0 : Vector<uInt> cEqIdx(2);
897 0 : Vector<Complex> obs(2);
898 0 : Vector<Float> wt(2);
899 0 : for (Int ich=0;ich<nChan;++ich) {
900 :
901 0 : solveCPar().reference(solveAllCPar()(Slice(),Slice(ich,1,1),Slice()));
902 0 : solveParOK().reference(solveAllParOK()(Slice(),Slice(ich,1,1),Slice()));
903 0 : solveParOK().set(false);
904 :
905 0 : LSQFit lsq(2*nAnt(),LSQComplex(),0);
906 :
907 0 : Int ng(0);
908 0 : for (Int irow=0;irow<sdb.nRows();++irow) {
909 0 : if (!sdb.flagRow()(irow)) {
910 :
911 : // Discern this row's antennas, avoid ACs
912 0 : Int a1=sdb.antenna1()(irow);
913 0 : Int a2=sdb.antenna2()(irow);
914 0 : if (a1!=a2) {
915 :
916 : // If not flagged...
917 0 : if (!sdb.flagCube()(0,ich,irow) &&
918 0 : !sdb.flagCube()(1,ich,irow) &&
919 0 : !sdb.flagCube()(2,ich,irow) &&
920 0 : !sdb.flagCube()(3,ich,irow) ) {
921 0 : ++ng;
922 :
923 : // simple-minded OK-setting, for now
924 0 : solveParOK().xyPlane(a1).set(true);
925 0 : solveParOK().xyPlane(a2).set(true);
926 :
927 : // Discern D indices for this baseline
928 0 : Int xi=2*a1;
929 0 : Int yi=2*a1+1;
930 0 : Int xj=2*a2;
931 0 : Int yj=2*a2+1;
932 : /*
933 : cout << "currCh="<<currCh
934 : << " irow="<<irow
935 : << " a1="<<a1 << " a2="<<a2
936 : << " XY: xi/yj*="<< xi<<"/"<<yj
937 : << " YX: yi/xj*="<< yi<<"/"<<xj
938 : << endl;
939 : */
940 :
941 0 : wt(0)=sdb.weightSpectrum()(1,ich,irow);
942 0 : wt(1)=sdb.weightSpectrum()(2,ich,irow);
943 :
944 0 : cEq(0)=sdb.visCubeModel()(0,ich,irow);
945 0 : cEq(1)=sdb.visCubeModel()(3,ich,irow);
946 0 : obs(0)=sdb.visCubeCorrected()(1,ich,irow)-sdb.visCubeModel()(1,ich,irow);
947 0 : obs(1)=sdb.visCubeCorrected()(2,ich,irow)-sdb.visCubeModel()(2,ich,irow);
948 :
949 0 : cEqIdx(0)=yj; cEqIdx(1)=xi;
950 0 : lsq.makeNorm(2,cEqIdx.data(),cEq.data(),wt(0),obs(0),LSQFit::Complex());
951 0 : cEqIdx(0)=yi; cEqIdx(1)=xj;
952 0 : lsq.makeNorm(2,cEqIdx.data(),cEq.data(),wt(1),conj(obs(1)),LSQFit::Complex());
953 :
954 : }
955 :
956 : } // !AC
957 : } // !flagRow
958 : } // irow
959 :
960 0 : if (ng>0) {
961 :
962 : // Turn the crank on the solver
963 : uInt rank;
964 0 : lsq.invert(rank,false); // true); // SVD?
965 0 : Cube<Complex> Dest(2,1,nAnt());
966 :
967 0 : lsq.solve(Dest.data());
968 :
969 : //if (Int(rank)<nAnt())
970 : //cout << " Need to reference!!";
971 :
972 : // Conjugate Dy's because we solved for their conjugates
973 0 : Dest(Slice(1,1,1),Slice(),Slice())=conj(Dest(Slice(1,1,1),Slice(),Slice()));
974 :
975 : // Store the result
976 0 : solveCPar()=Dest;
977 : //solveParOK().set(true); // we set this more 'carefully' above
978 0 : }
979 : //cout << endl;
980 :
981 0 : } // ich
982 :
983 0 : }
984 :
985 0 : void DllsJones::solveOne(SDBList& sdbs) {
986 :
987 : // In general, there will be many SDBs in the SDBList, usually
988 : // over time (parang). (If combine='spw', then there will be
989 : // multiple SDBs over spw; in this case, we have to be clearer
990 : // about what we mean by nChan below....)
991 :
992 : // The number of SDBs in the SDBList
993 0 : Int nSDB=sdbs.nSDB();
994 :
995 : //cout << "nSDB=" << nSDB << endl;
996 :
997 : // This will insist on uniformity over SDBs (usually the same spw)
998 0 : Int nChan=sdbs.nChannels();
999 :
1000 0 : solveAllCPar()=Complex(0.0);
1001 0 : solveAllParOK()=false;
1002 :
1003 : // Solve per chan
1004 0 : Vector<Complex> cEq(2);
1005 0 : Vector<uInt> cEqIdx(2);
1006 0 : Vector<Complex> obs(2);
1007 0 : Vector<Float> wt(2);
1008 0 : for (Int ich=0;ich<nChan;++ich) {
1009 :
1010 0 : solveCPar().reference(solveAllCPar()(Slice(),Slice(ich,1,1),Slice()));
1011 0 : solveParOK().reference(solveAllParOK()(Slice(),Slice(ich,1,1),Slice()));
1012 0 : solveParOK().set(false);
1013 :
1014 0 : LSQFit lsq(2*nAnt(),LSQComplex(),0);
1015 :
1016 0 : Int ng(0);
1017 0 : Float wtsum(0.0);
1018 0 : for (Int isdb=0;isdb<nSDB;++isdb) {
1019 0 : SolveDataBuffer& sdb(sdbs(isdb));
1020 0 : for (Int irow=0;irow<sdb.nRows();++irow) {
1021 0 : if (!sdb.flagRow()(irow)) {
1022 :
1023 : // Discern this row's antennas, avoid ACs
1024 0 : Int a1=sdb.antenna1()(irow);
1025 0 : Int a2=sdb.antenna2()(irow);
1026 0 : if (a1!=a2) {
1027 :
1028 : // If not flagged...
1029 0 : if (!sdb.flagCube()(0,ich,irow) &&
1030 0 : !sdb.flagCube()(1,ich,irow) &&
1031 0 : !sdb.flagCube()(2,ich,irow) &&
1032 0 : !sdb.flagCube()(3,ich,irow) ) {
1033 0 : ++ng;
1034 :
1035 : // simple-minded OK-setting, for now
1036 0 : solveParOK().xyPlane(a1).set(true);
1037 0 : solveParOK().xyPlane(a2).set(true);
1038 :
1039 : // Discern D indices for this baseline
1040 0 : Int xi=2*a1;
1041 0 : Int yi=2*a1+1;
1042 0 : Int xj=2*a2;
1043 0 : Int yj=2*a2+1;
1044 : /*
1045 : cout << "currCh="<<currCh
1046 : << " irow="<<irow
1047 : << " a1="<<a1 << " a2="<<a2
1048 : << " XY: xi/yj*="<< xi<<"/"<<yj
1049 : << " YX: yi/xj*="<< yi<<"/"<<xj
1050 : << endl;
1051 : */
1052 :
1053 0 : wt(0)=sdb.weightSpectrum()(1,ich,irow);
1054 0 : wt(1)=sdb.weightSpectrum()(2,ich,irow);
1055 :
1056 0 : cEq(0)=sdb.visCubeModel()(0,ich,irow);
1057 0 : cEq(1)=sdb.visCubeModel()(3,ich,irow);
1058 0 : obs(0)=sdb.visCubeCorrected()(1,ich,irow)-sdb.visCubeModel()(1,ich,irow);
1059 0 : obs(1)=sdb.visCubeCorrected()(2,ich,irow)-sdb.visCubeModel()(2,ich,irow);
1060 :
1061 0 : cEqIdx(0)=yj; cEqIdx(1)=xi;
1062 0 : lsq.makeNorm(2,cEqIdx.data(),cEq.data(),wt(0),obs(0),LSQFit::Complex());
1063 0 : cEqIdx(0)=yi; cEqIdx(1)=xj;
1064 0 : lsq.makeNorm(2,cEqIdx.data(),cEq.data(),wt(1),conj(obs(1)),LSQFit::Complex());
1065 :
1066 0 : wtsum+=(wt(0)+wt(1));
1067 :
1068 : }
1069 :
1070 : } // !AC
1071 : } // !flagRow
1072 : } // irow
1073 : } // isdb
1074 :
1075 : //cout << "wtsum=" << wtsum << endl;
1076 :
1077 :
1078 0 : if (ng>0) {
1079 :
1080 : // Turn the crank on the solver
1081 : uInt rank;
1082 0 : lsq.invert(rank,false); // true); // SVD?
1083 0 : Cube<Complex> Dest(2,1,nAnt());
1084 :
1085 0 : lsq.solve(Dest.data());
1086 :
1087 : //if (Int(rank)<nAnt())
1088 : //cout << " Need to reference!!";
1089 :
1090 : // Conjugate Dy's because we solved for their conjugates
1091 0 : Dest(Slice(1,1,1),Slice(),Slice())=conj(Dest(Slice(1,1,1),Slice(),Slice()));
1092 :
1093 : // Store the result
1094 0 : solveCPar()=Dest;
1095 : //solveParOK().set(true); // we set this more 'carefully' above
1096 0 : }
1097 : //cout << endl;
1098 :
1099 0 : } // ich
1100 :
1101 0 : }
1102 :
1103 :
1104 : // **********************************************************
1105 : // DfllsJones Implementations
1106 : //
1107 :
1108 : // Constructor
1109 0 : DfllsJones::DfllsJones(VisSet& vs) :
1110 : VisCal(vs), // virtual base
1111 : VisMueller(vs), // virtual base
1112 0 : DllsJones(vs) // immediate parent
1113 : {
1114 0 : if (prtlev()>2) cout << "Dflls::Dflls(vs)" << endl;
1115 0 : }
1116 :
1117 0 : DfllsJones::DfllsJones(String msname,Int MSnAnt,Int MSnSpw) :
1118 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
1119 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
1120 0 : DllsJones(msname,MSnAnt,MSnSpw) // immediate parent
1121 : {
1122 0 : if (prtlev()>2) cout << "Dflls::Dflls(msname,MSnAnt,MSnSpw)" << endl;
1123 0 : }
1124 :
1125 0 : DfllsJones::DfllsJones(const MSMetaInfoForCal& msmc) :
1126 : VisCal(msmc), // virtual base
1127 : VisMueller(msmc), // virtual base
1128 0 : DllsJones(msmc) // immediate parent
1129 : {
1130 0 : if (prtlev()>2) cout << "Dflls::Dflls(msmc)" << endl;
1131 0 : }
1132 :
1133 0 : DfllsJones::DfllsJones(const Int& nAnt) :
1134 : VisCal(nAnt),
1135 : VisMueller(nAnt),
1136 0 : DllsJones(nAnt)
1137 : {
1138 0 : if (prtlev()>2) cout << "Dflls::Dflls(nAnt)" << endl;
1139 0 : }
1140 :
1141 0 : DfllsJones::~DfllsJones() {
1142 0 : if (prtlev()>2) cout << "Dflls::~Dflls()" << endl;
1143 0 : }
1144 :
1145 : /*
1146 :
1147 : // **********************************************************
1148 : // XMueller: positiona angle for circulars
1149 : //
1150 :
1151 : XMueller::XMueller(VisSet& vs) :
1152 : VisCal(vs), // virtual base
1153 : VisMueller(vs), // virtual base
1154 : SolvableVisMueller(vs) // immediate parent
1155 : {
1156 : if (prtlev()>2) cout << "X::X(vs)" << endl;
1157 : }
1158 :
1159 : XMueller::XMueller(String msname,Int MSnAnt,Int MSnSpw) :
1160 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
1161 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
1162 : SolvableVisMueller(msname,MSnAnt,MSnSpw) // immediate parent
1163 : {
1164 : if (prtlev()>2) cout << "X::X(msname,MSnAnt,MSnSpw)" << endl;
1165 : }
1166 :
1167 : XMueller::XMueller(const MSMetaInfoForCal& msmc) :
1168 : VisCal(msmc), // virtual base
1169 : VisMueller(msmc), // virtual base
1170 : SolvableVisMueller(msmc) // immediate parent
1171 : {
1172 : if (prtlev()>2) cout << "X::X(msmc)" << endl;
1173 : }
1174 :
1175 : XMueller::XMueller(const Int& nAnt) :
1176 : VisCal(nAnt),
1177 : VisMueller(nAnt),
1178 : SolvableVisMueller(nAnt)
1179 : {
1180 : if (prtlev()>2) cout << "X::X(nAnt)" << endl;
1181 : }
1182 :
1183 : XMueller::~XMueller() {
1184 : if (prtlev()>2) cout << "X::~X()" << endl;
1185 : }
1186 :
1187 : void XMueller::setApply(const Record& apply) {
1188 :
1189 : SolvableVisCal::setApply(apply);
1190 :
1191 : // Force calwt to false
1192 : calWt()=false;
1193 :
1194 : }
1195 :
1196 :
1197 : void XMueller::setSolve(const Record& solvepar) {
1198 :
1199 : cout << "XMueller: parType() = " << this->parType() << endl;
1200 :
1201 : SolvableVisCal::setSolve(solvepar);
1202 :
1203 : // Force calwt to false
1204 : calWt()=false;
1205 :
1206 : // For X insist preavg is meaningful (5 minutes or user-supplied)
1207 : if (preavg()<0.0)
1208 : preavg()=300.0;
1209 :
1210 :
1211 : cout << "ct_ = " << ct_ << endl;
1212 :
1213 :
1214 : }
1215 :
1216 : void XMueller::newselfSolve(VisSet& vs, VisEquation& ve) {
1217 :
1218 : if (prtlev()>4) cout << " M::selfSolve(ve)" << endl;
1219 :
1220 : // Inform logger/history
1221 : logSink() << "Solving for " << typeName()
1222 : << LogIO::POST;
1223 :
1224 : // Initialize the svc according to current VisSet
1225 : // (this counts intervals, sizes CalSet)
1226 : Vector<Int> nChunkPerSol;
1227 : Int nSol = sizeUpSolve(vs,nChunkPerSol);
1228 :
1229 : // Create the Caltable
1230 : createMemCalTable();
1231 :
1232 : // The iterator, VisBuffer
1233 : VisIter& vi(vs.iter());
1234 : VisBuffer vb(vi);
1235 :
1236 : // cout << "nSol = " << nSol << endl;
1237 : // cout << "nChunkPerSol = " << nChunkPerSol << endl;
1238 :
1239 : Vector<Int> slotidx(vs.numberSpw(),-1);
1240 :
1241 : Int nGood(0);
1242 : vi.originChunks();
1243 : for (Int isol=0;isol<nSol && vi.moreChunks();++isol) {
1244 :
1245 : // Arrange to accumulate
1246 : VisBuffAccumulator vba(nAnt(),preavg(),false);
1247 :
1248 : for (Int ichunk=0;ichunk<nChunkPerSol(isol);++ichunk) {
1249 :
1250 : // Current _chunk_'s spw
1251 : Int spw(vi.spectralWindow());
1252 :
1253 : // Abort if we encounter a spw for which a priori cal not available
1254 : if (!ve.spwOK(spw))
1255 : throw(AipsError("Pre-applied calibration not available for at least 1 spw. Check spw selection carefully."));
1256 :
1257 :
1258 : // Collapse each timestamp in this chunk according to VisEq
1259 : // with calibration and averaging
1260 : for (vi.origin(); vi.more(); vi++) {
1261 :
1262 : // Force read of the field Id
1263 : vb.fieldId();
1264 :
1265 : // This forces the data/model/wt I/O, and applies
1266 : // any prior calibrations
1267 : ve.collapse(vb);
1268 :
1269 : // If permitted/required by solvable component, normalize
1270 : //if (normalizable())
1271 : // vb.normalize();
1272 :
1273 : // If this solve not freqdep, and channels not averaged yet, do so
1274 : if (!freqDepMat() && vb.nChannel()>1)
1275 : vb.freqAveCubes();
1276 :
1277 : // Accumulate collapsed vb in a time average
1278 : vba.accumulate(vb);
1279 : }
1280 : // Advance the VisIter, if possible
1281 : if (vi.moreChunks()) vi.nextChunk();
1282 :
1283 : }
1284 :
1285 : // Finalize the averged VisBuffer
1286 : vba.finalizeAverage();
1287 :
1288 : // The VisBuffer to solve with
1289 : VisBuffer& svb(vba.aveVisBuff());
1290 :
1291 : // Establish meta-data for this interval
1292 : // (some of this may be used _during_ solve)
1293 : // (this sets currSpw() in the SVC)
1294 : Bool vbOk=syncSolveMeta(svb,-1);
1295 :
1296 : Int thisSpw=spwMap()(svb.spectralWindow());
1297 : slotidx(thisSpw)++;
1298 :
1299 : // We are actually solving for all channels simultaneously
1300 : solveCPar().reference(solveAllCPar());
1301 : solveParOK().reference(solveAllParOK());
1302 : solveParErr().reference(solveAllParErr());
1303 : solveParSNR().reference(solveAllParSNR());
1304 :
1305 : // Fill solveCPar() with 1, nominally, and flagged
1306 : solveCPar()=Complex(1.0);
1307 : solveParOK()=false;
1308 :
1309 : if (vbOk && svb.nRow()>0) {
1310 :
1311 : // solve for the R-L phase term in the current VB
1312 : solveOneVB(svb);
1313 :
1314 : if (solveParOK()(0,0,0))
1315 : logSink() << "Position angle offset solution for "
1316 : << msmc().fieldName(currField())
1317 : << " (spw = " << currSpw() << ") = "
1318 : << arg(solveCPar()(0,0,0))*180.0/C::pi/2.0
1319 : << " deg."
1320 : << LogIO::POST;
1321 : else
1322 : logSink() << "Position angle offset solution for "
1323 : << msmc().fieldName(currField())
1324 : << " (spw = " << currSpw() << ") "
1325 : << " was not determined (insufficient data)."
1326 : << LogIO::POST;
1327 :
1328 : nGood++;
1329 : }
1330 :
1331 : keepNCT();
1332 :
1333 : }
1334 :
1335 : logSink() << " Found good "
1336 : << typeName() << " solutions in "
1337 : << nGood << " intervals."
1338 : << LogIO::POST;
1339 :
1340 : // Store whole of result in a caltable
1341 : if (nGood==0)
1342 : logSink() << "No output calibration table written."
1343 : << LogIO::POST;
1344 : else {
1345 :
1346 : // Do global post-solve tinkering (e.g., phase-only, normalization, etc.)
1347 : // TBD
1348 : // globalPostSolveTinker();
1349 :
1350 : // write the table
1351 : storeNCT();
1352 :
1353 : }
1354 :
1355 : }
1356 :
1357 : void XMueller::calcAllMueller() {
1358 :
1359 : // cout << "currMElem().shape() = " << currMElem().shape() << endl;
1360 :
1361 : // Put the phase factor into the cross-hand diagonals
1362 : // (1,0) for the para-hands
1363 : IPosition blc(3,0,0,0);
1364 : IPosition trc(3,0,nChanMat()-1,nElem()-1);
1365 : currMElem()(blc,trc)=Complex(1.0);
1366 :
1367 : blc(0)=trc(0)=1;
1368 : currMElem()(blc,trc)=currCPar()(0,0,0);
1369 : blc(0)=trc(0)=2;
1370 : currMElem()(blc,trc)=conj(currCPar()(0,0,0));
1371 :
1372 : blc(0)=trc(0)=3;
1373 : currMElem()(blc,trc)=Complex(1.0);
1374 :
1375 : currMElemOK()=true;
1376 :
1377 : }
1378 :
1379 :
1380 : void XMueller::solveOneVB(const VisBuffer& vb) {
1381 :
1382 : // This just a simple average of the cross-hand
1383 : // visbilities...
1384 :
1385 : Complex d,md;
1386 : Float wt,a;
1387 : DComplex rl(0.0),lr(0.0);
1388 : Double sumwt(0.0);
1389 : for (Int irow=0;irow<vb.nRow();++irow) {
1390 : if (!vb.flagRow()(irow) &&
1391 : vb.antenna1()(irow)!=vb.antenna2()(irow)) {
1392 :
1393 : for (Int ich=0;ich<vb.nChannel();++ich) {
1394 : if (!vb.flag()(ich,irow)) {
1395 :
1396 : // A common weight for both crosshands
1397 : // TBD: we should probably consider this carefully...
1398 : // (also in D::guessPar...)
1399 : wt=Double(vb.weightMat()(1,irow)+
1400 : vb.weightMat()(2,irow))/2.0;
1401 :
1402 : // correct weight for model normalization
1403 : a=abs(vb.modelVisCube()(1,ich,irow));
1404 : wt*=(a*a);
1405 :
1406 : if (wt>0.0) {
1407 : // Cross-hands only
1408 : for (Int icorr=1;icorr<3;++icorr) {
1409 : md=vb.modelVisCube()(icorr,ich,irow);
1410 : d=vb.visCube()(icorr,ich,irow);
1411 :
1412 : if (abs(d)>0.0) {
1413 :
1414 : if (icorr==1)
1415 : rl+=DComplex(Complex(wt)*d/md);
1416 : else
1417 : lr+=DComplex(Complex(wt)*d/md);
1418 :
1419 : sumwt+=Double(wt);
1420 :
1421 : } // abs(d)>0
1422 : } // icorr
1423 : } // wt>0
1424 : } // !flag
1425 : } // ich
1426 : } // !flagRow
1427 : } // row
1428 :
1429 :
1430 : // cout << "spw = " << currSpw() << endl;
1431 : // cout << " rl = " << rl << " " << arg(rl)*180.0/C::pi << endl;
1432 : // cout << " lr = " << lr << " " << arg(lr)*180.0/C::pi << endl;
1433 :
1434 :
1435 : // combine lr with rl
1436 : rl+=conj(lr);
1437 :
1438 : // Normalize to unit amplitude
1439 : // (note that the phase result is insensitive to sumwt)
1440 : Double amp=abs(rl);
1441 : if (sumwt>0 && amp>0.0) {
1442 : rl/=DComplex(amp);
1443 :
1444 : solveCPar()=Complex(rl);
1445 : solveParOK()=true;
1446 : }
1447 :
1448 : }
1449 :
1450 :
1451 :
1452 : // **********************************************************
1453 : // XJones: position angle for circulars (antenna-based
1454 : //
1455 :
1456 : XJones::XJones(VisSet& vs) :
1457 : VisCal(vs), // virtual base
1458 : VisMueller(vs), // virtual base
1459 : SolvableVisJones(vs) // immediate parent
1460 : {
1461 : if (prtlev()>2) cout << "X::X(vs)" << endl;
1462 : }
1463 :
1464 : XJones::XJones(String msname,Int MSnAnt,Int MSnSpw) :
1465 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
1466 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
1467 : SolvableVisJones(msname,MSnAnt,MSnSpw) // immediate parent
1468 : {
1469 : if (prtlev()>2) cout << "X::X(msname,MSnAnt,MSnSpw)" << endl;
1470 : }
1471 :
1472 : XJones::XJones(const MSMetaInfoForCal& msmc) :
1473 : VisCal(msmc), // virtual base
1474 : VisMueller(msmc), // virtual base
1475 : SolvableVisJones(msmc) // immediate parent
1476 : {
1477 : if (prtlev()>2) cout << "X::X(msmc)" << endl;
1478 : }
1479 :
1480 : XJones::XJones(const Int& nAnt) :
1481 : VisCal(nAnt),
1482 : VisMueller(nAnt),
1483 : SolvableVisJones(nAnt)
1484 : {
1485 : if (prtlev()>2) cout << "X::X(nAnt)" << endl;
1486 : }
1487 :
1488 : XJones::~XJones() {
1489 : if (prtlev()>2) cout << "X::~X()" << endl;
1490 : }
1491 :
1492 : void XJones::setApply(const Record& apply) {
1493 :
1494 : SolvableVisCal::setApply(apply);
1495 :
1496 : // Force calwt to false
1497 : calWt()=false;
1498 :
1499 : }
1500 :
1501 :
1502 : void XJones::setSolve(const Record& solvepar) {
1503 :
1504 : SolvableVisCal::setSolve(solvepar);
1505 :
1506 : // Force calwt to false
1507 : calWt()=false;
1508 :
1509 : // For X insist preavg is meaningful (5 minutes or user-supplied)
1510 : if (preavg()<0.0)
1511 : preavg()=300.0;
1512 :
1513 : // Force refant to none (==-1), because it is meaningless to
1514 : // apply a refant to an X solution
1515 : if (refant()>-1) {
1516 : logSink() << ". (Ignoring specified refant for "
1517 : << typeName() << " solve.)"
1518 : << LogIO::POST;
1519 : refantlist().resize(1);
1520 : refantlist()(0)=-1;
1521 : }
1522 :
1523 : }
1524 :
1525 : void XJones::newselfSolve(VisSet& vs, VisEquation& ve) {
1526 :
1527 : if (prtlev()>4) cout << " Xj::newselfSolve(ve)" << endl;
1528 :
1529 : // Inform logger/history
1530 : logSink() << "Solving for " << typeName()
1531 : << LogIO::POST;
1532 :
1533 : // Initialize the svc according to current VisSet
1534 : // (this counts intervals, sizes CalSet)
1535 : Vector<Int> nChunkPerSol;
1536 : Int nSol = sizeUpSolve(vs,nChunkPerSol);
1537 :
1538 : // Create the Caltable
1539 : createMemCalTable();
1540 :
1541 : // The iterator, VisBuffer
1542 : VisIter& vi(vs.iter());
1543 : VisBuffer vb(vi);
1544 :
1545 : // cout << "nSol = " << nSol << endl;
1546 : // cout << "nChunkPerSol = " << nChunkPerSol << endl;
1547 :
1548 : Vector<Int> slotidx(vs.numberSpw(),-1);
1549 :
1550 : Int nGood(0);
1551 : vi.originChunks();
1552 : for (Int isol=0;isol<nSol && vi.moreChunks();++isol) {
1553 :
1554 : // Arrange to accumulate
1555 : VisBuffAccumulator vba(nAnt(),preavg(),false);
1556 :
1557 : for (Int ichunk=0;ichunk<nChunkPerSol(isol);++ichunk) {
1558 :
1559 : // Current _chunk_'s spw
1560 : Int spw(vi.spectralWindow());
1561 :
1562 : // Abort if we encounter a spw for which a priori cal not available
1563 : if (!ve.spwOK(spw))
1564 : throw(AipsError("Pre-applied calibration not available for at least 1 spw. Check spw selection carefully."));
1565 :
1566 :
1567 : // Collapse each timestamp in this chunk according to VisEq
1568 : // with calibration and averaging
1569 : for (vi.origin(); vi.more(); vi++) {
1570 :
1571 : // Force read of the field Id
1572 : vb.fieldId();
1573 :
1574 : // This forces the data/model/wt I/O, and applies
1575 : // any prior calibrations
1576 : ve.collapse(vb);
1577 :
1578 : // If permitted/required by solvable component, normalize
1579 : if (normalizable())
1580 : vb.normalize();
1581 :
1582 : // If this solve not freqdep, and channels not averaged yet, do so
1583 : if (!freqDepMat() && vb.nChannel()>1)
1584 : vb.freqAveCubes();
1585 :
1586 : // Accumulate collapsed vb in a time average
1587 : vba.accumulate(vb);
1588 : }
1589 : // Advance the VisIter, if possible
1590 : if (vi.moreChunks()) vi.nextChunk();
1591 :
1592 : }
1593 :
1594 : // Finalize the averged VisBuffer
1595 : vba.finalizeAverage();
1596 :
1597 : // The VisBuffer to solve with
1598 : VisBuffer& svb(vba.aveVisBuff());
1599 :
1600 : // Establish meta-data for this interval
1601 : // (some of this may be used _during_ solve)
1602 : // (this sets currSpw() in the SVC)
1603 : Bool vbOk=syncSolveMeta(svb,-1);
1604 :
1605 : Int thisSpw=spwMap()(svb.spectralWindow());
1606 : slotidx(thisSpw)++;
1607 :
1608 : // We are actually solving for all channels simultaneously
1609 : solveCPar().reference(solveAllCPar());
1610 : solveParOK().reference(solveAllParOK());
1611 : solveParErr().reference(solveAllParErr());
1612 : solveParSNR().reference(solveAllParSNR());
1613 :
1614 : // Fill solveCPar() with 1, nominally, and flagged
1615 : solveCPar()=Complex(1.0);
1616 : solveParOK()=false;
1617 :
1618 : if (vbOk && svb.nRow()>0) {
1619 :
1620 : // solve for the R-L phase term in the current VB
1621 : solveOneVB(svb);
1622 :
1623 : if (ntrue(solveParOK())>0) {
1624 : Float ang=arg(sum(solveCPar()(solveParOK()))/Float(ntrue(solveParOK())))*90.0/C::pi;
1625 :
1626 :
1627 : logSink() << "Mean position angle offset solution for "
1628 : << msmc().fieldName(currField())
1629 : << " (spw = " << currSpw() << ") = "
1630 : << ang
1631 : << " deg."
1632 : << LogIO::POST;
1633 : }
1634 : else
1635 : logSink() << "Position angle offset solution for "
1636 : << msmc().fieldName(currField())
1637 : << " (spw = " << currSpw() << ") "
1638 : << " was not determined (insufficient data)."
1639 : << LogIO::POST;
1640 :
1641 : nGood++;
1642 : }
1643 :
1644 : keepNCT();
1645 :
1646 : }
1647 :
1648 : logSink() << " Found good "
1649 : << typeName() << " solutions in "
1650 : << nGood << " intervals."
1651 : << LogIO::POST;
1652 :
1653 : // Store whole of result in a caltable
1654 : if (nGood==0)
1655 : logSink() << "No output calibration table written."
1656 : << LogIO::POST;
1657 : else {
1658 :
1659 : // Do global post-solve tinkering (e.g., phase-only, normalization, etc.)
1660 : // TBD
1661 : // globalPostSolveTinker();
1662 :
1663 : // write the table
1664 : storeNCT();
1665 : }
1666 :
1667 : }
1668 :
1669 :
1670 : void XJones::calcAllJones() {
1671 :
1672 : // cout << "currJElem().shape() = " << currJElem().shape() << endl;
1673 :
1674 : // put the par in the first position on the diagonal
1675 : // [p 0]
1676 : // [0 1]
1677 :
1678 :
1679 : // Set first element to the parameter
1680 : IPosition blc(3,0,0,0);
1681 : IPosition trc(3,0,nChanMat()-1,nElem()-1);
1682 : currJElem()(blc,trc)=currCPar();
1683 : currJElemOK()(blc,trc)=currParOK();
1684 :
1685 : // Set second diag element to one
1686 : blc(0)=trc(0)=1;
1687 : currJElem()(blc,trc)=Complex(1.0);
1688 : currJElemOK()(blc,trc)=currParOK();
1689 :
1690 : }
1691 :
1692 :
1693 : void XJones::solveOneVB(const VisBuffer& vb) {
1694 :
1695 : // This just a simple average of the cross-hand
1696 : // visbilities...
1697 :
1698 : // We are actually solving for all channels simultaneously
1699 : solveCPar().reference(solveAllCPar());
1700 : solveParOK().reference(solveAllParOK());
1701 : solveParErr().reference(solveAllParErr());
1702 : solveParSNR().reference(solveAllParSNR());
1703 :
1704 : // Fill solveCPar() with 1, nominally, and flagged
1705 : solveCPar()=Complex(1.0);
1706 : solveParOK()=false;
1707 :
1708 : Int nChan=vb.nChannel();
1709 :
1710 : Complex d;
1711 : Float wt;
1712 : Vector<DComplex> rl(nChan,0.0),lr(nChan,0.0);
1713 : Double sumwt(0.0);
1714 : for (Int irow=0;irow<vb.nRow();++irow) {
1715 : if (!vb.flagRow()(irow) &&
1716 : vb.antenna1()(irow)!=vb.antenna2()(irow)) {
1717 :
1718 : for (Int ich=0;ich<nChan;++ich) {
1719 : if (!vb.flag()(ich,irow)) {
1720 :
1721 : // A common weight for both crosshands
1722 : // TBD: we should probably consider this carefully...
1723 : // (also in D::guessPar...)
1724 : wt=Double(vb.weightMat()(1,irow)+
1725 : vb.weightMat()(2,irow))/2.0;
1726 :
1727 : // correct weight for model normalization
1728 : // a=abs(vb.modelVisCube()(1,ich,irow));
1729 : // wt*=(a*a);
1730 :
1731 : if (wt>0.0) {
1732 : // Cross-hands only
1733 : for (Int icorr=1;icorr<3;++icorr) {
1734 : // md=vb.modelVisCube()(icorr,ich,irow);
1735 : d=vb.visCube()(icorr,ich,irow);
1736 :
1737 : if (abs(d)>0.0) {
1738 :
1739 : if (icorr==1)
1740 : rl(ich)+=DComplex(Complex(wt)*d);
1741 : // rl(ich)+=DComplex(Complex(wt)*d/md);
1742 : else
1743 : lr(ich)+=DComplex(Complex(wt)*d);
1744 : // lr(ich)+=DComplex(Complex(wt)*d/md);
1745 :
1746 : sumwt+=Double(wt);
1747 :
1748 : } // abs(d)>0
1749 : } // icorr
1750 : } // wt>0
1751 : } // !flag
1752 : } // ich
1753 : } // !flagRow
1754 : } // row
1755 :
1756 :
1757 : // cout << "spw = " << currSpw() << endl;
1758 : // cout << " rl = " << rl << " " << phase(rl)*180.0/C::pi << endl;
1759 : // cout << " lr = " << lr << " " << phase(lr)*180.0/C::pi << endl;
1760 :
1761 : // Record results
1762 : solveCPar()=Complex(1.0);
1763 : solveParOK()=false;
1764 : for (Int ich=0;ich<nChan;++ich) {
1765 : // combine lr with rl
1766 : rl(ich)+=conj(lr(ich));
1767 :
1768 : // Normalize to unit amplitude
1769 : // (note that the phase result is insensitive to sumwt)
1770 : Double amp=abs(rl(ich));
1771 : // For now, all antennas get the same solution
1772 : IPosition blc(3,0,0,0);
1773 : IPosition trc(3,0,0,nElem()-1);
1774 : if (sumwt>0 && amp>0.0) {
1775 : rl(ich)/=DComplex(amp);
1776 : blc(1)=trc(1)=ich;
1777 : solveCPar()(blc,trc)=Complex(rl(ich));
1778 : solveParOK()(blc,trc)=true;
1779 : }
1780 : }
1781 :
1782 :
1783 : if (ntrue(solveParOK())>0) {
1784 : Float ang=arg(sum(solveCPar()(solveParOK()))/Float(ntrue(solveParOK())))*90.0/C::pi;
1785 :
1786 :
1787 : logSink() << "Mean position angle offset solution for "
1788 : << msmc().fieldName(currField())
1789 : << " (spw = " << currSpw() << ") = "
1790 : << ang
1791 : << " deg."
1792 : << LogIO::POST;
1793 : }
1794 : else
1795 : logSink() << "Position angle offset solution for "
1796 : << msmc().fieldName(currField())
1797 : << " (spw = " << currSpw() << ") "
1798 : << " was not determined (insufficient data)."
1799 : << LogIO::POST;
1800 :
1801 : }
1802 :
1803 : void XJones::solveOneSDB(SolveDataBuffer& sdb) {
1804 :
1805 : // This just a simple average of the cross-hand
1806 : // visbilities...
1807 :
1808 : // We are actually solving for all channels simultaneously
1809 : solveCPar().reference(solveAllCPar());
1810 : solveParOK().reference(solveAllParOK());
1811 : solveParErr().reference(solveAllParErr());
1812 : solveParSNR().reference(solveAllParSNR());
1813 :
1814 : // Fill solveCPar() with 1, nominally, and flagged
1815 : solveCPar()=Complex(1.0);
1816 : solveParOK()=false;
1817 :
1818 : Int nChan=sdb.nChannels();
1819 :
1820 : Complex d;
1821 : Float wt;
1822 : Vector<DComplex> rl(nChan,0.0),lr(nChan,0.0);
1823 : Double sumwt(0.0);
1824 : for (Int irow=0;irow<sdb.nRows();++irow) {
1825 : if (!sdb.flagRow()(irow) &&
1826 : sdb.antenna1()(irow)!=sdb.antenna2()(irow)) {
1827 :
1828 : for (Int ich=0;ich<nChan;++ich) {
1829 : if (!sdb.flagCube()(1,ich,irow) &&
1830 : !sdb.flagCube()(2,ich,irow)) {
1831 :
1832 : // A common weight for both crosshands
1833 : // TBD: we should probably consider this carefully...
1834 : // (also in D::guessPar...)
1835 : wt=Double(sdb.weightSpectrum()(1,ich,irow)+
1836 : sdb.weightSpectrum()(2,ich,irow))/2.0;
1837 :
1838 : // correct weight for model normalization
1839 : // a=abs(vb.modelVisCube()(1,ich,irow));
1840 : // wt*=(a*a);
1841 :
1842 : if (wt>0.0) {
1843 : // Cross-hands only
1844 : for (Int icorr=1;icorr<3;++icorr) {
1845 : d=sdb.visCubeCorrected()(icorr,ich,irow);
1846 :
1847 : if (abs(d)>0.0) {
1848 :
1849 : if (icorr==1)
1850 : rl(ich)+=DComplex(Complex(wt)*d);
1851 : else
1852 : lr(ich)+=DComplex(Complex(wt)*d);
1853 :
1854 : sumwt+=Double(wt);
1855 :
1856 : } // abs(d)>0
1857 : } // icorr
1858 : } // wt>0
1859 : } // !flag
1860 : } // ich
1861 : } // !flagRow
1862 : } // row
1863 :
1864 :
1865 : // cout << "spw = " << currSpw() << endl;
1866 : // cout << " rl = " << rl << " " << phase(rl)*180.0/C::pi << endl;
1867 : // cout << " lr = " << lr << " " << phase(lr)*180.0/C::pi << endl;
1868 :
1869 : // Record results
1870 : solveCPar()=Complex(1.0);
1871 : solveParOK()=false;
1872 : for (Int ich=0;ich<nChan;++ich) {
1873 : // combine lr with rl
1874 : rl(ich)+=conj(lr(ich));
1875 :
1876 : // Normalize to unit amplitude
1877 : // (note that the phase result is insensitive to sumwt)
1878 : Double amp=abs(rl(ich));
1879 : // For now, all antennas get the same solution
1880 : IPosition blc(3,0,0,0);
1881 : IPosition trc(3,0,0,nElem()-1);
1882 : if (sumwt>0 && amp>0.0) {
1883 : rl(ich)/=DComplex(amp);
1884 : blc(1)=trc(1)=ich;
1885 : solveCPar()(blc,trc)=Complex(rl(ich));
1886 : solveParOK()(blc,trc)=true;
1887 : }
1888 : }
1889 :
1890 :
1891 : if (ntrue(solveParOK())>0) {
1892 : Float ang=arg(sum(solveCPar()(solveParOK()))/Float(ntrue(solveParOK())))*90.0/C::pi;
1893 :
1894 :
1895 : logSink() << "Mean position angle offset solution for "
1896 : << msmc().fieldName(currField())
1897 : << " (spw = " << currSpw() << ") = "
1898 : << ang
1899 : << " deg."
1900 : << LogIO::POST;
1901 : }
1902 : else
1903 : logSink() << "Position angle offset solution for "
1904 : << msmc().fieldName(currField())
1905 : << " (spw = " << currSpw() << ") "
1906 : << " was not determined (insufficient data)."
1907 : << LogIO::POST;
1908 :
1909 : }
1910 :
1911 :
1912 : void XJones::solveOne(SDBList& sdbs) {
1913 :
1914 : // This just a simple average of the cross-hand
1915 : // visbilities...
1916 :
1917 : Int nSDB=sdbs.nSDB();
1918 :
1919 : //cout << "nSDB=" << nSDB << endl;
1920 :
1921 : // We are actually solving for all channels simultaneously
1922 : solveCPar().reference(solveAllCPar());
1923 : solveParOK().reference(solveAllParOK());
1924 : solveParErr().reference(solveAllParErr());
1925 : solveParSNR().reference(solveAllParSNR());
1926 :
1927 : // Fill solveCPar() with 1, nominally, and flagged
1928 : solveCPar()=Complex(1.0);
1929 : solveParOK()=false;
1930 :
1931 : Int nChan=sdbs.nChannels();
1932 :
1933 : Complex d;
1934 : Float wt;
1935 : Vector<DComplex> rl(nChan,0.0),lr(nChan,0.0);
1936 : Double sumwt(0.0);
1937 : for (Int isdb=0;isdb<nSDB;++isdb) {
1938 : SolveDataBuffer& sdb(sdbs(isdb));
1939 : for (Int irow=0;irow<sdb.nRows();++irow) {
1940 : if (!sdb.flagRow()(irow) &&
1941 : sdb.antenna1()(irow)!=sdb.antenna2()(irow)) {
1942 :
1943 : for (Int ich=0;ich<nChan;++ich) {
1944 : if (!sdb.flagCube()(1,ich,irow) &&
1945 : !sdb.flagCube()(2,ich,irow)) {
1946 :
1947 : // A common weight for both crosshands
1948 : // TBD: we should probably consider this carefully...
1949 : // (also in D::guessPar...)
1950 : wt=Double(sdb.weightSpectrum()(1,ich,irow)+
1951 : sdb.weightSpectrum()(2,ich,irow))/2.0;
1952 :
1953 : // correct weight for model normalization
1954 : // a=abs(vb.modelVisCube()(1,ich,irow));
1955 : // wt*=(a*a);
1956 :
1957 : if (wt>0.0) {
1958 : // Cross-hands only
1959 : for (Int icorr=1;icorr<3;++icorr) {
1960 : d=sdb.visCubeCorrected()(icorr,ich,irow);
1961 :
1962 : if (abs(d)>0.0) {
1963 :
1964 : if (icorr==1)
1965 : rl(ich)+=DComplex(Complex(wt)*d);
1966 : else
1967 : lr(ich)+=DComplex(Complex(wt)*d);
1968 :
1969 : sumwt+=Double(wt);
1970 :
1971 : } // abs(d)>0
1972 : } // icorr
1973 : } // wt>0
1974 : } // !flag
1975 : } // ich
1976 : } // !flagRow
1977 : } // row
1978 : } // isdb
1979 :
1980 : // cout << "spw = " << currSpw() << endl;
1981 : // cout << " rl = " << rl << " " << phase(rl)*180.0/C::pi << endl;
1982 : // cout << " lr = " << lr << " " << phase(lr)*180.0/C::pi << endl;
1983 :
1984 : // Record results
1985 : solveCPar()=Complex(1.0);
1986 : solveParOK()=false;
1987 : for (Int ich=0;ich<nChan;++ich) {
1988 : // combine lr with rl
1989 : rl(ich)+=conj(lr(ich));
1990 :
1991 : // Normalize to unit amplitude
1992 : // (note that the phase result is insensitive to sumwt)
1993 : Double amp=abs(rl(ich));
1994 : // For now, all antennas get the same solution
1995 : IPosition blc(3,0,0,0);
1996 : IPosition trc(3,0,0,nElem()-1);
1997 : if (sumwt>0 && amp>0.0) {
1998 : rl(ich)/=DComplex(amp);
1999 : blc(1)=trc(1)=ich;
2000 : solveCPar()(blc,trc)=Complex(rl(ich));
2001 : solveParOK()(blc,trc)=true;
2002 : }
2003 : }
2004 :
2005 :
2006 : if (ntrue(solveParOK())>0) {
2007 : Float ang=arg(sum(solveCPar()(solveParOK()))/Float(ntrue(solveParOK())))*90.0/C::pi;
2008 :
2009 :
2010 : logSink() << "Mean position angle offset solution for "
2011 : << msmc().fieldName(currField())
2012 : << " (spw = " << currSpw() << ") = "
2013 : << ang
2014 : << " deg."
2015 : << LogIO::POST;
2016 : }
2017 : else
2018 : logSink() << "Position angle offset solution for "
2019 : << msmc().fieldName(currField())
2020 : << " (spw = " << currSpw() << ") "
2021 : << " was not determined (insufficient data)."
2022 : << LogIO::POST;
2023 :
2024 : }
2025 :
2026 : // **********************************************************
2027 : // XfJones: CHANNELIZED position angle for circulars (antenna-based)
2028 : //
2029 :
2030 : XfJones::XfJones(VisSet& vs) :
2031 : VisCal(vs), // virtual base
2032 : VisMueller(vs), // virtual base
2033 : XJones(vs) // immediate parent
2034 : {
2035 : if (prtlev()>2) cout << "Xf::Xf(vs)" << endl;
2036 : }
2037 :
2038 : XfJones::XfJones(String msname,Int MSnAnt,Int MSnSpw) :
2039 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
2040 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
2041 : XJones(msname,MSnAnt,MSnSpw) // immediate parent
2042 : {
2043 : if (prtlev()>2) cout << "Xf::Xf(msname,MSnAnt,MSnSpw)" << endl;
2044 : }
2045 :
2046 : XfJones::XfJones(const MSMetaInfoForCal& msmc) :
2047 : VisCal(msmc), // virtual base
2048 : VisMueller(msmc), // virtual base
2049 : XJones(msmc) // immediate parent
2050 : {
2051 : if (prtlev()>2) cout << "Xf::Xf(msmc)" << endl;
2052 : }
2053 :
2054 : XfJones::XfJones(const Int& nAnt) :
2055 : VisCal(nAnt),
2056 : VisMueller(nAnt),
2057 : XJones(nAnt)
2058 : {
2059 : if (prtlev()>2) cout << "Xf::Xf(nAnt)" << endl;
2060 : }
2061 :
2062 : XfJones::~XfJones() {
2063 : if (prtlev()>2) cout << "Xf::~Xf()" << endl;
2064 : }
2065 :
2066 : void XfJones::initSolvePar() {
2067 :
2068 : SolvableVisJones::initSolvePar();
2069 : return;
2070 :
2071 : }
2072 :
2073 :
2074 :
2075 :
2076 : // **********************************************************
2077 : // GlinXphJones Implementations
2078 : //
2079 :
2080 : GlinXphJones::GlinXphJones(VisSet& vs) :
2081 : VisCal(vs), // virtual base
2082 : VisMueller(vs), // virtual base
2083 : GJones(vs), // immediate parent
2084 : QU_()
2085 : {
2086 : if (prtlev()>2) cout << "GlinXph::GlinXph(vs)" << endl;
2087 : }
2088 :
2089 : GlinXphJones::GlinXphJones(String msname,Int MSnAnt,Int MSnSpw) :
2090 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
2091 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
2092 : GJones(msname,MSnAnt,MSnSpw), // immediate parent
2093 : QU_()
2094 : {
2095 : if (prtlev()>2) cout << "GlinXph::GlinXph(msname,MSnAnt,MSnSpw)" << endl;
2096 : }
2097 :
2098 : GlinXphJones::GlinXphJones(const MSMetaInfoForCal& msmc) :
2099 : VisCal(msmc), // virtual base
2100 : VisMueller(msmc), // virtual base
2101 : GJones(msmc), // immediate parent
2102 : QU_()
2103 : {
2104 : if (prtlev()>2) cout << "GlinXph::GlinXph(msmc)" << endl;
2105 : }
2106 :
2107 : GlinXphJones::GlinXphJones(const Int& nAnt) :
2108 : VisCal(nAnt),
2109 : VisMueller(nAnt),
2110 : GJones(nAnt),
2111 : QU_()
2112 : {
2113 : if (prtlev()>2) cout << "GlinXph::GlinXph(nAnt)" << endl;
2114 : }
2115 :
2116 : GlinXphJones::~GlinXphJones() {
2117 : if (prtlev()>2) cout << "GlinXph::~GlinXph()" << endl;
2118 : }
2119 :
2120 : void GlinXphJones::setApply(const Record& apply) {
2121 :
2122 : GJones::setApply(apply);
2123 :
2124 : // Force calwt to false
2125 : calWt()=false;
2126 :
2127 : }
2128 :
2129 :
2130 :
2131 : void GlinXphJones::selfGatherAndSolve(VisSet& vs, VisEquation& ve) {
2132 :
2133 : if (prtlev()>4) cout << " GlnXph::selfGatherAndSolve(ve)" << endl;
2134 :
2135 : // Inform logger/history
2136 : logSink() << "Solving for " << typeName()
2137 : << LogIO::POST;
2138 :
2139 : // Initialize the svc according to current VisSet
2140 : // (this counts intervals, sizes CalSet)
2141 : Vector<Int> nChunkPerSol;
2142 : Int nSol = sizeUpSolve(vs,nChunkPerSol);
2143 :
2144 : // Create the Caltable
2145 : createMemCalTable();
2146 :
2147 : // The iterator, VisBuffer
2148 : VisIter& vi(vs.iter());
2149 : VisBuffer vb(vi);
2150 :
2151 : // cout << "nSol = " << nSol << endl;
2152 : // cout << "nChunkPerSol = " << nChunkPerSol << endl;
2153 :
2154 : Int nGood(0);
2155 : vi.originChunks();
2156 : for (Int isol=0;isol<nSol && vi.moreChunks();++isol) {
2157 :
2158 : // Arrange to accumulate
2159 : VisBuffAccumulator vba(nAnt(),preavg(),false);
2160 :
2161 : for (Int ichunk=0;ichunk<nChunkPerSol(isol);++ichunk) {
2162 :
2163 : // Current _chunk_'s spw
2164 : Int spw(vi.spectralWindow());
2165 :
2166 : // Abort if we encounter a spw for which a priori cal not available
2167 : if (!ve.spwOK(spw))
2168 : throw(AipsError("Pre-applied calibration not available for at least 1 spw. Check spw selection carefully."));
2169 :
2170 :
2171 : // Collapse each timestamp in this chunk according to VisEq
2172 : // with calibration and averaging
2173 : for (vi.origin(); vi.more(); vi++) {
2174 :
2175 : // Force read of the field Id
2176 : vb.fieldId();
2177 :
2178 : // This forces the data/model/wt I/O, and applies
2179 : // any prior calibrations
2180 : ve.collapse(vb);
2181 :
2182 : // If permitted/required by solvable component, normalize
2183 : if (normalizable())
2184 : vb.normalize();
2185 :
2186 : // If this solve not freqdep, and channels not averaged yet, do so
2187 : if (!freqDepMat() && vb.nChannel()>1)
2188 : vb.freqAveCubes();
2189 :
2190 : // Accumulate collapsed vb in a time average
2191 : vba.accumulate(vb);
2192 : }
2193 : // Advance the VisIter, if possible
2194 : if (vi.moreChunks()) vi.nextChunk();
2195 :
2196 : }
2197 :
2198 : // Finalize the averged VisBuffer
2199 : vba.finalizeAverage();
2200 :
2201 : // The VisBuffer to solve with
2202 : VisBuffer& svb(vba.aveVisBuff());
2203 :
2204 : // Establish meta-data for this interval
2205 : // (some of this may be used _during_ solve)
2206 : // (this sets currSpw() in the SVC)
2207 : Bool vbOk=syncSolveMeta(svb,-1);
2208 :
2209 : if (vbOk && svb.nRow()>0) {
2210 :
2211 : // solve for the X-Y phase term in the current VB
2212 : solveOneVB(svb);
2213 :
2214 : nGood++;
2215 : }
2216 :
2217 : keepNCT();
2218 :
2219 : }
2220 :
2221 : logSink() << " Found good "
2222 : << typeName() << " solutions in "
2223 : << nGood << " intervals."
2224 : << LogIO::POST;
2225 :
2226 : // Store whole of result in a caltable
2227 : if (nGood==0)
2228 : logSink() << "No output calibration table written."
2229 : << LogIO::POST;
2230 : else {
2231 :
2232 : // Do global post-solve tinkering (e.g., phase-only, normalization, etc.)
2233 : globalPostSolveTinker();
2234 :
2235 : // write the table
2236 : storeNCT();
2237 : }
2238 :
2239 : }
2240 :
2241 : // Handle trivial vbga
2242 : void GlinXphJones::selfSolveOne(VisBuffGroupAcc& vbga) {
2243 :
2244 : // Expecting only on VB in the vbga (with many times)
2245 : if (vbga.nBuf()!=1)
2246 : throw(AipsError("GlinXphJones can't process multi-vb vbga."));
2247 :
2248 : // Call single-VB specialized solver with the one vb
2249 : this->solveOneVB(vbga(0));
2250 :
2251 : }
2252 :
2253 : // SDBList (VI2) version
2254 : void GlinXphJones::selfSolveOne(SDBList& sdbs) {
2255 :
2256 : // Expecting multiple SDBs (esp. in time)
2257 : if (sdbs.nSDB()==1)
2258 : throw(AipsError("GlinXphJones needs multiple SDBs"));
2259 :
2260 : // Call single-VB specialized solver with the one vb
2261 : this->solveOne(sdbs);
2262 :
2263 : }
2264 :
2265 : // Solve for the X-Y phase from the cross-hand's slope in R/I
2266 : void GlinXphJones::solveOneVB(const VisBuffer& vb) {
2267 :
2268 : // ensure
2269 : if (QU_.shape()!=IPosition(2,2,nSpw())) {
2270 : QU_.resize(2,nSpw());
2271 : QU_.set(0.0);
2272 : }
2273 :
2274 : Int thisSpw=spwMap()(vb.spectralWindow());
2275 :
2276 : // We are actually solving for all channels simultaneously
2277 : solveCPar().reference(solveAllCPar());
2278 : solveParOK().reference(solveAllParOK());
2279 : solveParErr().reference(solveAllParErr());
2280 : solveParSNR().reference(solveAllParSNR());
2281 :
2282 : // Fill solveCPar() with 1, nominally, and flagged
2283 : solveCPar()=Complex(1.0);
2284 : solveParOK()=false;
2285 :
2286 : Int nChan=vb.nChannel();
2287 : // if (nChan>1)
2288 : // throw(AipsError("X-Y phase solution NYI for channelized data"));
2289 :
2290 : // Find number of timestamps in the VB
2291 : Vector<uInt> ord;
2292 : Int nTime=genSort(ord,vb.time(),Sort::Ascending,Sort::NoDuplicates);
2293 :
2294 : Matrix<Double> x(nTime,nChan,0.0),y(nTime,nChan,0.0),wt(nTime,nChan,0.0),sig(nTime,nChan,0.0);
2295 : Matrix<Bool> mask(nTime,nChan,false);
2296 :
2297 : mask.set(false);
2298 : Complex v(0.0);
2299 : Float wt0(0.0);
2300 : Int iTime(-1);
2301 : Double currtime(-1.0);
2302 : for (Int irow=0;irow<vb.nRow();++irow) {
2303 : if (!vb.flagRow()(irow) &&
2304 : vb.antenna1()(irow)!=vb.antenna2()(irow)) {
2305 :
2306 : // Advance time index when we see a new time
2307 : if (vb.time()(irow)!=currtime) {
2308 : ++iTime;
2309 : currtime=vb.time()(irow); // remember the new current time
2310 : }
2311 :
2312 : // Weights not yet chan-dep
2313 : wt0=(vb.weightMat()(1,irow)+vb.weightMat()(2,irow));
2314 : if (wt0>0.0) {
2315 :
2316 : for (Int ich=0;ich<nChan;++ich) {
2317 : if (!vb.flag()(ich,irow)) {
2318 : v=vb.visCube()(1,ich,irow)+conj(vb.visCube()(2,ich,irow));
2319 : x(iTime,ich)+=Double(wt0*real(v));
2320 : y(iTime,ich)+=Double(wt0*imag(v));
2321 : wt(iTime,ich)+=Double(wt0);
2322 : }
2323 : }
2324 : }
2325 : }
2326 : }
2327 :
2328 : // Normalize data by accumulated weights
2329 : for (Int itime=0;itime<nTime;++itime) {
2330 : for (Int ich=0;ich<nChan;++ich) {
2331 : if (wt(itime,ich)>0.0) {
2332 : x(itime,ich)/=wt(itime,ich);
2333 : y(itime,ich)/=wt(itime,ich);
2334 : sig(itime,ich)=sqrt(1.0/wt(itime,ich));
2335 : mask(itime,ich)=true;
2336 : }
2337 : else
2338 : sig(itime,ich)=DBL_MAX; // ~zero weight
2339 : }
2340 : }
2341 :
2342 : // Solve for each channel
2343 : Vector<Complex> Cph(nChan);
2344 : Cph.set(Complex(1.0,0.0));
2345 : Float currAmb(1.0);
2346 : Bool report(false);
2347 : for (Int ich=0;ich<nChan;++ich) {
2348 :
2349 : if (sum(wt.column(ich))>0.0) {
2350 : report=true;
2351 : LinearFit<Double> phfitter;
2352 : Polynomial<AutoDiff<Double> > line(1);
2353 : phfitter.setFunction(line);
2354 : Vector<Bool> m(mask.column(ich));
2355 :
2356 : // Fit shallow and steep, and always prefer shallow
2357 :
2358 : // Assumed shallow solve:
2359 : Vector<Double> solnA;
2360 : solnA.assign(phfitter.fit(x.column(ich),y.column(ich),sig.column(ich),&m));
2361 :
2362 : // Assumed steep solve:
2363 : Vector<Double> solnB;
2364 : solnB.assign(phfitter.fit(y.column(ich),x.column(ich),sig.column(ich),&m));
2365 :
2366 : Double Xph(0.0);
2367 : if (abs(solnA(1))<abs(solnB(1))) {
2368 : Xph=atan(solnA(1));
2369 : }
2370 : else {
2371 : Xph=atan(1.0/solnB(1));
2372 : }
2373 :
2374 : Cph(ich)=currAmb*Complex(DComplex(cos(Xph),sin(Xph)));
2375 :
2376 : // Watch for and remove ambiguity changes which can
2377 : // occur near Xph~= +/-90 deg (the atan above can jump)
2378 : // (NB: this does _not_ resolve the amb; it merely makes
2379 : // it consistent within the spw)
2380 : if (ich>0) {
2381 : // If Xph changes by more than pi/2, probably a ambig jump...
2382 : Float dang=abs(arg(Cph(ich)/Cph(ich-1)));
2383 : if (dang > (C::pi/2.)) {
2384 : Cph(ich)*=-1.0; // fix this one
2385 : currAmb*=-1.0; // reverse currAmb, so curr amb is carried forward
2386 : // cout << " Found XY phase ambiguity jump at chan=" << ich << " in spw=" << currSpw();
2387 : }
2388 : }
2389 :
2390 : // cout << " (" << currAmb << ")";
2391 : // cout << endl;
2392 :
2393 :
2394 : // Set all antennas with this X-Y phase (as a complex number)
2395 : solveCPar()(Slice(0,1,1),Slice(ich,1,1),Slice())=Cph(ich);
2396 : solveParOK()(Slice(),Slice(ich,1,1),Slice())=true;
2397 : }
2398 : else {
2399 : solveCPar()(Slice(0,1,1),Slice(ich,1,1),Slice())=Complex(1.0);
2400 : solveParOK()(Slice(),Slice(ich,1,1),Slice())=false;
2401 : }
2402 : }
2403 :
2404 : if (report)
2405 : cout << endl
2406 : << "Spw = " << thisSpw
2407 : << " (ich=" << nChan/2 << "/" << nChan << "): " << endl
2408 : << " X-Y phase = " << arg(Cph[nChan/2])*180.0/C::pi << " deg." << endl;
2409 :
2410 :
2411 : // Now fit for the source polarization
2412 : {
2413 :
2414 : Vector<Double> wtf(nTime,0.0),sigf(nTime,0.0),xf(nTime,0.0),yf(nTime,0.0);
2415 : Vector<Bool> maskf(nTime,false);
2416 : Float wt0;
2417 : Complex v;
2418 : Double currtime(-1.0);
2419 : Int iTime(-1);
2420 : for (Int irow=0;irow<vb.nRow();++irow) {
2421 : if (!vb.flagRow()(irow) &&
2422 : vb.antenna1()(irow)!=vb.antenna2()(irow)) {
2423 :
2424 : if (vb.time()(irow)!=currtime) {
2425 : // Advance time index when we see a new time
2426 : ++iTime;
2427 : currtime=vb.time()(irow); // remember the new current time
2428 : }
2429 :
2430 : // Weights not yet chan-dep
2431 : wt0=(vb.weightMat()(1,irow)+vb.weightMat()(2,irow));
2432 : if (wt0>0.0) {
2433 : for (Int ich=0;ich<nChan;++ich) {
2434 :
2435 : if (!vb.flag()(ich,irow)) {
2436 : // Correct x-hands for xy-phase and add together
2437 : v=vb.visCube()(1,ich,irow)/Cph(ich)+vb.visCube()(2,ich,irow)/conj(Cph(ich));
2438 : xf(iTime)+=Double(wt0*2.0*(vb.feed_pa(vb.time()(irow))(0)));
2439 : yf(iTime)+=Double(wt0*real(v)/2.0);
2440 : wtf(iTime)+=Double(wt0);
2441 : }
2442 : }
2443 : }
2444 : }
2445 : }
2446 :
2447 : // Normalize data by accumulated weights
2448 : for (Int itime=0;itime<nTime;++itime) {
2449 : if (wtf(itime)>0.0) {
2450 : xf(itime)/=wtf(itime);
2451 : yf(itime)/=wtf(itime);
2452 : sigf(itime)=sqrt(1.0/wtf(itime));
2453 : maskf(itime)=true;
2454 : }
2455 : else
2456 : sigf(itime)=DBL_MAX; // ~zero weight
2457 : }
2458 :
2459 : // p0=Q, p1=U, p2 = real part of net instr pol offset
2460 : // x is TWICE the parallactic angle
2461 : CompiledFunction<AutoDiff<Double> > fn;
2462 : fn.setFunction("-p0*sin(x) + p1*cos(x) + p2");
2463 :
2464 : LinearFit<Double> fitter;
2465 : fitter.setFunction(fn);
2466 :
2467 : Vector<Double> soln=fitter.fit(xf,yf,sigf,&maskf);
2468 :
2469 : QU_(0,thisSpw) = soln(0);
2470 : QU_(1,thisSpw) = soln(1);
2471 :
2472 : cout << " Fractional Poln: "
2473 : << "Q = " << QU_(0,thisSpw) << ", "
2474 : << "U = " << QU_(1,thisSpw) << "; "
2475 : << "P = " << sqrt(soln(0)*soln(0)+soln(1)*soln(1)) << ", "
2476 : << "X = " << atan2(soln(1),soln(0))*90.0/C::pi << "deg."
2477 : << endl;
2478 : cout << " Net (over baselines) instrumental polarization: "
2479 : << soln(2) << endl;
2480 :
2481 : }
2482 :
2483 : }
2484 :
2485 : // Solve for the X-Y phase from the cross-hand's slope in R/I
2486 : void GlinXphJones::solveOne(SDBList& sdbs) {
2487 :
2488 : // ensure
2489 : if (QU_.shape()!=IPosition(2,2,nSpw())) {
2490 : QU_.resize(2,nSpw());
2491 : QU_.set(0.0);
2492 : }
2493 :
2494 : Int thisSpw=sdbs.aggregateSpw();
2495 :
2496 : // We are actually solving for all channels simultaneously
2497 : solveCPar().reference(solveAllCPar());
2498 : solveParOK().reference(solveAllParOK());
2499 : solveParErr().reference(solveAllParErr());
2500 : solveParSNR().reference(solveAllParSNR());
2501 :
2502 : // Fill solveCPar() with 1, nominally, and flagged
2503 : solveCPar()=Complex(1.0);
2504 : solveParOK()=false;
2505 :
2506 : Int nChan=sdbs.nChannels();
2507 :
2508 : // Number of datapoints in fit is the number of SDBs
2509 : Int nSDB=sdbs.nSDB();
2510 :
2511 : Matrix<Double> x(nSDB,nChan,0.0),y(nSDB,nChan,0.0),wt(nSDB,nChan,0.0),sig(nSDB,nChan,0.0);
2512 : Matrix<Bool> mask(nSDB,nChan,false);
2513 :
2514 : mask.set(false);
2515 : Complex v(0.0);
2516 : Float wt0(0.0);
2517 :
2518 : for (Int isdb=0;isdb<nSDB;++isdb) {
2519 : SolveDataBuffer& sdb(sdbs(isdb));
2520 :
2521 : for (Int irow=0;irow<sdb.nRows();++irow) {
2522 : if (!sdb.flagRow()(irow) &&
2523 : sdb.antenna1()(irow)!=sdb.antenna2()(irow)) {
2524 :
2525 : for (Int ich=0;ich<nChan;++ich) {
2526 : if (!sdb.flagCube()(1,ich,irow)) {
2527 : wt0=sdb.weightSpectrum()(1,ich,irow);
2528 : v=sdb.visCubeCorrected()(1,ich,irow);
2529 : x(isdb,ich)+=Double(wt0*real(v));
2530 : y(isdb,ich)+=Double(wt0*imag(v));
2531 : wt(isdb,ich)+=Double(wt0);
2532 : }
2533 : if (!sdb.flagCube()(2,ich,irow)) {
2534 : wt0=sdb.weightSpectrum()(2,ich,irow);
2535 : v=conj(sdb.visCubeCorrected()(2,ich,irow));
2536 : x(isdb,ich)+=Double(wt0*real(v));
2537 : y(isdb,ich)+=Double(wt0*imag(v));
2538 : wt(isdb,ich)+=Double(wt0);
2539 : }
2540 : }
2541 : }
2542 : }
2543 : }
2544 :
2545 : // Normalize data by accumulated weights
2546 : for (Int isdb=0;isdb<nSDB;++isdb) {
2547 : for (Int ich=0;ich<nChan;++ich) {
2548 : if (wt(isdb,ich)>0.0) {
2549 : x(isdb,ich)/=wt(isdb,ich);
2550 : y(isdb,ich)/=wt(isdb,ich);
2551 : sig(isdb,ich)=sqrt(1.0/wt(isdb,ich));
2552 : mask(isdb,ich)=true;
2553 : }
2554 : else
2555 : sig(isdb,ich)=DBL_MAX; // ~zero weight
2556 : }
2557 : }
2558 :
2559 : // Solve for each channel
2560 : Vector<Complex> Cph(nChan);
2561 : Cph.set(Complex(1.0,0.0));
2562 : Float currAmb(1.0);
2563 : Bool report(false);
2564 : for (Int ich=0;ich<nChan;++ich) {
2565 :
2566 : if (sum(wt.column(ich))>0.0) {
2567 : report=true;
2568 : LinearFit<Double> phfitter;
2569 : Polynomial<AutoDiff<Double> > line(1);
2570 : phfitter.setFunction(line);
2571 : Vector<Bool> m(mask.column(ich));
2572 :
2573 : // Fit shallow and steep, and always prefer shallow
2574 :
2575 : // Assumed shallow solve:
2576 : Vector<Double> solnA;
2577 : solnA.assign(phfitter.fit(x.column(ich),y.column(ich),sig.column(ich),&m));
2578 :
2579 : // Assumed steep solve:
2580 : Vector<Double> solnB;
2581 : solnB.assign(phfitter.fit(y.column(ich),x.column(ich),sig.column(ich),&m));
2582 :
2583 : Double Xph(0.0);
2584 : if (abs(solnA(1))<abs(solnB(1))) {
2585 : Xph=atan(solnA(1));
2586 : }
2587 : else {
2588 : Xph=atan(1.0/solnB(1));
2589 : }
2590 :
2591 : Cph(ich)=currAmb*Complex(DComplex(cos(Xph),sin(Xph)));
2592 :
2593 : // Watch for and remove ambiguity changes which can
2594 : // occur near Xph~= +/-90 deg (the atan above can jump)
2595 : // (NB: this does _not_ resolve the amb; it merely makes
2596 : // it consistent within the spw)
2597 : if (ich>0) {
2598 : // If Xph changes by more than pi/2, probably a ambig jump...
2599 : Float dang=abs(arg(Cph(ich)/Cph(ich-1)));
2600 : if (dang > (C::pi/2.)) {
2601 : Cph(ich)*=-1.0; // fix this one
2602 : currAmb*=-1.0; // reverse currAmb, so curr amb is carried forward
2603 : // cout << " Found XY phase ambiguity jump at chan=" << ich << " in spw=" << currSpw();
2604 : }
2605 : }
2606 :
2607 : // cout << " (" << currAmb << ")";
2608 : // cout << endl;
2609 :
2610 :
2611 : // Set all antennas with this X-Y phase (as a complex number)
2612 : solveCPar()(Slice(0,1,1),Slice(ich,1,1),Slice())=Cph(ich);
2613 : solveParOK()(Slice(),Slice(ich,1,1),Slice())=true;
2614 : }
2615 : else {
2616 : solveCPar()(Slice(0,1,1),Slice(ich,1,1),Slice())=Complex(1.0);
2617 : solveParOK()(Slice(),Slice(ich,1,1),Slice())=false;
2618 : }
2619 : }
2620 :
2621 : if (report)
2622 : cout << endl
2623 : << "Spw = " << thisSpw
2624 : << " (ich=" << nChan/2 << "/" << nChan << "): " << endl
2625 : << " X-Y phase = " << arg(Cph[nChan/2])*180.0/C::pi << " deg." << endl;
2626 :
2627 :
2628 : // Now fit for the source polarization
2629 : {
2630 :
2631 : Vector<Double> wtf(nSDB,0.0),sigf(nSDB,0.0),xf(nSDB,0.0),yf(nSDB,0.0);
2632 : Vector<Bool> maskf(nSDB,false);
2633 : Float wt0;
2634 : Complex v;
2635 : for (Int isdb=0;isdb<nSDB;++isdb) {
2636 : SolveDataBuffer& sdb(sdbs(isdb));
2637 :
2638 : for (Int irow=0;irow<sdb.nRows();++irow) {
2639 : if (!sdb.flagRow()(irow) &&
2640 : sdb.antenna1()(irow)!=sdb.antenna2()(irow)) {
2641 :
2642 : Float fpa(sdb.feedPa()(0)); // assumes same for all antennas!
2643 :
2644 : for (Int ich=0;ich<nChan;++ich) {
2645 :
2646 : if (!sdb.flagCube()(1,ich,irow)) {
2647 : // Correct x-hands for xy-phase and add together
2648 : wt0=sdb.weightSpectrum()(1,ich,irow);
2649 : v=sdb.visCubeCorrected()(1,ich,irow)/Cph(ich);
2650 : xf(isdb)+=Double(wt0*2.0*(fpa));
2651 : yf(isdb)+=Double(wt0*real(v));
2652 : wtf(isdb)+=Double(wt0);
2653 : }
2654 : if (!sdb.flagCube()(2,ich,irow)) {
2655 : // Correct x-hands for xy-phase and add together
2656 : wt0=sdb.weightSpectrum()(2,ich,irow);
2657 : v=sdb.visCubeCorrected()(2,ich,irow)/conj(Cph(ich));
2658 : xf(isdb)+=Double(wt0*2.0*(fpa));
2659 : yf(isdb)+=Double(wt0*real(v));
2660 : wtf(isdb)+=Double(wt0);
2661 : }
2662 : }
2663 : }
2664 : }
2665 : }
2666 :
2667 : // Normalize data by accumulated weights
2668 : for (Int isdb=0;isdb<nSDB;++isdb) {
2669 : if (wtf(isdb)>0.0) {
2670 : xf(isdb)/=wtf(isdb);
2671 : yf(isdb)/=wtf(isdb);
2672 : sigf(isdb)=sqrt(1.0/wtf(isdb));
2673 : maskf(isdb)=true;
2674 : }
2675 : else
2676 : sigf(isdb)=DBL_MAX; // ~zero weight
2677 : }
2678 :
2679 : // p0=Q, p1=U, p2 = real part of net instr pol offset
2680 : // x is TWICE the parallactic angle
2681 : CompiledFunction<AutoDiff<Double> > fn;
2682 : fn.setFunction("-p0*sin(x) + p1*cos(x) + p2");
2683 :
2684 : LinearFit<Double> fitter;
2685 : fitter.setFunction(fn);
2686 :
2687 : Vector<Double> soln=fitter.fit(xf,yf,sigf,&maskf);
2688 :
2689 : srcPolPar().resize(2);
2690 : srcPolPar()(0)=soln(0);
2691 : srcPolPar()(1)=soln(1);
2692 :
2693 : QU_(0,thisSpw) = soln(0);
2694 : QU_(1,thisSpw) = soln(1);
2695 :
2696 : cout << " Fractional Poln: "
2697 : << "Q = " << QU_(0,thisSpw) << ", "
2698 : << "U = " << QU_(1,thisSpw) << "; "
2699 : << "P = " << sqrt(soln(0)*soln(0)+soln(1)*soln(1)) << ", "
2700 : << "X = " << atan2(soln(1),soln(0))*90.0/C::pi << "deg."
2701 : << endl;
2702 : cout << " Net (over baselines) instrumental polarization: "
2703 : << soln(2) << endl;
2704 :
2705 : }
2706 :
2707 : }
2708 :
2709 : void GlinXphJones::globalPostSolveTinker() {
2710 :
2711 : // Add QU info the the keywords
2712 : TableRecord& tr(ct_->rwKeywordSet());
2713 : Record qu;
2714 : qu.define("QU",QU_);
2715 : tr.defineRecord("QU",qu);
2716 :
2717 : }
2718 :
2719 :
2720 : // **********************************************************
2721 : // GlinXphfJones Implementations
2722 : //
2723 :
2724 : // Constructor
2725 : GlinXphfJones::GlinXphfJones(VisSet& vs) :
2726 : VisCal(vs), // virtual base
2727 : VisMueller(vs), // virtual base
2728 : GlinXphJones(vs) // immediate parent
2729 : {
2730 : if (prtlev()>2) cout << "GlinXphf::GlinXphf(vs)" << endl;
2731 : }
2732 :
2733 : GlinXphfJones::GlinXphfJones(String msname,Int MSnAnt,Int MSnSpw) :
2734 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
2735 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
2736 : GlinXphJones(msname,MSnAnt,MSnSpw) // immediate parent
2737 : {
2738 : if (prtlev()>2) cout << "GlinXphf::GlinXphf(msname,MSnAnt,MSnSpw)" << endl;
2739 : }
2740 :
2741 : GlinXphfJones::GlinXphfJones(const MSMetaInfoForCal& msmc) :
2742 : VisCal(msmc), // virtual base
2743 : VisMueller(msmc), // virtual base
2744 : GlinXphJones(msmc) // immediate parent
2745 : {
2746 : if (prtlev()>2) cout << "GlinXphf::GlinXphf(msmc)" << endl;
2747 : }
2748 :
2749 : GlinXphfJones::GlinXphfJones(const Int& nAnt) :
2750 : VisCal(nAnt),
2751 : VisMueller(nAnt),
2752 : GlinXphJones(nAnt)
2753 : {
2754 : if (prtlev()>2) cout << "GlinXphf::GlinXphf(nAnt)" << endl;
2755 : }
2756 :
2757 : GlinXphfJones::~GlinXphfJones() {
2758 : if (prtlev()>2) cout << "GlinXphf::~GlinXphf()" << endl;
2759 : }
2760 :
2761 : */
2762 :
2763 :
2764 :
2765 : } //# NAMESPACE CASA - END
|
