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