Line data Source code
1 : //# BJonesPoly.cc: Implementation of BJonesPoly.h
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id: BJonesPoly.cc,v 19.18 2006/02/03 00:29:52 gmoellen Exp $
27 :
28 : #include <synthesis/MeasurementComponents/BPoly.h>
29 : #include <synthesis/MeasurementComponents/MSMetaInfoForCal.h>
30 : #include <synthesis/MeasurementEquations/VisEquation.h>
31 : #include <casacore/casa/Logging/LogIO.h>
32 : #include <casacore/casa/Utilities/Assert.h>
33 : #include <casacore/casa/Exceptions/Error.h>
34 : #include <casacore/casa/Arrays/ArrayMath.h>
35 : #include <casacore/casa/Arrays/ArrayLogical.h>
36 : #include <casacore/ms/MSOper/MSMetaData.h>
37 : #include <casacore/ms/MeasurementSets/MSColumns.h>
38 : #include <msvis/MSVis/VisBuffAccumulator.h>
39 : #include <msvis/MSVis/VisBuffGroupAcc.h>
40 : #include <sstream>
41 : #include <cmath>
42 : #include <casacore/casa/OS/Memory.h>
43 : #include <casacore/casa/System/PGPlotter.h>
44 : #include <synthesis/CalTables/BJonesMBuf.h>
45 : #include <synthesis/CalTables/BJonesMCol.h>
46 : #include <synthesis/CalTables/NewCalTable.h>
47 : #include <casacore/ms/MSSel/MSSpWindowIndex.h>
48 : #include <casacore/tables/Tables/TableUtil.h>
49 :
50 : using namespace casacore;
51 : namespace casa { //# NAMESPACE CASA - BEGIN
52 :
53 : // Define external CLIC solvers
54 : #define NEED_UNDERSCORES
55 : #if defined(NEED_UNDERSCORES)
56 : #define polyant polyant_
57 : #define splinant splinant_
58 : #define getbspl getbspl_
59 : #define phaseant phaseant_
60 : #define ampliant ampliant_
61 : #define cheb cheb_
62 : #endif
63 :
64 : extern "C" {
65 : void polyant(Int*, Int*, Int*, Int*, Int*, Int*, Int*, Int*,
66 : Double*, Double*, Double*, Double*, Double*, Double*,
67 : Double*, Double* );
68 : void splinant(Int*, Int*, Int*, Int*, Int*, Int*, Int*, Int*,
69 : Double*, Double*, Double*, Double*, Double*, Double*,
70 : Double*, Double*, Double* );
71 :
72 : void getbspl(Int*, Double*, Double*, Double*, Double*, Int*);
73 :
74 : void phaseant(Int*, Int*, Int*, Int*, Double*, Double*, Double*, Double*,
75 : Int*, Double*, Double*);
76 :
77 : void ampliant(Int*, Int*, Int*, Int*, Double*, Double*, Double*, Double*,
78 : Int*, Double*, Double*);
79 :
80 : void cheb(Int*, Double*, Double*, Int*);
81 : }
82 :
83 : //----------------------------------------------------------------------------
84 :
85 0 : BJonesPoly::BJonesPoly (VisSet& vs) :
86 : VisCal(vs), // virtual base
87 : VisMueller(vs), // virtual base
88 : BJones(vs), // immediate base
89 : vs_p(&vs),
90 : degamp_p(3),
91 : degphase_p(3),
92 : visnorm_p(false),
93 : maskcenter_p(1),
94 : maskedge_p(5.0),
95 : maskcenterHalf_p(0),
96 : maskedgeFrac_p(0.05),
97 : solTimeStamp(0.0),
98 : solSpwId(-1),
99 0 : solFldId(-1)
100 : {
101 : // Construct from a visibility set
102 : // Input:
103 : // vs VisSet& Visibility set
104 : // Output to private data:
105 : // vs_p VisSet& Visibility set pointer (needed locally)
106 : // degamp_p Int Polynomial degree in amplitude
107 : // degphase_p Int Polynomial degree in phase
108 : // visnorm const Bool& true if pre-normalization of the
109 : // visibility data over frequency is
110 : // required before solving.
111 : // maskcenter const Int& No. of central channels to mask
112 : // during the solution
113 : // maskedge const Float& Fraction of spectrum to mask at
114 : // either edge during solution
115 : // maskcenterHalf_p Int Central mask half-width
116 : // maskedgeFrac_p Float Fractional edge mask
117 :
118 : // Neither solved nor applied at this point
119 0 : setSolved(false);
120 0 : setApplied(false);
121 0 : };
122 :
123 : //----------------------------------------------------------------------------
124 :
125 :
126 :
127 0 : BJonesPoly::BJonesPoly (const MSMetaInfoForCal& msmc) :
128 : VisCal(msmc), // virtual base
129 : VisMueller(msmc), // virtual base
130 : BJones(msmc), // immediate base
131 : vs_p(NULL),
132 : degamp_p(3),
133 : degphase_p(3),
134 : visnorm_p(false),
135 : maskcenter_p(1),
136 : maskedge_p(5.0),
137 : maskcenterHalf_p(0),
138 : maskedgeFrac_p(0.05),
139 : solTimeStamp(0.0),
140 : solSpwId(-1),
141 0 : solFldId(-1)
142 : {
143 : // Construct from a visibility set
144 : // Input:
145 : // msmc MSMetaInfoForCal& MS Meta info server
146 : // Output to private data:
147 : // vs_p VisSet& Visibility set pointer (needed locally)
148 : // degamp_p Int Polynomial degree in amplitude
149 : // degphase_p Int Polynomial degree in phase
150 : // visnorm const Bool& true if pre-normalization of the
151 : // visibility data over frequency is
152 : // required before solving.
153 : // maskcenter const Int& No. of central channels to mask
154 : // during the solution
155 : // maskedge const Float& Fraction of spectrum to mask at
156 : // either edge during solution
157 : // maskcenterHalf_p Int Central mask half-width
158 : // maskedgeFrac_p Float Fractional edge mask
159 :
160 : // Neither solved nor applied at this point
161 0 : setSolved(false);
162 0 : setApplied(false);
163 0 : };
164 :
165 :
166 :
167 :
168 :
169 0 : void BJonesPoly::setSolve(const Record& solvepar)
170 : {
171 : // Set the solver parameters
172 : // Input:
173 : // solvepar const Record& Solver parameters
174 : // Output to private data:
175 : // degamp_p Int Polynomial degree in amplitude
176 : // degphase_p Int Polynomial degree in phase
177 : // visnorm const Bool& true if pre-normalization of the
178 : // visibility data over frequency is
179 : // required before solving.
180 : // maskcenter const Int& No. of central channels to mask
181 : // during the solution
182 : // maskedge const Float& Fraction of spectrum to mask at
183 : // either edge during solution
184 : // maskcenterHalf_p Int Central mask half-width
185 : // maskedgeFrac_p Float Fractional edge mask
186 : //
187 :
188 : // Call parent
189 0 : SolvableVisCal::setSolve(solvepar);
190 :
191 : // Delete calTableName() if it exists and we are not appending
192 : // TBD: move to SVC?
193 0 : if (!append()) {
194 0 : if (TableUtil::canDeleteTable(calTableName())) {
195 0 : TableUtil::deleteTable(calTableName());
196 : }
197 : // else
198 : // throw(AipsError(calTableName()+" exists and can't delete! (plotting or browsing it?)"));
199 : }
200 :
201 : // Extract the BPOLY-specific solve parameters
202 0 : Vector<Int> degree;
203 0 : if (solvepar.isDefined("degree")) degree = solvepar.asArrayInt("degree");
204 0 : degamp_p = degree(0);
205 0 : degphase_p = degree(1);
206 0 : if (solvepar.isDefined("visnorm")) visnorm_p = solvepar.asBool("visnorm");
207 0 : if (solvepar.isDefined("maskcenter")) maskcenter_p =
208 0 : solvepar.asInt("maskcenter");
209 0 : if (solvepar.isDefined("maskedge")) maskedge_p = solvepar.asFloat("maskedge");
210 :
211 : // Compute derived mask parameters
212 0 : maskcenterHalf_p = maskcenter_p / 2;
213 0 : maskedgeFrac_p = maskedge_p / 100.0;
214 :
215 0 : preavg()=DBL_MAX;
216 :
217 0 : };
218 :
219 : //----------------------------------------------------------------------------
220 :
221 0 : void BJonesPoly::setApply(const Record& applypar)
222 : {
223 : // Set the applypar parameters
224 :
225 : // Call parent (loadMemCalTable is now overloaded)
226 :
227 : // nChanParList()=vs_p->numberChan();
228 : // startChanList()=vs_p->startChan();
229 :
230 0 : BJones::setApply(applypar);
231 :
232 : // The old BPOLY never used calWt=true; preserve
233 : // this behavior for now
234 0 : calWt()=false;
235 :
236 : /*
237 :
238 : // Extract the parameters
239 : if (applypar.isDefined("table"))
240 : calTableName() = applypar.asString("table");
241 : if (applypar.isDefined("select"))
242 : calTableSelect() = applypar.asString("select");
243 : if (applypar.isDefined("t"))
244 : interval() = applypar.asDouble("t");
245 : if (applypar.isDefined("interp"))
246 : tInterpType()=applypar.asString("interp");
247 :
248 : // Protect against non-specific interp
249 : if (tInterpType()=="")
250 : tInterpType()="linear";
251 :
252 : indgen(spwMap());
253 : if (applypar.isDefined("spwmap")) {
254 : Vector<Int> spwmap(applypar.asArrayInt("spwmap"));
255 : if (allGE(spwmap,0)) {
256 : // User has specified a valid spwmap
257 : if (spwmap.nelements()==1)
258 : spwMap()=spwmap(0);
259 : else
260 : spwMap()(IPosition(1,0),IPosition(1,spwmap.nelements()-1))=spwmap;
261 : // TBD: Report non-trivial spwmap to logger.
262 : // cout << "Note: spwMap() = " << spwMap() << endl;
263 : }
264 : }
265 : AlwaysAssert(allGE(spwMap(),0),AipsError);
266 :
267 : // cout << "BPOLY: spwMap() = " << spwMap()<< endl;
268 :
269 : // Follow the selected data, for the moment
270 : // NB: this requires setdata prior to setapply!
271 : nChanParList()=vs_p->numberChan();
272 : startChanList()=vs_p->startChan();
273 :
274 : // If selection is non-trivial, complain for now
275 : if (calTableSelect()!="") {
276 : LogIO os (LogOrigin("BJonesPoly", "setApply()", WHERE));
277 :
278 : os << LogIO::WARN
279 : << "Selection on BPOLY caltables not yet supported."
280 : << LogIO::POST;
281 : }
282 :
283 : // Fill the bandpass correction cache from the applied cal. table
284 : load(calTableName());
285 :
286 : // Signal apply context
287 : setApplied(true);
288 : setSolved(false);
289 :
290 : */
291 :
292 0 : };
293 :
294 : //----------------------------------------------------------------------------
295 :
296 :
297 0 : void BJonesPoly::selfSolveOne(VisBuffGroupAcc& vbga)
298 : {
299 : // Solver for the polynomial bandpass solution
300 :
301 : // TODO:
302 : // 1. Make pointers private, make delete function and use it
303 : // 2. Use antenna names
304 : //
305 :
306 :
307 0 : LogIO os (LogOrigin("BJonesPoly", "selfSolveOne()", WHERE));
308 :
309 : os << LogIO::NORMAL
310 : << "THIS IS THE NEW MULTI-SPW-FLEXIBLE VERSION"
311 0 : << LogIO::POST;
312 :
313 : os << LogIO::NORMAL
314 : << "Fitting bandpass amplitude and phase polynomials."
315 0 : << LogIO::POST;
316 : os << LogIO::NORMAL
317 : << "Polynomial degree for amplitude is " << degamp_p
318 0 : << LogIO::POST;
319 : os << LogIO::NORMAL
320 : << "Polynomial degree for phase is " << degphase_p
321 0 : << LogIO::POST;
322 :
323 : // Bool to short-circuit operation
324 0 : Bool ok(true);
325 :
326 : // Initialize the baseline index
327 0 : Vector<Int> ant1(nBln(), -1);
328 0 : Vector<Int> ant2(nBln(), -1);
329 0 : Vector<Int> numFreqChan(nSpw(), 0);
330 :
331 : // make gridded freq array
332 : // The minimum number of frequency channels required for a solution
333 : // is the number of coefficients in the fit. For a gridded spectrum,
334 : // filled from irregularly spaced input spectral windows, it is possible
335 : // that only very few channels get filled. So, we will be conservative
336 : // and make a gridded spectrum with 2*ncoeff channels, where ncoeff is
337 : // the maximum of the number of coefficients requested for the phase and
338 : // amp solutions. We will then check to make sure that a sufficient
339 : // number of gridded slots will be filled by the input frequency channels.
340 :
341 :
342 0 : Int nPH(0);
343 0 : Double minfreq(DBL_MAX), maxfreq(0.0), maxdf(0.0);
344 0 : PtrBlock<Vector<Double>*> freq; freq.resize(nSpw()); freq=0;
345 :
346 0 : for (Int ibuf=0;ibuf<vbga.nBuf();++ibuf) {
347 :
348 0 : CalVisBuffer& cvb(vbga(ibuf));
349 :
350 0 : Int spwid=cvb.spectralWindow();
351 0 : numFreqChan(spwid) = cvb.nChannel();
352 0 : freq[spwid] = new Vector<Double>;
353 0 : (*freq[spwid])=cvb.frequency();
354 0 : Double df2=abs((*freq[spwid])(1)-(*freq[spwid])(0))/2.0;
355 0 : maxdf=max(maxdf,2.0*df2);
356 0 : minfreq=min(minfreq,min((*freq[spwid])));
357 0 : maxfreq=max(maxfreq,max((*freq[spwid])));
358 :
359 0 : nPH= max(nPH,min(cvb.nCorr(),2));
360 : }
361 0 : minfreq=minfreq-maxdf/2.0;
362 0 : maxfreq=maxfreq+maxdf/2.0;
363 :
364 : // minfreq is the low edge of the lowest channel
365 : // maxfreq is the high edge of the highest channel
366 : // nPH is the number of parallel-hand correlations present
367 :
368 : // cout << "freq, corr: " << minfreq << " " << maxfreq << " " << nPH << endl;
369 :
370 : Double freqstep;
371 : Int nFreqGrid;
372 :
373 : // Derive grid spacing/number from requested poly degree
374 : // FOR NOW, it is less error-prone to use input spacing
375 : // nFreqGrid=2*(max(degamp_p,degphase_p)+1);
376 : // nFreqGrid=max(nFreqGrid,16); // no fewer than 16, in any case.
377 : // freqstep=((maxfreq-minfreq)/Double(nFreqGrid));
378 :
379 : // Grid spacing is (multiple of?) maximum input channel spacing
380 0 : freqstep=maxdf;
381 0 : nFreqGrid = Int ((maxfreq-minfreq)/freqstep+0.5);
382 :
383 : // Fill the gridded frequency list
384 0 : Vector<Double> totalFreq(nFreqGrid,0.0);
385 0 : for (Int i=0;i<nFreqGrid;i++) {
386 0 : totalFreq(i)=minfreq + freqstep*(Double(i)+0.5);
387 : }
388 :
389 : // Populate the frequency grid with the input frequency channels,
390 : // and demand that enough (?) will get filled.
391 0 : Vector<Bool> freqGridOk(nFreqGrid,false);
392 0 : for (Int ispw=0;ispw<nSpw();ispw++) {
393 0 : if (freq[ispw]) {
394 0 : for (Int ichan=0;ichan<numFreqChan(ispw);ichan++) {
395 0 : Int chanidx=(Int) floor(((*freq[ispw])(ichan)-minfreq)/freqstep);
396 0 : if(chanidx >= nFreqGrid){
397 0 : cout << "spw " << ispw <<" " << chanidx << endl;
398 : }
399 0 : freqGridOk(chanidx)=true;
400 : }
401 : }
402 : }
403 :
404 : // Sanity check on polynomial order and grid count
405 0 : Int nok=ntrue(freqGridOk);
406 0 : if (nok < (degamp_p+1) ) {
407 : os << LogIO::SEVERE
408 : << "Selected spectral window(s) nominally fill only " << nok << " grid points."
409 0 : << LogIO::POST;
410 : os << LogIO::SEVERE
411 0 : << "Reduce degamp by at least " << degamp_p+1-nok << " and try again."
412 0 : << LogIO::POST;
413 0 : ok=false;
414 : }
415 0 : if (nok < (degphase_p+1) ) {
416 : os << LogIO::SEVERE
417 : << "Selected spectral window(s) nominally fill only " << nok << " grid points."
418 0 : << LogIO::POST;
419 : os << LogIO::SEVERE
420 0 : << "Reduce degphase by at least " << degphase_p+1-nok << " and try again."
421 0 : << LogIO::POST;
422 0 : ok=false;
423 : }
424 : // If either degree vs nGrid test failed, quit here
425 0 : if (!ok) throw(AipsError("Invalid polynomial degree specification"));
426 :
427 : // Report spectral information
428 : os << LogIO::NORMAL
429 : << "Spectral grid for fit will have " << nFreqGrid
430 : << " points spaced by " << freqstep/1000.0 << " kHz."
431 0 : << LogIO::POST;
432 :
433 : os << LogIO::NORMAL
434 : << "Polynomial solution will be valid over frequency range: "
435 0 : << totalFreq(0) << "-" << totalFreq(nFreqGrid-1) << " Hz."
436 0 : << LogIO::POST;
437 :
438 : os << LogIO::NORMAL
439 : << "Total bandwidth: "
440 0 : << freqstep*Double(nFreqGrid)/1000.0 << " kHz."
441 0 : << LogIO::POST;
442 :
443 :
444 : // We must keep track of good ants and baselines
445 0 : Vector<Bool> antok(nAnt(),false);
446 0 : Matrix<Bool> bslok(nAnt(),nAnt(),false);
447 0 : Int nGoodBasl=0;
448 0 : Matrix<Int> bslidx(nAnt(),nAnt(),-1);
449 0 : Matrix<Int> antOkChan(nFreqGrid, nAnt(),0);
450 0 : Vector<Int> ant1num, ant2num; // clic solver antenna numbers (1-based, contiguous)
451 0 : ant1num.resize(nBln());
452 0 : ant2num.resize(nBln());
453 0 : Vector<Int> ant1idx, ant2idx; // MS.ANTENNA indices (for plot, storage)
454 0 : ant1idx.resize(nBln());
455 0 : ant2idx.resize(nBln());
456 :
457 0 : PtrBlock<Matrix<Complex>*> accumVis(nPH,NULL);
458 0 : PtrBlock<Matrix<Double>*> accumWeight(nPH,NULL);
459 0 : PtrBlock<Vector<Complex>*> normVis(nPH,NULL);
460 0 : PtrBlock<Vector<Double>*> normWeight(nPH,NULL);
461 :
462 0 : for (Int i=0;i<nPH;++i) {
463 0 : accumVis[i] = new Matrix<Complex>(nFreqGrid, nBln()); (*accumVis[i])=Complex(0.0,0.0);
464 0 : accumWeight[i] = new Matrix<Double>(nFreqGrid, nBln()); (*accumWeight[i])=0.0;
465 0 : normVis[i] = new Vector<Complex>(nBln()); (*normVis[i])=Complex(0.0,0.0);
466 0 : normWeight[i] = new Vector<Double>(nBln()); (*normWeight[i])=0.0;
467 : }
468 :
469 0 : Vector<Int> indexSpw;
470 :
471 : // By constraint, this solver should see data only from one sideband.
472 : // Pick up the frequency group name from the current spectral window
473 : // accordingly.
474 0 : String freqGroup("");
475 :
476 0 : solTimeStamp=refTime();
477 0 : solFldId=currField();
478 0 : solSpwId=currSpw();
479 :
480 : // Filter data from VBs to FORTRAN arrays for the solver
481 0 : for (Int ibuf=0;ibuf<vbga.nBuf();++ibuf) {
482 :
483 0 : CalVisBuffer& cvb(vbga(ibuf));
484 :
485 : // Extract the current visibility buffer spectral window id.
486 : // and number for frequency channels
487 :
488 0 : Int spwid = cvb.spectralWindow();
489 0 : Int nChan = cvb.nChannel();
490 0 : Int nCorr = cvb.nCorr();
491 0 : freqGroup = freqGrpName(spwid);
492 :
493 0 : Vector<Int> polidx(2,0);
494 0 : polidx(1)=nCorr-1; // TBD: should be pol-sensitive!
495 :
496 : // Data and model values
497 0 : Complex data;
498 :
499 : // Compute the amplitude and phase spectrum for this spectral window
500 0 : for (Int row=0; row < cvb.nRow(); row++) {
501 : // Antenna indices
502 0 : Int a1 = cvb.antenna1()(row);
503 0 : Int a2 = cvb.antenna2()(row);
504 0 : Vector<Double> rowwt(2,0.0);
505 0 : for (Int iph=0;iph<nPH;++iph)
506 0 : rowwt(iph) = cvb.weightMat()(polidx(iph),row);
507 :
508 : // Reject auto-correlation data, zero weights, fully-flagged spectra
509 0 : if ( (a1 != a2) &&
510 0 : (sum(rowwt) > 0) &&
511 0 : nfalse(cvb.flag().column(row)) > 0 ) {
512 :
513 : // These ants, this baseline is ok (there is at least some good data here)
514 0 : antok(a1)=true;
515 0 : antok(a2)=true;
516 0 : if (!bslok(a1,a2)) { // first visit to this baseline
517 0 : bslok(a1,a2)=true;
518 0 : bslidx(a1,a2)=nGoodBasl++;
519 0 : ant1idx(bslidx(a1,a2))=a1;
520 0 : ant2idx(bslidx(a1,a2))=a2;
521 : }
522 :
523 : // Loop over the frequency channels
524 0 : for (Int ichan=0; ichan < nChan; ichan++) {
525 : // Reject flagged/masked channels
526 0 : if (!cvb.flag()(ichan,row) && !maskedChannel(ichan, nChan)) {
527 :
528 0 : for (Int iph=0;iph<nPH;++iph) {
529 :
530 0 : data = cvb.visCube()(polidx(iph),ichan,row);
531 :
532 : // accumulate accumVis, accumWeight
533 0 : if (abs(data) > 0 && rowwt(iph)>0.0) {
534 0 : Vector<Double> freq1=cvb.frequency();
535 0 : Int chanidx=(Int) floor((freq1(ichan)-minfreq)/freqstep);
536 0 : (*accumVis[iph])(chanidx,bslidx(a1,a2))+=(Complex(rowwt(iph),0.0)*data);
537 0 : (*accumWeight[iph])(chanidx,bslidx(a1,a2))+=rowwt(iph);
538 :
539 : // count this channel good for these ants
540 0 : antOkChan(chanidx,a1)++;
541 0 : antOkChan(chanidx,a2)++;
542 :
543 : // Accumulate data normalizing factor (visnorm)
544 0 : (*normVis[iph])(bslidx(a1,a2))+=(Complex(rowwt(iph),0.0)*data);
545 0 : (*normWeight[iph])(bslidx(a1,a2))+=rowwt(iph);
546 :
547 : // cout << row << " " << ichan << " " << chanidx << " "
548 : // << data << " " << rowwt << " "
549 : // << accumVis(chanidx,bslidx(a1,a2)) << " "
550 : // << accumWeight(chanidx,bslidx(a1,a2)) << endl;
551 : }
552 :
553 : }
554 :
555 : };
556 : };
557 : };
558 : };
559 : }; // for (ibuf...)
560 :
561 : // cout << "solFldId = " << solFldId << endl;
562 : // cout << "solSpwId = " << solSpwId << endl;
563 : // cout << "solTimeStamp = " << MVTime(solTimeStamp/C::day).string( MVTime::YMD,7) << endl;
564 :
565 : // Delete freq PtrBlock
566 0 : for (Int ispw=0;ispw<nSpw();ispw++) {
567 0 : if (freq[ispw]) { delete freq[ispw]; freq[ispw]=0; }
568 : }
569 :
570 : // Form a contiguous 1-based antenna index list
571 0 : Vector<Int> antnum(nAnt(),-1);
572 0 : Int antcount=0;
573 0 : Int refantenna(-1);
574 0 : for (Int iant=0;iant<nAnt();iant++) {
575 0 : if (antok(iant)) {
576 0 : antnum(iant)=(++antcount);
577 :
578 : // detect reference antenna with revised counting
579 0 : if (iant==refant())
580 0 : refantenna = antnum(iant);
581 : }
582 : }
583 :
584 : // Use first antenna as refant if requested one not ok
585 0 : if (refantenna<1)
586 0 : refantenna=1;
587 :
588 : // cout << boolalpha << "antok = " << antok << endl;
589 : // cout << "antnum = " << antnum << endl;
590 :
591 0 : for (Int ibl=0;ibl<nGoodBasl;ibl++) {
592 0 : ant1num(ibl)=antnum(ant1idx(ibl));
593 0 : ant2num(ibl)=antnum(ant2idx(ibl));
594 : }
595 0 : ant1num.resize(nGoodBasl,true);
596 0 : ant2num.resize(nGoodBasl,true);
597 0 : ant1idx.resize(nGoodBasl,true);
598 0 : ant2idx.resize(nGoodBasl,true);
599 :
600 :
601 0 : Int nGoodAnt=Int(ntrue(antok));
602 :
603 : os << LogIO::NORMAL
604 : << "Found data for " << nGoodBasl
605 : << " baselines among " << nGoodAnt
606 : << " antennas."
607 0 : << LogIO::POST;
608 :
609 : // cout << "antOkChan = " << antOkChan << endl;
610 :
611 : // BPoly is nominally dual-pol
612 0 : Matrix<Double> ampCoeff(nAnt(), 2*(degamp_p+1),0.0); // solutions stored here later
613 0 : Matrix<Double> phaseCoeff(nAnt(), 2*(degphase_p+1), 0.0); // solutions stored here later
614 :
615 0 : PtrBlock<Matrix<Double>*> totalWeight(nPH,NULL), totalPhase(nPH,NULL), totalAmp(nPH,NULL);
616 :
617 0 : for (Int iph=0;iph<nPH;++iph) {
618 :
619 0 : totalAmp[iph] = new Matrix<Double>(nFreqGrid, nGoodBasl); (*totalAmp[iph])=0.0;
620 0 : totalPhase[iph] = new Matrix<Double>(nFreqGrid, nGoodBasl); (*totalPhase[iph])=0.0;
621 0 : totalWeight[iph] = new Matrix<Double>(nFreqGrid, nGoodBasl); (*totalWeight[iph])=0.0;
622 :
623 0 : for (Int ibl=0;ibl<nGoodBasl;ibl++) {
624 0 : ant1num(ibl)=antnum(ant1idx(ibl));
625 0 : ant2num(ibl)=antnum(ant2idx(ibl));
626 0 : if ( (*normWeight[iph])(ibl)> 0.0)
627 0 : (*normVis[iph])(ibl)/=(static_cast<Complex>((*normWeight[iph])(ibl)));
628 0 : for (Int ichan=0;ichan<nFreqGrid;ichan++) {
629 0 : Double &wt=(*accumWeight[iph])(ichan,ibl);
630 : // insist at least 2 baselines with good data for these antennas in this channel
631 0 : if (wt > 0 &&
632 0 : antOkChan(ichan,ant1idx(ibl)) > 1 &&
633 0 : antOkChan(ichan,ant2idx(ibl)) > 1 ) {
634 0 : (*accumVis[iph])(ichan,ibl)/= (static_cast<Complex>(wt));
635 :
636 : // If requested, normalize the data, if possible
637 0 : if (visnorm_p) {
638 0 : if (abs((*normVis[iph])(ibl))>0.0 )
639 0 : (*accumVis[iph])(ichan,ibl)/=(*normVis[iph])(ibl);
640 : else
641 0 : (*accumVis[iph])(ichan,ibl)=0.0; // causes problems in polyant?
642 : }
643 :
644 0 : (*totalAmp[iph])(ichan,ibl)=static_cast<Double>(log(abs((*accumVis[iph])(ichan,ibl))));
645 : // Note phase sign convention!
646 0 : (*totalPhase[iph])(ichan,ibl)=static_cast<Double>(-1.0*arg((*accumVis[iph])(ichan,ibl)));
647 0 : (*totalWeight[iph])(ichan,ibl)=wt;
648 : }
649 : }
650 :
651 : }
652 :
653 0 : if ( accumVis[iph] ) delete accumVis[iph];
654 0 : if ( accumWeight[iph] ) delete accumWeight[iph];
655 0 : if ( normVis[iph] ) delete normVis[iph];
656 0 : if ( normWeight[iph] ) delete normWeight[iph];
657 0 : accumVis[iph]=NULL;
658 0 : accumWeight[iph]=NULL;
659 0 : normVis[iph]=NULL;
660 0 : normWeight[iph]=NULL;
661 :
662 : // First fit the bandpass amplitude
663 : os << LogIO::NORMAL
664 : << "Fitting amplitude polynomial."
665 0 : << LogIO::POST;
666 :
667 0 : Int degree = degamp_p + 1;
668 0 : Int iy = 1;
669 : Bool dum;
670 0 : Vector<Double> rmsAmpFit2(nGoodBasl,0.0);
671 0 : Matrix<Double> ampCoeff2(nGoodAnt, degree, 1.0);
672 :
673 : {
674 : // Create work arrays
675 0 : Vector<Double> wk1(degree);
676 0 : Vector<Double> wk2(degree*degree*nGoodAnt*nGoodAnt);
677 0 : Vector<Double> wk3(degree*nGoodAnt);
678 :
679 : // Call the CLIC solver for amplitude
680 0 : polyant(&iy,
681 : &nFreqGrid,
682 : &nGoodBasl,
683 : ant1num.getStorage(dum),
684 : ant2num.getStorage(dum),
685 : &refantenna,
686 : °ree,
687 : &nGoodAnt,
688 : totalFreq.getStorage(dum),
689 0 : totalAmp[iph]->getStorage(dum),
690 0 : totalWeight[iph]->getStorage(dum),
691 : wk1.getStorage(dum),
692 : wk2.getStorage(dum),
693 : wk3.getStorage(dum),
694 : rmsAmpFit2.getStorage(dum),
695 : ampCoeff2.getStorage(dum));
696 : }
697 :
698 0 : if (totalAmp[iph]) delete totalAmp[iph];
699 0 : totalAmp[iph]=NULL;
700 :
701 : os << LogIO::NORMAL
702 : << "Per-baseline RMS log(Amp) statistics: (min/mean/max) = "
703 : << min(rmsAmpFit2) << "/" << mean(rmsAmpFit2) << "/" << max(rmsAmpFit2)
704 0 : << LogIO::POST;
705 :
706 : // Now fit the bandpass phase
707 : os << LogIO::NORMAL
708 : << "Fitting phase polynomial."
709 0 : << LogIO::POST;
710 :
711 : // Call the CLIC solver for phase
712 0 : degree = degphase_p + 1;
713 0 : iy = 2;
714 0 : Vector<Double> rmsPhaseFit2(nGoodBasl,0.0);
715 0 : Matrix<Double> phaseCoeff2(nGoodAnt, degree, 123.0);
716 :
717 : {
718 : // Create work arrays
719 0 : Vector<Double> wk6(degree);
720 0 : Vector<Double> wk7(degree*degree*nGoodAnt*nGoodAnt);
721 0 : Vector<Double> wk8(degree*nGoodAnt);
722 :
723 0 : polyant(&iy,
724 : &nFreqGrid,
725 : &nGoodBasl,
726 : ant1num.getStorage(dum),
727 : ant2num.getStorage(dum),
728 : &refantenna,
729 : °ree,
730 : &nGoodAnt,
731 : totalFreq.getStorage(dum),
732 0 : totalPhase[iph]->getStorage(dum),
733 0 : totalWeight[iph]->getStorage(dum),
734 : wk6.getStorage(dum),
735 : wk7.getStorage(dum),
736 : wk8.getStorage(dum),
737 : rmsPhaseFit2.getStorage(dum),
738 : phaseCoeff2.getStorage(dum));
739 : }
740 :
741 0 : if (totalPhase[iph]) delete totalPhase[iph];
742 0 : if (totalWeight[iph]) delete totalWeight[iph];
743 0 : totalPhase[iph]=NULL;
744 0 : totalWeight[iph]=NULL;
745 :
746 : os << LogIO::NORMAL
747 : << "Per baseline RMS phase (deg) statistics: (min/mean/max) = "
748 0 : << min(rmsPhaseFit2)*180.0/C::pi << "/"
749 0 : << mean(rmsPhaseFit2)*180.0/C::pi << "/"
750 0 : << max(rmsPhaseFit2)*180.0/C::pi
751 0 : << LogIO::POST;
752 :
753 : // Expand solutions into full antenna list
754 0 : IPosition iplo(1,iph*(degphase_p+1));
755 0 : IPosition iphi(1,(iph+1)*(degphase_p+1)-1);
756 0 : IPosition ialo(1,iph*(degamp_p+1));
757 0 : IPosition iahi(1,(iph+1)*(degamp_p+1)-1);
758 :
759 0 : for (Int iant=0;iant<nAnt();iant++) {
760 0 : if (antok(iant)) {
761 0 : ampCoeff.row(iant)(ialo,iahi)=ampCoeff2.row(antnum(iant)-1);
762 0 : phaseCoeff.row(iant)(iplo,iphi)=phaseCoeff2.row(antnum(iant)-1);
763 : }
764 : }
765 :
766 : // plot amplitude and phase baseline data/solutions
767 : // plotsolve2(totalFreq, totalAmp, totalPhase, totalWeight, ant1idx, ant2idx,
768 : // rmsAmpFit2, ampCoeff, rmsPhaseFit2, phaseCoeff);
769 :
770 : } // iph
771 :
772 0 : Int nChanTotal=nFreqGrid;
773 0 : Double meanfreq=mean(totalFreq);
774 :
775 : // Compute the reference frequency and reference phasor
776 0 : Vector<Complex> refAmp;
777 0 : Vector<MFrequency> refFreq;
778 0 : refAmp.resize(nAnt());
779 0 : refFreq.resize(nAnt());
780 0 : for (Int k=0; k < nAnt(); k++) {
781 0 : refAmp(k) = Complex(1.0,0);
782 0 : refFreq(k) = MFrequency(Quantity(meanfreq, "Hz"));
783 : };
784 :
785 : // Frequency range
786 0 : Double loFreq = totalFreq(0);
787 0 : Double hiFreq = totalFreq(nChanTotal-1);
788 :
789 0 : Matrix<Double> validDomain(nAnt(), 2);
790 0 : validDomain.column(0) = loFreq;
791 0 : validDomain.column(1) = hiFreq;
792 : // TBD:
793 : // Double edgefreq(maskedge_p*(hiFreq-loFreq));
794 : // validDomain.column(0) = loFreq + edgefreq;
795 : // validDomain.column(1) = hiFreq - edgefreq;
796 :
797 : // Normalize the output calibration solutions if required
798 0 : Vector<Complex> scaleFactor(nAnt(), Complex(1,0));
799 :
800 0 : if (solnorm()) {
801 : os << LogIO::NORMAL
802 : << "Normalizing antenna-based solutions."
803 0 : << LogIO::POST;
804 :
805 0 : Vector<Double> Ca, Cp;
806 0 : for (Int iph=0;iph<nPH;++iph) {
807 :
808 0 : IPosition iplo(1,iph*(degphase_p+1));
809 0 : IPosition iphi(1,(iph+1)*(degphase_p+1)-1);
810 0 : IPosition ialo(1,iph*(degamp_p+1));
811 0 : IPosition iahi(1,(iph+1)*(degamp_p+1)-1);
812 :
813 : // Loop over antenna
814 0 : for (Int iant=0; iant < nAnt(); iant++) {
815 :
816 0 : Ca.reference(ampCoeff.row(iant)(ialo,iahi));
817 0 : Cp.reference(phaseCoeff.row(iant)(iplo,iphi));
818 :
819 :
820 : // Normalize mean phasor across the spectrum
821 0 : Int nChAve(0);
822 0 : Complex meanPha(0,0);
823 0 : Double meanAmp(0);
824 0 : for (Int chan=0; chan < nChanTotal; chan++) {
825 : // only use channels that participated in the fit
826 0 : if (antOkChan(chan,iant) > 3) {
827 :
828 0 : nChAve++;
829 0 : Double ampval = getChebVal(Ca, loFreq, hiFreq,
830 0 : totalFreq(chan));
831 0 : Double phaseval = getChebVal(Cp, loFreq, hiFreq,
832 0 : totalFreq(chan));
833 :
834 0 : meanPha += Complex(cos(phaseval), sin(phaseval));
835 0 : meanAmp += exp(ampval);
836 :
837 : };
838 : }
839 : // Normalize by adjusting the zeroth order term
840 0 : if (nChAve>0) {
841 0 : meanPha/=Complex(nChAve);
842 0 : meanAmp/=Double(nChAve);
843 :
844 : // cout << "mean B = " << meanAmp << " " << arg(meanPha)*180.0/C::pi << endl;
845 :
846 0 : Ca(0)-=2.0*log(meanAmp);
847 0 : Cp(0)-=2.0*arg(meanPha);
848 : }
849 : };
850 : }
851 : };
852 :
853 : // Update the output calibration table
854 0 : Vector<Int> antId(nAnt());
855 0 : indgen(antId);
856 0 : Vector<String> polyType(nAnt(), "CHEBYSHEV");
857 0 : Vector<String> phaseUnits(nAnt(), "RADIANS");
858 0 : Vector<Int> refAnt(nAnt(), refant());
859 :
860 0 : updateCalTable (freqGroup, antId, polyType, scaleFactor, validDomain,
861 : ampCoeff, phaseCoeff, phaseUnits, refAmp, refFreq, refAnt);
862 :
863 : // After first call, append must be true in case we have more to
864 : // write.
865 0 : append()=true;
866 :
867 0 : };
868 :
869 :
870 : // "Calculate" current parameters
871 0 : void BJonesPoly::calcPar() {
872 :
873 : // Currently, we only support a single bandpass solution
874 : // (i.e., no time-dep), so if any currParOK() (for currSpw),
875 : // then we have a good solution
876 :
877 : // If any
878 : // if (sum(currParOK())>0 )
879 : // validateP();
880 : // else
881 : // throw(AipsError("No calibration available for current Spw!"));
882 :
883 : // Call parent
884 0 : BJones::calcPar();
885 :
886 0 : }
887 :
888 :
889 : //----------------------------------------------------------------------------
890 :
891 0 : void BJonesPoly::updateCalTable (const String& freqGrpName,
892 : const Vector<Int>& antennaId,
893 : const Vector<String>& polyType,
894 : const Vector<Complex>& scaleFactor,
895 : const Matrix<Double>& validDomain,
896 : const Matrix<Double>& polyCoeffAmp,
897 : const Matrix<Double>& polyCoeffPhase,
898 : const Vector<String>& phaseUnits,
899 : const Vector<Complex>& sideBandRef,
900 : const Vector<MFrequency>& refFreq,
901 : const Vector<Int>& refAnt)
902 : {
903 : // Update the output calibration table to include the current soln. parameters
904 : // Input:
905 : // freqGrpName const String& Freq. group name
906 : // antennaId const Vector<Int>& Antenna id. for each soln.
907 : // polyType const Vector<String>& Polynomial type (per soln.)
908 : // scaleFactor const Vector<Complex>& Scale factor (per soln.)
909 : // validDomain const Matrix<Double>& Valid domain [x_0, x_1]
910 : // (per solution)
911 : // polyCoeffAmp const Matrix<Double>& Polynomial amplitude
912 : // coefficients (per soln.)
913 : // polyCoeffPhase const Matrix<Double>& Polynomial phase coefficients
914 : // (per solution)
915 : // phaseUnits const Vector<String>& Phase units (per solution)
916 : // sideBandRef const Vector<Complex>& Sideband reference phasor
917 : // (per solution)
918 : // refFreq const Vec<MFrequency>& Sideband reference frequency
919 : // (per solution)
920 : // refAnt const Vector<Int>& Reference antenna (per soln.)
921 : // Input from private data:
922 : // solveTable_p String Output calibration table name
923 : //
924 :
925 0 : LogIO os (LogOrigin("BJonesPoly", "updateCalTable()", WHERE));
926 :
927 : // Open the caltable
928 0 : Table::TableOption accessMode = Table::New;
929 0 : if (append() && Table::isWritable(calTableName())) {
930 0 : accessMode = Table::Update;
931 : }
932 0 : BJonesPolyTable calTable(calTableName(), accessMode);
933 :
934 0 : Int nDesc=calTable.nRowDesc();
935 0 : Int idesc=0;
936 0 : Int thisDesc=-1;
937 0 : while (idesc<nDesc && thisDesc<0) {
938 0 : Vector<Int> spw;
939 0 : CalDescRecord calDescRec(calTable.getRowDesc(idesc));
940 0 : calDescRec.getSpwId(spw);
941 :
942 0 : if (spw(0)==solSpwId)
943 0 : thisDesc=idesc;
944 :
945 0 : ++idesc;
946 : }
947 0 : if (thisDesc<0)
948 0 : thisDesc=nDesc;
949 :
950 : // Fill the bandpass solution parameters to a BJonesPoly calibration
951 : // buffer spanning the antenna id.'s
952 0 : Vector<Int> key(1, MSCalEnums::ANTENNA1);
953 0 : Block<Vector<Int> > keyvals(1, antennaId);
954 0 : BJonesPolyMBuf buffer(key, keyvals);
955 :
956 : // Add each solution to the calibration buffer
957 0 : for (uInt k=0; k < antennaId.nelements(); k++) {
958 0 : buffer.putAntGain(antennaId(k), freqGrpName, polyType(k), scaleFactor(k),
959 0 : validDomain.row(k),
960 : // polyCoeffAmp.row(k).nelements(),
961 : // polyCoeffPhase.row(k).nelements(),
962 0 : degamp_p+1,degphase_p+1,
963 0 : polyCoeffAmp.row(k), polyCoeffPhase.row(k),
964 0 : phaseUnits(k), sideBandRef(k), refFreq(k), refAnt(k));
965 : };
966 :
967 0 : buffer.fieldId().set(solFldId);
968 0 : buffer.calDescId().set(thisDesc);
969 0 : buffer.timeMeas().set(MEpoch(MVEpoch(solTimeStamp/86400.0)));
970 :
971 : os << LogIO::NORMAL
972 0 : << "Storing calibration in table " << calTableName()
973 0 : << LogIO::POST;
974 :
975 : // Append the calibration buffer to the output calibration table
976 0 : buffer.append(calTable);
977 :
978 : // Only update CAL_DESC if a new row required
979 0 : if (thisDesc==nDesc) {
980 0 : CalDescRecord cdr;
981 0 : cdr.defineSpwId(Vector<Int>(1,solSpwId));
982 0 : cdr.defineMSName(Path(msName()).baseName());
983 0 : calTable.putRowDesc(thisDesc,cdr);
984 : }
985 :
986 0 : return;
987 : }
988 :
989 : //----------------------------------------------------------------------------
990 :
991 0 : Double BJonesPoly::getChebVal (const Vector<Double>& coeff,
992 : const Double& xinit, const Double& xfinal,
993 : const Double& x)
994 : {
995 : // Compute a Chebyshev polynomial value using the CLIC library
996 : // Input:
997 : // coeff const Vector<Double>& Chebyshev coefficients
998 : // xinit const Double& Domain start
999 : // xfinal const Double& Domain end
1000 : // x const Double& x-ordinate
1001 : // Output:
1002 : // getChebVal Double Chebyshev polynomial value
1003 : //
1004 : // Re-scale x-ordinate
1005 0 : Double xcap=0;
1006 0 : xcap=((x-xinit)-(xfinal-x))/(xfinal-xinit);
1007 :
1008 : // Compute polynomial
1009 0 : Int deg=coeff.shape().asVector()(0);
1010 0 : Vector<Double> val(deg);
1011 : Bool check;
1012 : Int checkval;
1013 0 : cheb(°,
1014 : &xcap,
1015 : val.getStorage(check),
1016 : &checkval);
1017 :
1018 0 : Double soly=0.0;
1019 0 : for (Int mm=0; mm< deg; mm++){
1020 0 : soly+= coeff[mm]*val[mm];
1021 : }
1022 :
1023 0 : return soly;
1024 : }
1025 :
1026 : //----------------------------------------------------------------------------
1027 :
1028 0 : Bool BJonesPoly::maskedChannel (const Int& chan, const Int& nChan)
1029 : {
1030 : // Check if a given channel is masked or not
1031 : // Input:
1032 : // chan const Int& Channel number
1033 : // nChan const Int& No. of channels in spectrum
1034 : // Output:
1035 : // maskedChannel Bool Returns true if channel lies
1036 : // in edge or center mask
1037 : //
1038 : // Initialization
1039 0 : Bool masked = false;
1040 0 : Int loChan = nChan / 2 - maskcenterHalf_p;
1041 0 : Int hiChan = loChan + maskcenter_p;
1042 :
1043 : // Check mask at center of spectrum
1044 0 : if ((chan >= loChan) && (chan < hiChan)) {
1045 0 : masked = true;
1046 : };
1047 :
1048 : // Check mask at edge of spectrum
1049 0 : if ((chan < (nChan*maskedgeFrac_p)) || (chan > nChan*(1-maskedgeFrac_p))) {
1050 0 : masked = true;
1051 : };
1052 :
1053 0 : return masked;
1054 : };
1055 :
1056 : //----------------------------------------------------------------------------
1057 :
1058 0 : void BJonesPoly::loadMemCalTable (String applyTable,String /*field*/)
1059 : {
1060 :
1061 : // Load and cache the polynomial bandpass corrections
1062 : // Input:
1063 : // applyTable const String& Cal. table to be applied
1064 : // Output to protected data:
1065 : // Antenna and interferometer bandpass correction caches
1066 : // in the BJones parent class.
1067 : //
1068 :
1069 : // NB: For the moment we will read the Chebychev params and calculate
1070 : // a sampled (complex) bandpass on a per-channel basis. Downstream,
1071 : // BPOLY will thus behave as if it were an ordinary B (except for the
1072 : // fact that currently it CANNOT be time-dependent).
1073 :
1074 : // Generate a NCT in memory to hold the BPOLY as a B
1075 0 : ct_=new NewCalTable("BpolyAsB.temp",VisCalEnum::COMPLEX,"B Jones",msName(),false);
1076 :
1077 : // Open the BJonesPoly calibration table
1078 0 : BJonesPolyTable calTable(applyTable, Table::Update);
1079 :
1080 : // Ensure sort on TIME, so CalSet is filled in order
1081 0 : Block <String> sortCol(3);
1082 0 : sortCol[0]="CAL_DESC_ID";
1083 0 : sortCol[1]="TIME";
1084 0 : sortCol[2]="ANTENNA1";
1085 0 : calTable.sort2(sortCol);
1086 :
1087 : // Attach a calibration table columns accessor
1088 0 : BJonesPolyMCol col(calTable);
1089 :
1090 : // Fill the bandpass correction cache
1091 0 : Int nrows = calTable.nRowMain();
1092 0 : Int nDesc = calTable.nRowDesc();
1093 :
1094 : // A matrix to show which spws to be calibrated by each caldesc
1095 0 : Vector<Int> spwmap(nDesc,-1);
1096 0 : for (Int idesc=0;idesc<nDesc;++idesc) {
1097 :
1098 : // This cal desc
1099 0 : CalDescRecord calDescRec(calTable.getRowDesc(idesc));
1100 :
1101 : // Get this spw ID
1102 0 : Vector<Int> spwId;
1103 0 : calDescRec.getSpwId(spwId);
1104 0 : Int nSpw; spwId.shape(nSpw);
1105 0 : if (nSpw > 1) {}; // ERROR!!! Should only be one spw per cal desc!
1106 0 : spwmap(idesc)=spwId(0);
1107 :
1108 0 : currSpw()=spwId(0);
1109 0 : Vector<Double> freq = freqAxis(currSpw());
1110 :
1111 0 : startChan()=0;
1112 0 : nChanPar()=freq.nelements(); // currSpw
1113 :
1114 : // Set SPW subtable freqs
1115 : //IPosition blc(1,startChan()), trc(1,startChan()+nChanPar()-1);
1116 0 : ct_->setSpwFreqs(currSpw(),freq); // (blc,trc));
1117 :
1118 : }
1119 :
1120 : // We will fill in a _sampled_ bandpass from the Cheby coeffs
1121 0 : for (Int ispw=0; ispw<nSpw(); ispw++) {
1122 0 : currSpw()=ispw;
1123 : // currCPar().resize(nPar(),nChanPar(),nAnt());
1124 : // currCPar()=Complex(1.0);
1125 : // currParOK().resize(nPar(),nChanPar(),nAnt());
1126 : // currParOK()=false;
1127 0 : invalidateP();
1128 0 : invalidateCalMat();
1129 : }
1130 :
1131 : // Inflate the solve arrays, so we can fill them
1132 0 : initSolvePar();
1133 :
1134 0 : for (Int row=0; row < nrows; row++) {
1135 :
1136 0 : Int idesc = col.calDescId().asInt(row);
1137 :
1138 0 : Double thisTime=col.time().asdouble(row);
1139 0 : Int thisSpw=spwmap(idesc);
1140 :
1141 : // Antenna id.
1142 0 : Int antennaId = col.antenna1().asInt(row);
1143 :
1144 0 : currSpw()=thisSpw;
1145 0 : refTime()=thisTime;
1146 :
1147 : // Frequency group name
1148 0 : String freqGrpName = col.freqGrpName().asString(row);
1149 :
1150 : // Extract the polynomial scale factor
1151 0 : Complex factor = col.scaleFactor().asComplex(row);
1152 :
1153 : // Extract the valid domain for the polynomial
1154 0 : Vector<Double> freqDomain(2);
1155 0 : col.validDomain().get(row, freqDomain);
1156 :
1157 : // Extract the polynomial coefficients in amplitude and phase
1158 0 : Int nAmp = col.nPolyAmp().asInt(row);
1159 0 : Int nPhase = col.nPolyPhase().asInt(row);
1160 0 : Matrix<Double> ampCoeff(nAmp,2);
1161 0 : Matrix<Double> phaseCoeff(nPhase,2);
1162 0 : Array<Double> ampCoeffArray, phaseCoeffArray;
1163 0 : col.polyCoeffAmp().get(row, ampCoeffArray);
1164 0 : col.polyCoeffPhase().get(row, phaseCoeffArray);
1165 :
1166 0 : IPosition ampPos = ampCoeffArray.shape();
1167 0 : ampPos = 0;
1168 0 : for (Int k=0; k < 2*nAmp; k++) {
1169 0 : ampPos.setLast(IPosition(1,k));
1170 0 : ampCoeff(k%nAmp,k/nAmp) = ampCoeffArray(ampPos);
1171 : };
1172 :
1173 0 : IPosition phasePos = phaseCoeffArray.shape();
1174 0 : phasePos = 0;
1175 0 : for (Int k=0; k < 2*nPhase; k++) {
1176 0 : phasePos.setLast(IPosition(1,k));
1177 0 : phaseCoeff(k%nPhase,k/nPhase) = phaseCoeffArray(phasePos);
1178 : };
1179 :
1180 : // Loop over all spectral window id.'s in this frequency group
1181 : // Vector<Int> spwIds = spwIdsInGroup(freqGrpName);
1182 :
1183 : // This cal desc
1184 0 : CalDescRecord calDescRec(calTable.getRowDesc(idesc));
1185 :
1186 : // Get this spw ID
1187 0 : Vector<Int> spwIds;
1188 0 : calDescRec.getSpwId(spwIds);
1189 :
1190 :
1191 0 : Vector<Double> freq = freqAxis(currSpw());
1192 0 : Int nChan=freq.nelements();
1193 :
1194 0 : Double x1 = freqDomain(0);
1195 0 : Double x2 = freqDomain(1);
1196 :
1197 0 : for (Int ipol=0;ipol<2;ipol++) {
1198 :
1199 0 : Vector<Double> ac(ampCoeff.column(ipol));
1200 0 : Vector<Double> pc(phaseCoeff.column(ipol));
1201 :
1202 : // Only do non-triv calc if coeffs are non-triv
1203 0 : if (anyNE(ac,Double(0.0)) ||
1204 0 : anyNE(pc,Double(0.0)) ) {
1205 :
1206 : // Loop over frequency channel
1207 0 : for (Int chan=0; chan < nChan; ++chan) {
1208 :
1209 0 : Double ampval(1.0),phaseval(0.0);
1210 : // only if in domain, calculate Cheby
1211 0 : Double thisfreq(freq(chan));
1212 0 : if (thisfreq >=x1 && thisfreq <= x2) {
1213 0 : ampval = getChebVal(ac, x1, x2, thisfreq);
1214 0 : phaseval = getChebVal(pc, x1, x2, thisfreq);
1215 0 : solveAllCPar()(ipol,chan,antennaId)= factor *
1216 0 : Complex(exp(ampval)) * Complex(cos(phaseval),sin(phaseval));
1217 0 : solveAllParOK()(ipol,chan,antennaId)= true;
1218 : }
1219 : else {
1220 : // Unflagged unit calibration for now
1221 0 : solveAllCPar()(ipol,chan,antennaId)= Complex(1.0);
1222 0 : solveAllParOK()(ipol,chan,antennaId)= true;
1223 : }
1224 : }
1225 : }
1226 :
1227 : } // ipol
1228 :
1229 : // Every nAnt rows, store the result
1230 0 : if ((row+1)%nAnt()==0) {
1231 0 : ct_->fillAntBasedMainRows(nAnt(),refTime(),0.0,-1,currSpw(),
1232 0 : -1,Vector<Int>(),-1,
1233 0 : solveAllCPar(),!solveAllParOK(),
1234 0 : solveAllParErr(),solveAllParSNR());
1235 :
1236 : // reset arrays
1237 0 : solveAllCPar().set(Complex(1.0));
1238 0 : solveAllParOK().set(false);
1239 0 : solveAllParErr().set(0.0);
1240 0 : solveAllParSNR().set(1.0);
1241 :
1242 : }
1243 :
1244 : } // rows
1245 :
1246 : /*
1247 : String cTN;
1248 : cTN=calTableName();
1249 : calTableName()=calTableName()+".BpolyAsB";
1250 : storeNCT();
1251 : calTableName()=cTN;
1252 : */
1253 :
1254 0 : return;
1255 : }
1256 :
1257 : //----------------------------------------------------------------------------
1258 :
1259 0 : Double BJonesPoly::meanFrequency (const Vector<Int>& spwid)
1260 : {
1261 : // Compute the bandwidth-weighted average frequency of a set of spw id.'s
1262 : // Input:
1263 : // spwid const Vector<Int>& Spectral window id.'s
1264 : // Input from private data:
1265 : // vs_p VisSet* Current visibility set
1266 : // Output:
1267 : // meanFrequency Double Mean frequency (as Double)
1268 : //
1269 :
1270 0 : if (!vs_p) throw(AipsError("Error in BJonesPoly::meanFrequency"));
1271 :
1272 0 : const MSColumns& mscol(vs_p->iter().msColumns());
1273 0 : const MSSpWindowColumns& spwcol(mscol.spectralWindow());
1274 0 : const ArrayColumn<Double>& frequencies(spwcol.chanFreq());
1275 0 : const ScalarColumn<Double>& totalbw(spwcol.totalBandwidth());
1276 :
1277 0 : Int numspw=spwid.shape().asVector()(0);
1278 :
1279 0 : Double meanFreq=0.0;
1280 0 : Double sumbw=0.0;
1281 0 : for (Int k=0; k<numspw; k++){
1282 0 : meanFreq = meanFreq +
1283 0 : sum(frequencies(spwid(k)))*totalbw(spwid(k)) /
1284 0 : (frequencies(spwid(k)).nelements());
1285 0 : sumbw = sumbw + totalbw(spwid(k));
1286 : }
1287 : // Compute the mean frequency
1288 0 : meanFreq=meanFreq/sumbw;
1289 0 : return meanFreq;
1290 : }
1291 :
1292 : //----------------------------------------------------------------------------
1293 :
1294 0 : String BJonesPoly::freqGrpName (const Int& spwId)
1295 : {
1296 : // Return the frequency group name for a given spectral window id.
1297 : // Input:
1298 : // spwId const Int& Spectral window id.
1299 : // Input from private data:
1300 : // vs_p VisSet* Current visibility set
1301 : // Output:
1302 : // freqGrpName String Frequency group name
1303 : //
1304 :
1305 0 : if (!vs_p) throw(AipsError("Error in BJonesPoly::freqGrpName"));
1306 :
1307 0 : const MSColumns& mscol(vs_p->iter().msColumns());
1308 0 : const MSSpWindowColumns& spwCol(mscol.spectralWindow());
1309 :
1310 0 : return spwCol.freqGroupName().asString(spwId);
1311 : }
1312 :
1313 : //----------------------------------------------------------------------------
1314 :
1315 0 : Vector<Int> BJonesPoly::spwIdsInGroup (const String& freqGrpName)
1316 : {
1317 : // Return the spw. id.'s in a freq. group of a given name
1318 : // Input:
1319 : // freqGrpName const String& Frequency group name
1320 : // Input from private data:
1321 : // vs_p VisSet* Current visibility set
1322 : // Output:
1323 : // spwIdsInGroup Vector<Int> Spw. id.'s in freq. group
1324 : //
1325 : // Open a SPECTRAL_WINDOW sub-table index
1326 0 : if (!vs_p) throw(AipsError("Error in BJones::spwIdsInGroup"));
1327 0 : MeasurementSet ms(vs_p->msName().before("/SELECTED"));
1328 0 : MSSpWindowIndex spwIndex(ms.spectralWindow());
1329 0 : return spwIndex.matchFreqGrpName(freqGrpName);
1330 : }
1331 :
1332 : //----------------------------------------------------------------------------
1333 :
1334 0 : Vector<Double> BJonesPoly::freqAxis (const Int& spwId)
1335 : {
1336 : // Return the frequency axis values for a given spectral window id.
1337 : // Input:
1338 : // spwId const Int& Spectral window id.
1339 : // Input from private data:
1340 : // vs_p VisSet* Current visibility set
1341 : // Output:
1342 : // freqAxis Vector<Double> Frequency axis values
1343 : //
1344 :
1345 0 : Vector<Double> freqVal;
1346 0 : if (vs_p) {
1347 0 : const MSColumns& mscol(vs_p->iter().msColumns());
1348 0 : const MSSpWindowColumns& spwCol(mscol.spectralWindow());
1349 :
1350 0 : spwCol.chanFreq().get(spwId, freqVal);
1351 : }
1352 : else {
1353 : // Try msmc...
1354 0 : const MSColumns& mscol(*(msmc().msmd().getMS()));
1355 0 : const MSSpWindowColumns& spwCol(mscol.spectralWindow());
1356 0 : spwCol.chanFreq().get(spwId, freqVal);
1357 : }
1358 :
1359 0 : return freqVal;
1360 : }
1361 :
1362 :
1363 : //----------------------------------------------------------------------------
1364 0 : void BJonesPoly::plotsolve2(const Vector<Double>& x,
1365 : const Matrix<Double>& ampdata,
1366 : const Matrix<Double>& phadata,
1367 : const Matrix<Double>& wtdata,
1368 : const Vector<Int>& ant1idx, const Vector<Int>& ant2idx,
1369 : const Vector<Double>& amperr, Matrix<Double>& ampcoeff,
1370 : const Vector<Double>& phaerr, Matrix<Double>& phacoeff)
1371 :
1372 : {
1373 : // Function to plot and compare Bandpass data and solution
1374 : // Input:
1375 : // x - vector of frequecies
1376 : // ampdata - matrix of data - shape is yall(freqindex, baselineindex)
1377 : // phadata - matrix of data - shape is yall(freqindex, baselineindex)
1378 : // wtdata - matrix of weight; same shape as yall
1379 : // ant1idx - indices for antenna1
1380 : // ant2idx - indices for antenna2
1381 : // amperr - baseline-based amplitude fit errors
1382 : // ampcoeff - antenna-based amplitude poly coefficients
1383 : // phaerr - baseline-based phase fit errors
1384 : // phacoeff - antenna-based phase poly coefficients
1385 :
1386 0 : LogIO os (LogOrigin("BJonesPoly", "plotsolve2()", WHERE));
1387 :
1388 0 : String device;
1389 0 : device=calTableName() + ".ps/cps";
1390 :
1391 : os << LogIO::NORMAL
1392 : << "Generating plot file:" << device
1393 0 : << LogIO::POST;
1394 :
1395 0 : PGPlotter pg(device);
1396 0 : pg.subp(4,3);
1397 :
1398 0 : Int ndatapts= max(x.shape().asVector());
1399 0 : Int nmodelpts=4*(ndatapts-2)-1;
1400 0 : Int num_valid_points=0;
1401 :
1402 0 : Int numbasl=ampdata.shape()(1);
1403 :
1404 0 : Vector<Float> x1(ndatapts);
1405 0 : Vector<Float> amp1(ndatapts);
1406 0 : Vector<Float> pha1(ndatapts);
1407 0 : Vector<Float> amperr1(ndatapts);
1408 0 : Vector<Float> phaerr1(ndatapts);
1409 0 : Vector<Double> weight;
1410 :
1411 0 : Vector<Double> xval(nmodelpts);
1412 0 : Vector<Float> xfloatval(nmodelpts);
1413 :
1414 0 : Vector<Float> amp2a(nmodelpts);
1415 0 : Vector<Float> amp2b(nmodelpts);
1416 0 : Vector<Float> amp2(nmodelpts);
1417 0 : Vector<Float> pha2a(nmodelpts);
1418 0 : Vector<Float> pha2b(nmodelpts);
1419 0 : Vector<Float> pha2(nmodelpts);
1420 :
1421 : Double minfreq,maxfreq;
1422 0 : minfreq=min(x);
1423 0 : maxfreq=max(x);
1424 :
1425 : Double dmodfreq;
1426 0 : dmodfreq = (maxfreq-minfreq)/Double(ndatapts-1)/4;
1427 :
1428 0 : Vector<Double> ant1ampcoeff, ant2ampcoeff;
1429 0 : Vector<Double> ant1phacoeff, ant2phacoeff;
1430 :
1431 0 : Vector<Float> meanamp(numbasl,0.0);
1432 0 : Vector<Float> meanampmod(numbasl,0.0);
1433 0 : Vector<Float> meanpha(numbasl,0.0);
1434 0 : Vector<Float> meanphamod(numbasl,0.0);
1435 0 : for (Int ibl=0; ibl<numbasl; ibl++) {
1436 :
1437 : // refresh all plotting variables
1438 0 : x1.resize(ndatapts);
1439 0 : amp1.resize(ndatapts);
1440 0 : pha1.resize(ndatapts);
1441 0 : amperr1.resize(ndatapts);
1442 0 : phaerr1.resize(ndatapts);
1443 0 : weight.resize();
1444 0 : weight=wtdata.column(ibl);
1445 :
1446 0 : num_valid_points=0;
1447 :
1448 :
1449 0 : for(Int k=0; k < ndatapts; k++){
1450 0 : if (weight(k)>0){
1451 0 : x1(num_valid_points)=(x(k)-minfreq)/1e3;
1452 :
1453 0 : amp1(num_valid_points)=exp(ampdata.column(ibl)(k));
1454 0 : amperr1(num_valid_points)=Float(amperr(ibl))*amp1(num_valid_points);
1455 :
1456 0 : pha1(num_valid_points)=remainder(phadata.column(ibl)(k), 2*C::pi)/C::pi*180.0;
1457 0 : phaerr1(num_valid_points)=Float(phaerr(ibl))/C::pi*180.0;
1458 :
1459 0 : num_valid_points+=1;
1460 : }
1461 : }
1462 :
1463 0 : if (num_valid_points > 0){
1464 :
1465 : // resize data arrays according to num_valid_points
1466 0 : x1.resize(num_valid_points, true);
1467 0 : amp1.resize(num_valid_points, true);
1468 0 : amperr1.resize(num_valid_points, true);
1469 0 : pha1.resize(num_valid_points, true);
1470 0 : phaerr1.resize(num_valid_points, true);
1471 :
1472 0 : ant1ampcoeff=ampcoeff.row(ant1idx(ibl));
1473 0 : ant2ampcoeff=ampcoeff.row(ant2idx(ibl));
1474 0 : ant1phacoeff=phacoeff.row(ant1idx(ibl));
1475 0 : ant2phacoeff=phacoeff.row(ant2idx(ibl));
1476 :
1477 0 : for(Int k=0; k < nmodelpts; k++){
1478 :
1479 0 : xval[k]= Double(k-1)*dmodfreq + x[1];
1480 :
1481 : // Float version, offset & scaled, for plotting
1482 0 : xfloatval(k)=Float((xval(k)-minfreq)/1.0e3);
1483 :
1484 0 : pha2(k)=0;
1485 0 : pha2a(k)=getChebVal(ant1phacoeff,x(0),x(ndatapts-1), xval[k]);
1486 0 : pha2b(k)=getChebVal(ant2phacoeff,x(0),x(ndatapts-1), xval[k]);
1487 0 : pha2(k)=pha2b(k)-pha2a(k);
1488 :
1489 0 : pha2a(k)=remainder(pha2a(k),2*C::pi)*180.0/C::pi;
1490 0 : pha2b(k)=remainder(pha2b(k),2*C::pi)*180.0/C::pi;
1491 0 : pha2(k)=remainder(pha2(k),2*C::pi)*180.0/C::pi;
1492 :
1493 :
1494 0 : amp2a(k)=getChebVal(ant1ampcoeff,x(0),x(ndatapts-1), xval(k));
1495 0 : amp2b(k)=getChebVal(ant2ampcoeff,x(0),x(ndatapts-1), xval(k));
1496 0 : amp2(k)=amp2b(k) + amp2a(k);
1497 :
1498 0 : amp2a(k)=exp(amp2a(k));
1499 0 : amp2b(k)=exp(amp2b(k));
1500 0 : amp2(k)=exp(amp2(k));
1501 :
1502 : }
1503 :
1504 0 : meanamp(ibl) = mean(amp1);
1505 0 : meanampmod(ibl) = mean(amp2);
1506 0 : meanpha(ibl) = mean(pha1);
1507 0 : meanphamod(ibl) = mean(pha2);
1508 :
1509 :
1510 : /*
1511 : cout << ant1idx(ibl) << "-" << ant2idx(ibl) << " Means: "
1512 : << mean(amp1) << " "
1513 : << mean(amp2) << " "
1514 : << mean(pha1) << " "
1515 : << mean(pha2) << endl;
1516 : cout << "means = "
1517 : << mean(meanamp) << " "
1518 : << mean(meanampmod) << " "
1519 : << mean(meanpha) << " "
1520 : << mean(meanphamod) << " "
1521 : << endl;
1522 : */
1523 :
1524 : //cout << " Mean of sol " << mean(soly1) << endl;
1525 :
1526 0 : Float min_x= 0.0;
1527 0 : Float max_x= Float((maxfreq-minfreq)/1.0e3);
1528 :
1529 : Float min_amp, max_amp, min_pha, max_pha;
1530 0 : min_amp=min( (min(amp1)-3.0*max(amperr1)), min(amp2) );
1531 0 : max_amp=max( (max(amp1)+3.0*max(amperr1)), max(amp2) );
1532 0 : min_pha=min( (min(pha1)-3.0*max(phaerr1)), min(pha2) );
1533 0 : max_pha=max( (max(pha1)+3.0*max(phaerr1)), max(pha2) );
1534 :
1535 0 : min_amp=min(amp1)-1.0*max(amperr1);
1536 0 : max_amp=max(amp1)+1.0*max(amperr1);
1537 0 : min_pha=min(pha1)-1.0*max(phaerr1);
1538 0 : max_pha=max(pha1)+1.0*max(phaerr1);
1539 :
1540 0 : Float damp=0.1*(max_amp-min_amp);
1541 0 : Float dpha=0.1*(max_pha-min_pha);
1542 0 : min_amp-=damp;
1543 0 : max_amp+=damp;
1544 0 : min_pha-=dpha;
1545 0 : max_pha+=dpha;
1546 :
1547 :
1548 0 : pg.page();
1549 :
1550 : // arrange labels
1551 0 : ostringstream titlelab, xlab, ampdeg, phadeg, pharms, amprms;
1552 0 : xlab.precision(12);
1553 0 : xlab << "Frequency in kHz (-" << minfreq/1.0e9 << " GHz)";
1554 0 : titlelab << "B polynomial for baseline "
1555 0 : << ant1idx(ibl) << " & " << ant2idx(ibl);
1556 0 : ampdeg << "degree = " << degamp_p;
1557 0 : phadeg << "degree = " << degphase_p;
1558 0 : pharms << "rms = " << phaerr1(0);
1559 0 : amprms << "rms = " << amperr1(0);
1560 :
1561 :
1562 : // plot phase in upper half
1563 0 : pg.sci(1); pg.sch(1.5);
1564 : // pg.svp(0.13,0.87,0.5,0.85);
1565 0 : pg.svp(0.10,0.90,0.525,0.90);
1566 0 : pg.swin(min_x,max_x,min_pha,max_pha);
1567 0 : pg.box("BCST",0.0,0,"BCNST",0.0,0);
1568 :
1569 0 : pg.sch(0.75);
1570 0 : pg.sci(2);
1571 0 : pg.errb(6, x1, pha1, phaerr1, 2.0);
1572 0 : pg.sci(3);
1573 0 : pg.pt(x1, pha1, 17);
1574 0 : pg.sci(4); pg.sls(1);
1575 0 : pg.line(xfloatval, pha2);
1576 : /*
1577 : pg.sci(4); pg.sls(4); // plot each antenna soln
1578 : pg.line(xfloatval, pha2a);
1579 : pg.line(xfloatval, pha2b);
1580 : pg.sls(1);
1581 : */
1582 0 : pg.sci(1); pg.sch(1.5);
1583 0 : pg.lab(" ", "Phase (deg)", " ");
1584 0 : pg.mtxt("T",0.75,0.5,0.5,titlelab);
1585 0 : pg.sch(1.2);
1586 0 : pg.mtxt("T",-1.5,0.95,1.0,phadeg);
1587 0 : pg.mtxt("B",-1.5,0.95,1.0,pharms);
1588 :
1589 :
1590 : // plot amp in upper half
1591 0 : pg.sci(1); pg.sch(1.5);
1592 :
1593 0 : pg.svp(0.10,0.90,0.15,0.525);
1594 0 : pg.swin(min_x,max_x,min_amp,max_amp);
1595 0 : pg.box("BCNST",0.0,0,"BCNST",0.0,0);
1596 0 : pg.sch(0.75);
1597 0 : pg.sci(2);
1598 0 : pg.errb(6, x1, amp1, amperr1, 2.0);
1599 0 : pg.sci(3);
1600 0 : pg.pt(x1, amp1, 17);
1601 0 : pg.sci(4); pg.sls(1);
1602 0 : pg.line(xfloatval, amp2);
1603 : /*
1604 : pg.sci(4); pg.sls(4); // plot each antenna soln
1605 : pg.line(xfloatval, amp2a);
1606 : pg.line(xfloatval, amp2b);
1607 : pg.sls(1);
1608 : */
1609 0 : pg.sci(1) ; pg.sch(1.5);
1610 0 : pg.lab(xlab, "Amp", " ");
1611 0 : pg.sch(1.2);
1612 0 : pg.mtxt("T",-1.5,0.95,1.0,ampdeg);
1613 0 : pg.mtxt("B",-1.5,0.95,1.0,amprms);
1614 :
1615 :
1616 : } else {
1617 : // shouldn't get here because this wouldn't be one of the good baselines
1618 0 : pg.sci(1);
1619 0 : pg.env(0.0,1.0,0.0,1.0,0,-2);
1620 0 : pg.sch(2.0);
1621 0 : ostringstream oos;
1622 0 : oos << "No data for baseline " << ant1idx(ibl) << " & " << ant2idx(ibl);
1623 0 : pg.ptxt(0.5,0.5,0.0,0.5,oos);
1624 0 : pg.sch(1.0);
1625 : }
1626 : }
1627 0 : }
1628 :
1629 : //---------------------------------------------------------------------------
1630 :
1631 : } //# NAMESPACE CASA - END
1632 :
|
