Line data Source code
1 : //# KJones.cc: Implementation of KJones
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003,2011
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 :
27 : #include <synthesis/MeasurementComponents/KJones.h>
28 :
29 : #include <msvis/MSVis/VisBuffer.h>
30 : #include <msvis/MSVis/VisBuffer2.h>
31 : #include <msvis/MSVis/VisBuffAccumulator.h>
32 : #include <casacore/ms/MeasurementSets/MSColumns.h>
33 : #include <casacore/ms/MSOper/MSMetaData.h>
34 : #include <synthesis/MeasurementEquations/VisEquation.h> // *
35 : #include <synthesis/MeasurementComponents/SolveDataBuffer.h>
36 : #include <synthesis/MeasurementComponents/MSMetaInfoForCal.h>
37 : #include <casacore/lattices/Lattices/ArrayLattice.h>
38 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
39 : #include <casacore/scimath/Mathematics/FFTServer.h>
40 : #include <casacore/casa/Quanta/QuantumHolder.h>
41 :
42 : #include <casacore/casa/Arrays/ArrayMath.h>
43 : #include <casacore/casa/Arrays/MatrixMath.h>
44 : #include <casacore/casa/BasicSL/String.h>
45 : #include <casacore/casa/Utilities/Assert.h>
46 : #include <casacore/casa/Exceptions/Error.h>
47 : #include <casacore/casa/System/Aipsrc.h>
48 :
49 : #include <sstream>
50 :
51 : #include <casacore/measures/Measures/MCBaseline.h>
52 : #include <casacore/measures/Measures/MDirection.h>
53 : #include <casacore/measures/Measures/MEpoch.h>
54 : #include <casacore/measures/Measures/MeasTable.h>
55 :
56 : #include <casacore/casa/Logging/LogMessage.h>
57 : #include <casacore/casa/Logging/LogSink.h>
58 :
59 :
60 : using namespace casa::vi;
61 : using namespace casacore;
62 :
63 : namespace casa { //# NAMESPACE CASA - BEGIN
64 :
65 :
66 : // **********************************************************
67 : // DelayFFT Implementations
68 : //
69 :
70 : // Construct from freq info and a data-like Cube<Complex>
71 0 : DelayFFT::DelayFFT(Double f0, Double df, Double padBW,
72 0 : Cube<Complex> V,Cube<Float> wt) :
73 : f0_(f0),
74 : df_(df),
75 : padBW_(padBW),
76 0 : nCorr_(V.shape()(0)),
77 0 : nPadChan_(Int(0.5+padBW/abs(df))),
78 0 : nElem_(V.shape()(2)),
79 : refant_(-1), // ok?
80 : Vpad_(),
81 : delay_(),
82 0 : flag_()
83 : {
84 :
85 : // cout << "pad ch=" << padBW/df_ << " " << nPadChan_*df << endl;
86 :
87 0 : IPosition ip1(3,nCorr_,nPadChan_,nElem_);
88 0 : Vpad_.resize(ip1);
89 0 : Vpad_.set(0.0);
90 :
91 0 : Slicer sl1(Slice(),Slice(0,V.shape()(1),1),Slice());
92 0 : if (V.shape()==wt.shape())
93 0 : Vpad_(sl1)=(V*wt);
94 : else
95 0 : Vpad_(sl1)=V;
96 :
97 : //this->state();
98 0 : }
99 :
100 :
101 : // Construct from freq info and shape, w/ initialization
102 12 : DelayFFT::DelayFFT(Double f0, Double df, Double padBW,
103 12 : Int nCorr, Int nElem, Int refant, Complex v0, Float wt) :
104 : f0_(f0),
105 : df_(df),
106 : padBW_(padBW),
107 : nCorr_(nCorr),
108 24 : nPadChan_(Int(0.5+padBW/abs(df))),
109 : nElem_(nElem),
110 : refant_(refant),
111 : Vpad_(),
112 : delay_(),
113 12 : flag_()
114 : {
115 :
116 : // cout << "ctor0: " << f0 << " " << df << " " << padBW << " " << nPadChan_ << endl;
117 : // cout << "ctor0: " << f0_ << " " << df_ << " " << padBW_ << " " << nPadChan_ << endl;
118 : // cout << "pad ch=" << padBW/df_ << " " << nPadChan_*df << endl;
119 :
120 12 : Vpad_.resize(nCorr_,nPadChan_,nElem_);
121 12 : Vpad_.set(v0);
122 12 : Vpad_*=wt; // scale by supplied per-channel weight
123 :
124 : // this->state();
125 :
126 12 : }
127 :
128 :
129 0 : DelayFFT::DelayFFT(const VisBuffer& vb,Double padBW,Int refant) :
130 0 : f0_(vb.frequency()(0)/1.e9), // GHz
131 0 : df_(vb.frequency()(1)/1.e9-f0_),
132 : padBW_(padBW),
133 : nCorr_(0), // set in body
134 0 : nPadChan_(Int(0.5+padBW/abs(df_))),
135 0 : nElem_(vb.numberAnt()), // antenna-based
136 : refant_(refant),
137 : Vpad_(),
138 : delay_(),
139 0 : flag_()
140 : {
141 :
142 : // cout << "DelayFFT(vb)..." << endl;
143 : // cout << "ctor1: " << f0_ << " " << df_ << " " << padBW_ << " " << nPadChan_ << endl;
144 :
145 : // VB facts
146 0 : Int nCorr=vb.nCorr();
147 0 : Int nChan=vb.nChannel();
148 0 : Int nRow=vb.nRow();
149 :
150 : // Discern effective shapes
151 0 : nCorr_=(nCorr>1 ? 2 : 1); // number of p-hands
152 0 : Int sC=(nCorr>2 ? 3 : 1); // step for p-hands
153 :
154 : // Shape the data
155 0 : IPosition ip1(3,nCorr_,nPadChan_,nElem_);
156 0 : Vpad_.resize(ip1);
157 0 : Vpad_.set(Complex(0.0));
158 :
159 : // this->state();
160 :
161 : // Fill the relevant data
162 0 : Int iant=0;
163 0 : for (Int irow=0;irow<nRow;++irow) {
164 0 : Int a1(vb.antenna1()(irow)), a2(vb.antenna2()(irow));
165 0 : if (!vb.flagRow()(irow) && a1!=a2) {
166 :
167 0 : if (a1==refant)
168 0 : iant=a2;
169 0 : else if (a2==refant)
170 0 : iant=a1;
171 : else
172 0 : continue; // an irrelevant baseline
173 :
174 0 : Slicer sl0(Slice(0,nCorr_,sC),Slice(),Slice(irow,1,1)); // this visCube slice
175 0 : Slicer sl1(Slice(),Slice(0,nChan,1),Slice(iant,1,1)); // this Vpad_ slice
176 :
177 0 : Cube<Complex> vC(vb.visCube()(sl0));
178 :
179 : // Divide by non-zero amps
180 0 : Cube<Float> vCa(amplitude(vC));
181 0 : vCa(vCa<FLT_EPSILON)=1.0;
182 0 : vC/=vCa;
183 :
184 : // Zero flagged channels
185 0 : Slicer fsl0(Slice(),Slice(irow,1,1));
186 0 : Array<Bool> fl(vb.flag()(fsl0).nonDegenerate(IPosition(1,0)));
187 0 : for (Int icor=0;icor<nCorr_;++icor)
188 0 : vC(Slice(icor,1,1),Slice(),Slice()).nonDegenerate(IPosition(1,1))(fl)=Complex(0.0);
189 :
190 : // TBD: apply weights
191 : // Matrix<Float> wt(vb.weightMat()(Slice,Slice(irow,1,1)));
192 :
193 : // Acquire this baseline for solving
194 0 : Vpad_(sl1)=vC;
195 : }
196 : }
197 : //cout << "...end DelayFFT(vb)" << endl;
198 :
199 0 : }
200 24 : DelayFFT::DelayFFT(SolveDataBuffer& sdb,Double padBW,Int refant,Int nElem) :
201 24 : f0_(sdb.freqs()(0)/1e9), // GHz
202 24 : df_(sdb.freqs()(1)/1e9-f0_),
203 : padBW_(padBW),
204 : nCorr_(0), // set in body
205 48 : nPadChan_(Int(0.5+padBW_/abs(df_))),
206 : nElem_(nElem), // antenna-based
207 : refant_(refant),
208 : Vpad_(),
209 : delay_(),
210 48 : flag_()
211 : {
212 :
213 : // VB facts
214 24 : Int nCorr=sdb.nCorrelations();
215 24 : Int nChan=sdb.nChannels();
216 24 : Int nRow=sdb.nRows();
217 :
218 : // Discern effective shapes
219 24 : nCorr_=(nCorr>1 ? 2 : 1); // number of p-hands
220 24 : Int sC=(nCorr>2 ? 3 : 1); // step for p-hands
221 :
222 : // Shape the data
223 48 : IPosition ip1(3,nCorr_,nPadChan_,nElem_);
224 24 : Vpad_.resize(ip1);
225 24 : Vpad_.set(Complex(0.0));
226 :
227 : // this->state();
228 :
229 : // Fill the relevant data
230 24 : Int iant=0;
231 1104 : for (Int irow=0;irow<nRow;++irow) {
232 1080 : Int a1(sdb.antenna1()(irow)), a2(sdb.antenna2()(irow));
233 1080 : if (!sdb.flagRow()(irow) && a1!=a2) {
234 :
235 1080 : if (a1==refant)
236 216 : iant=a2;
237 864 : else if (a2==refant)
238 0 : iant=a1;
239 : else
240 864 : continue; // an irrelevant baseline
241 :
242 432 : Slicer sl0(Slice(0,nCorr_,sC),Slice(),Slice(irow,1,1)); // this visCube slice
243 432 : Slicer sl1(Slice(),Slice(0,nChan,1),Slice(iant,1,1)); // this Vpad_ slice
244 :
245 432 : Cube<Complex> vC(sdb.visCubeCorrected()(sl0));
246 :
247 : // Divide by non-zero amps
248 : //Cube<Float> vCa(amplitude(vC));
249 : //vCa(vCa<FLT_EPSILON)=1.0;
250 : //vC/=vCa;
251 :
252 : // Zero flagged channels
253 432 : Cube<Bool> fl(sdb.flagCube()(sl0));
254 216 : Cube<Float> wt(sdb.weightSpectrum()(sl0));
255 216 : vC(fl)=Complex(0.0);
256 216 : wt(fl)=0.0f;
257 :
258 : // Scale visibilities by weights
259 216 : vC*=wt;
260 :
261 : // Acquire this baseline for solving
262 216 : Vpad_(sl1)=vC;
263 : }
264 : }
265 : //cout << "...end DelayFFT(sdb)" << endl;
266 :
267 24 : }
268 :
269 24 : void DelayFFT::FFT() {
270 :
271 : // We always transform only the chan axis (1)
272 48 : Vector<Bool> ax(3,false);
273 24 : ax(1)=true;
274 :
275 : // Transform
276 : // TBD: can Vpad_ be an ArrayLattice?
277 48 : ArrayLattice<Complex> c(Vpad_);
278 24 : LatticeFFT::cfft0(c,ax,true);
279 :
280 24 : }
281 :
282 0 : void DelayFFT::shift(Double f) {
283 :
284 : // Shift _this_ at f0_, to _that_ at f
285 :
286 0 : Double shift=-(f0_-f)/df_; // samples
287 0 : Vector<Double> ph(nPadChan_); indgen(ph); // indices
288 : // ph-=Double(nPadChan_/2); // centered
289 0 : ph/=Double(nPadChan_); // cycles/sample
290 0 : ph*=shift; // cycles
291 0 : ph*=(C::_2pi); // rad
292 0 : Vector<Double> fsh(nPadChan_*2);
293 0 : fsh(Slice(0,nPadChan_,2))=cos(ph);
294 0 : fsh(Slice(1,nPadChan_,2))=sin(ph);
295 0 : Vector<DComplex> csh(nPadChan_);
296 0 : RealToComplex(csh,fsh);
297 0 : Vector<Complex> sh(nPadChan_);
298 0 : convertArray(sh,csh); // downcovert to Complex
299 :
300 : // Apply to each elem, corr
301 0 : for (Int ielem=0;ielem<nElem_;++ielem)
302 0 : for (Int icorr=0;icorr<nCorr_;++icorr) {
303 0 : Vector<Complex> v(Vpad_.xyPlane(ielem).row(icorr));
304 0 : v*=sh;
305 : }
306 :
307 0 : }
308 :
309 :
310 :
311 0 : void DelayFFT::add(const DelayFFT& other) {
312 :
313 : // cout << "DelayFFT::add(x)..." << endl;
314 :
315 0 : AlwaysAssert( (other.nCorr_==nCorr_), AipsError);
316 0 : AlwaysAssert( (other.nPadChan_<=nPadChan_), AipsError);
317 0 : AlwaysAssert( (other.nElem_==nElem_), AipsError);
318 :
319 0 : Int oNchan=other.nPadChan_;
320 0 : Int lo(0);
321 0 : while (lo < nPadChan_) {
322 0 : Int nch=min(oNchan,nPadChan_-lo); // ensure we don't overrun (impossible anyway?)
323 0 : Slicer sl0(Slice(),Slice(0,nch,1),Slice());
324 0 : Slicer sl1(Slice(),Slice(lo,nch,1),Slice());
325 0 : Cube<Complex> v1(Vpad_(sl1)), v0(other.Vpad_(sl0));
326 0 : v1+=v0;
327 0 : lo+=oNchan;
328 : }
329 :
330 0 : }
331 :
332 :
333 24 : void DelayFFT::addWithDupAndShift(const DelayFFT& other) {
334 :
335 : // cout << "DelayFFT::addWithDupAndShift(x)..." << endl;
336 :
337 : // Adds other to this, with phase gradient for freq-space shift
338 : // Includes duplication of other (if necessary), to account for resolution/range differences
339 : // (as long as they are congruent: other must have channel bandwidth that is integral multiple of this's)
340 :
341 :
342 24 : const Int& oNPadChan(other.nPadChan_);
343 :
344 : // Verify shapes match adequately
345 24 : AlwaysAssert( (other.nCorr_==nCorr_), AipsError); // Same on corr axis
346 24 : AlwaysAssert( (other.nElem_==nElem_), AipsError); // Same on baseline axis
347 24 : AlwaysAssert( (oNPadChan<=nPadChan_), AipsError); // fewer samples, or equal
348 24 : AlwaysAssert( ((nPadChan_%oNPadChan)==0), AipsError); // Ensures congruency
349 :
350 : // All shapes seem ok, so proceed....
351 :
352 : // Work space for other, dupe'd (if necessary) and shifted
353 48 : Cube<Complex> oDupShifted(Vpad_.shape(),Complex(0.0));
354 :
355 : // Duplicate other into workspace
356 24 : Int lo(0);
357 48 : while (lo < nPadChan_) {
358 48 : Slicer sl0(Slice(),Slice(0,oNPadChan,1),Slice());
359 48 : Slicer sl1(Slice(),Slice(lo,oNPadChan,1),Slice());
360 72 : Cube<Complex> v1(oDupShifted(sl1)), v0(other.Vpad_(sl0));
361 24 : v1+=v0;
362 24 : lo+=oNPadChan;
363 : }
364 :
365 : /*
366 : // Add incoming to low and (if nec.) high end of accumulator
367 : Slicer sl0(Slice(),Slice(0,oNPadChan,1),Slice());
368 : Cube<Complex> v1(oDupShifted(sl0)), v0(other.Vpad_(sl0));
369 : v1+=v0;
370 : if (oNPadChan<nPadChan_) {
371 : Slicer sl1(Slice(),Slice(nPadChan_-oNPadChan,oNPadChan,1),Slice());
372 : Cube<Complex> v2(oDupShifted(sl1));
373 : v2+=v0;
374 : }
375 : */
376 :
377 : // Form and apply shift
378 :
379 : // Generate delay-dep phase shift Vector
380 48 : Vector<Double> ph(nPadChan_); indgen(ph); // delay index
381 24 : Int nPC2(nPadChan_/2);
382 48 : Vector<Double> ph2(ph(Slice(nPC2,nPC2,1))); // (2nd half is negative)
383 24 : ph2-=Double(nPadChan_); // unfolded
384 :
385 24 : ph/=Double(oNPadChan); // cycles/sample for INCOMING delay function!
386 :
387 24 : Double shift=-(other.f0_-f0_); // Shift _other_ to this at f0_
388 24 : shift/=other.df_; // Shift in samples (for INCOMING!)
389 24 : ph*=shift; // cycles
390 24 : ph*=(C::_2pi); // rad
391 :
392 48 : Vector<Double> fsh(nPadChan_*2);
393 24 : fsh(Slice(0,nPadChan_,2))=cos(ph);
394 24 : fsh(Slice(1,nPadChan_,2))=sin(ph);
395 48 : Vector<DComplex> csh(nPadChan_);
396 24 : RealToComplex(csh,fsh);
397 48 : Vector<Complex> sh(nPadChan_);
398 24 : convertArray(sh,csh); // downcovert from DComplex to Complex
399 :
400 : // Apply to each elem, corr of shifted
401 264 : for (Int ielem=0;ielem<nElem_;++ielem) {
402 720 : for (Int icorr=0;icorr<nCorr_;++icorr) {
403 960 : Vector<Complex> v(oDupShifted.xyPlane(ielem).row(icorr));
404 480 : v*=sh;
405 : }
406 : }
407 :
408 : // Finally, add it into this
409 24 : Vpad_+=oDupShifted;
410 :
411 24 : }
412 :
413 12 : void DelayFFT::searchPeak() {
414 :
415 12 : delay_.resize(nCorr_,nElem_);
416 12 : delay_.set(0.0);
417 12 : flag_.resize(nCorr_,nElem_);
418 12 : flag_.set(true); // all flagged
419 24 : Vector<Float> amp,pha;
420 : Int ipk;
421 : Float alo,amax,ahi,fpk;
422 132 : for (Int ielem=0;ielem<nElem_;++ielem) {
423 120 : if (ielem==refant_)
424 12 : continue;
425 324 : for (Int icorr=0;icorr<nCorr_;++icorr) {
426 216 : amp=amplitude(Vpad_(Slice(icorr,1,1),Slice(),Slice(ielem,1,1)));
427 : //pha= phase(Vpad_(Slice(icorr,1,1),Slice(),Slice(ielem,1,1)))*(180.0/C::pi);
428 216 : amax=-1.0;
429 216 : ipk=0;
430 15768 : for (Int ich=0;ich<nPadChan_;++ich) {
431 15552 : if (amp[ich]>amax) {
432 630 : ipk=ich;
433 630 : amax=amp[ich];
434 : }
435 : }
436 :
437 216 : Int ilo(ipk>0 ? ipk-1 : (nPadChan_-1));
438 216 : Int ihi(ipk<(nPadChan_-1) ? ipk+1 : 0);
439 216 : alo=amp(ilo);
440 216 : ahi=amp(ihi);
441 :
442 216 : Float denom=(alo-2.*amax+ahi);
443 216 : if (amax>0.0 && abs(denom)>0.0) {
444 216 : fpk=Float(ipk)+0.5-(ahi-amax)/denom;
445 216 : Float delay=fpk/Float(nPadChan_); // cycles/sample
446 216 : if (delay>0.5) delay-=1.0; // fold
447 216 : delay/=df_; // nsec
448 :
449 : /*
450 : if (ielem<8 && icorr>-1) {
451 : cout << endl << ielem << "-" << icorr << " nPadChan_=" << nPadChan_ << " df_=" << df_ << endl
452 : << "peak indices=[" << ilo << ", " << ipk << ", " << ihi << "]-->"
453 : << fpk << endl
454 : << " delay=["
455 : << Float(ilo)/Float(nPadChan_)/df_ << ", "
456 : << Float(ipk)/Float(nPadChan_)/df_ << ", "
457 : << Float(ihi)/Float(nPadChan_)/df_ << "]-->"
458 : << delay << " " << endl
459 : << "amps=" << alo << " " << amax << " " << ahi << " (rel="
460 : << alo/amax << " " << amax/amax << " " << ahi/amax << ")"
461 : << endl
462 : << "phases=" << pha(ilo) << "..." << pha(ipk) << "..." << pha(ihi)
463 : << endl;
464 : if (ipk==0)
465 : cout << amp(Slice(nPadChan_-5,5,1)) << amp(Slice(0,5,1)) << endl;
466 : }
467 : */
468 :
469 216 : delay_(icorr,ielem)=delay;
470 216 : flag_(icorr,ielem)=false;
471 : }
472 : }
473 : }
474 :
475 : // delays for ants>refant must be negated!
476 12 : if (refant_>-1 && refant_<(nElem_-1)) { // at least 1 such ant
477 12 : Matrix<Float> d2(delay_(Slice(),Slice(refant_+1,nElem_-refant_-1,1)));
478 12 : d2*=-1.0f;
479 : }
480 :
481 : // refant (if specified) is unflagged
482 12 : if (refant_>-1)
483 12 : flag_(Slice(),Slice(refant_,1,1))=false;
484 :
485 :
486 12 : }
487 :
488 0 : void DelayFFT::state() {
489 :
490 0 : cout << "DelayFFT::state()..." << endl << " ";
491 0 : cout << " f0_=" << f0_
492 0 : << " df_=" << df_
493 0 : << " padBW_=" << padBW_ << endl << " "
494 0 : << " nCorr_=" << nCorr_
495 0 : << " nPadChan_=" << nPadChan_
496 0 : << " nElem_=" << nElem_
497 0 : << " refant_=" << refant_ << endl << " "
498 0 : << " Vpad_.shape()=" << Vpad_.shape()
499 0 : << " delay_.shape()=" << delay_.shape()
500 0 : << " flag_.shape()=" << flag_.shape()
501 0 : << endl;
502 :
503 0 : }
504 :
505 : // **********************************************************
506 : // CrossDelayFFT Implementations
507 : //
508 :
509 : // Construct from freq info and a data-like Cube<Complex>
510 0 : CrossDelayFFT::CrossDelayFFT(Double f0, Double df, Double padBW,
511 0 : Cube<Complex> V) :
512 0 : DelayFFT(f0,df,padBW,V)
513 0 : {}
514 :
515 :
516 : // Construct from freq info and shape, w/ initialization
517 0 : CrossDelayFFT::CrossDelayFFT(Double f0, Double df, Double padBW) :
518 0 : DelayFFT(f0,df,padBW,1,1,-1,Complex(0.0))
519 0 : {}
520 :
521 : // Construct from a (single) SDB
522 : // (init DelayFFT via freq info/shape ctor, and fill data _here_
523 0 : CrossDelayFFT::CrossDelayFFT(SolveDataBuffer& sdb,Double padBW) :
524 0 : DelayFFT(sdb.freqs()(0)/1e9,
525 0 : (sdb.freqs()(1)-sdb.freqs()(0))/1e9,
526 : padBW,
527 0 : 1,1,-1,Complex(0.0))
528 : {
529 :
530 0 : AlwaysAssert(Vpad_.shape()(0)==1,AipsError); // Working array only 1 corr
531 0 : AlwaysAssert(Vpad_.shape()(2)==1,AipsError); // Working array only 1 elem
532 :
533 : // SDB facts
534 0 : Int nCorr=sdb.nCorrelations();
535 0 : AlwaysAssert(nCorr==4,AipsError); // Must have full set of correlations!
536 0 : Int nChan=sdb.nChannels();
537 0 : Int nRow=sdb.nRows();
538 :
539 : // this->state();
540 :
541 : // Fill the relevant data
542 0 : Vpad_.set(Complex(0.0));
543 :
544 : // Zero weights of flagged samples
545 : // (first, balance cross-hand flags)
546 0 : Slicer pq(Slice(1,1,1),Slice(),Slice()); // PQ
547 0 : Slicer qp(Slice(2,1,1),Slice(),Slice()); // QP
548 0 : sdb.flagCube()(pq)(sdb.flagCube()(qp))=true;
549 0 : sdb.flagCube()(qp)(sdb.flagCube()(pq))=true;
550 0 : sdb.weightSpectrum()(sdb.flagCube())=0.0f;
551 0 : Cube<Complex> Vsum(2,sdb.nChannels(),1,Complex(0.0));
552 0 : Cube<Float> wtsum(2,sdb.nChannels(),1,0.0f);
553 0 : for (Int irow=0;irow<nRow;++irow) {
554 0 : Int a1(sdb.antenna1()(irow)), a2(sdb.antenna2()(irow));
555 0 : if (!sdb.flagRow()(irow) && a1!=a2) { // not flagged, not an AC
556 :
557 : // Accumulate cross-hands
558 0 : Slicer xh(Slice(1,2,1),Slice(),Slice(irow,1,1));
559 0 : Cube<Complex> Vxh(sdb.visCubeCorrected()(xh));
560 0 : Cube<Float> Wxh(sdb.weightSpectrum()(xh));
561 0 : Vsum=Vsum+Vxh*Wxh;
562 0 : wtsum=wtsum+Wxh;
563 : }
564 : }
565 :
566 : // Combine cross-hands into the padded Cube
567 0 : if (sum(wtsum)>0.0f) {
568 0 : Slicer sl1(Slice(),Slice(0,nChan,1),Slice()); // add to this Vpad_ slice
569 0 : Vpad_(sl1)=Vpad_(sl1)+(Vsum(Slice(0,1,1),Slice(),Slice())+conj(Vsum(Slice(1,1,1),Slice(),Slice())));
570 : }
571 :
572 0 : }
573 :
574 :
575 :
576 :
577 : // **********************************************************
578 : // KJones Implementations
579 : //
580 :
581 0 : KJones::KJones(VisSet& vs) :
582 : VisCal(vs), // virtual base
583 : VisMueller(vs), // virtual base
584 0 : GJones(vs) // immediate parent
585 : {
586 0 : if (prtlev()>2) cout << "K::K(vs)" << endl;
587 :
588 : // Extract per-spw ref Freq for phase(delay) calculation
589 : // TBD: these should be in the caltable!!
590 0 : MSSpectralWindow msSpw(vs.spectralWindowTableName());
591 0 : MSSpWindowColumns msCol(msSpw);
592 0 : msCol.refFrequency().getColumn(KrefFreqs_,true);
593 0 : KrefFreqs_/=1.0e9; // in GHz
594 :
595 0 : }
596 :
597 0 : KJones::KJones(String msname,Int MSnAnt,Int MSnSpw) :
598 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
599 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
600 0 : GJones(msname,MSnAnt,MSnSpw) // immediate parent
601 : {
602 0 : if (prtlev()>2) cout << "K::K(msname,MSnAnt,MSnSpw)" << endl;
603 :
604 : // Extract per-spw ref Freq for phase(delay) calculation
605 : // TBD: these should be in the caltable!!
606 :
607 : /* Apparently deprecated; we get it from the ct_ upon setApply...
608 :
609 : NEEDED IN SOLVE CONTEXT?
610 :
611 : MSSpectralWindow msSpw(vs.spectralWindowTableName());
612 : MSSpWindowColumns msCol(msSpw);
613 : msCol.refFrequency().getColumn(KrefFreqs_,true);
614 : KrefFreqs_/=1.0e9; // in GHz
615 : */
616 :
617 :
618 0 : }
619 :
620 15 : KJones::KJones(const MSMetaInfoForCal& msmc) :
621 : VisCal(msmc), // virtual base
622 : VisMueller(msmc), // virtual base
623 15 : GJones(msmc) // immediate parent
624 : {
625 15 : if (prtlev()>2) cout << "K::K(msmc)" << endl;
626 15 : }
627 :
628 0 : KJones::KJones(const Int& nAnt) :
629 : VisCal(nAnt),
630 : VisMueller(nAnt),
631 0 : GJones(nAnt)
632 : {
633 0 : if (prtlev()>2) cout << "K::K(nAnt)" << endl;
634 0 : }
635 :
636 25 : KJones::~KJones() {
637 15 : if (prtlev()>2) cout << "K::~K()" << endl;
638 25 : }
639 :
640 4 : void KJones::setApply(const Record& apply) {
641 :
642 : // Call parent to do conventional things
643 4 : GJones::setApply(apply);
644 :
645 4 : if (calWt())
646 4 : logSink() << " (" << this->typeName() << ": Enforcing calWt()=false for phase/delay-like terms)" << LogIO::POST;
647 :
648 : // Enforce calWt() = false for delays
649 4 : calWt()=false;
650 :
651 : // Extract per-spw ref Freq for phase(delay) calculation
652 : // from the CalTable
653 8 : MSSpectralWindow ctSpw(ct_->spectralWindow());
654 8 : MSSpWindowColumns ctSpwCol(ctSpw);
655 4 : Int nCalSpws(ctSpw.nrow());
656 :
657 8 : String ctvers=ct_->CASAvers();
658 4 : if (ctvers==String("Unknown") || // pre-5.3.0-80 (no version recorded in table)
659 8 : ctvers==String("5.3.0-100") || // a few pre-release versions with reverted behavior
660 8 : ctvers==String("5.3.0-101") ||
661 8 : ctvers==String("5.3.0-102") ||
662 8 : ctvers==String("5.3.0-103") ||
663 8 : ctvers==String("5.3.0-104") ||
664 16 : ctvers==String("5.3.0-105") ||
665 8 : ctvers==String("5.3.0-106") ) {
666 : // Old-fashioned; use spw edge freq
667 0 : ctSpwCol.refFrequency().getColumn(KrefFreqs_,true);
668 0 : if (typeName()!=String("KMBD Jones") &&
669 0 : typeName()!=String("KAntPos Jones") )
670 : logSink() << LogIO::WARN
671 : << " Found pre-5.3.0 CASA delay cal table; using spw REF_FREQUENCY pivot (usually the edge) for phase(freq) calculation."
672 0 : << LogIO::POST;
673 : }
674 : else {
675 : // Use the "physical" (centroid) frequency, per spw
676 8 : Vector<Double> chanfreq;
677 : // Nominally, there should be nSpw() (MS) reference frequencies,
678 : // but caltable may have more (or less)
679 4 : KrefFreqs_.resize(max(nSpw(),nCalSpws)); KrefFreqs_.set(0.0);
680 : // Fill KrefFreqs_ with as many as caltable can support (nCalSpws)
681 : // assuming identity with MS spw ids, for now (spwmap applied below)
682 : // (if nCalSpws>nSpw(), maybe spwmap will need more than nSpw() spws)
683 16 : for (Int ispw=0;ispw<nCalSpws;++ispw) {
684 12 : ctSpwCol.chanFreq().get(ispw,chanfreq,true); // reshape, if nec.
685 12 : Int nch=chanfreq.nelements();
686 12 : KrefFreqs_(ispw)=chanfreq(nch/2);
687 : }
688 : }
689 :
690 4 : KrefFreqs_/=1.0e9; // in GHz
691 :
692 : // Catch spwmap indices not available in the caltable
693 4 : if (spwMap().nelements()>0 && max(spwMap())>=nCalSpws)
694 0 : throw(AipsError("Specified spwmap includes calibration spws not available in the caltable ("+ct_->tableName()+")"));
695 :
696 : /// Re-assign KrefFreq_ according spwmap (if any)
697 4 : if (spwMap().nelements()>0) {
698 8 : Vector<Double> tmpfreqs;
699 4 : tmpfreqs.assign(KrefFreqs_);
700 4 : AlwaysAssert(spwMap().nelements()<=uInt(nSpw()),AipsError); // Should be guaranteed by now
701 12 : for (uInt ispw=0;ispw<spwMap().nelements();++ispw)
702 8 : if (spwMap()(ispw)>-1)
703 6 : KrefFreqs_(ispw)=tmpfreqs(spwMap()(ispw));
704 : }
705 :
706 :
707 :
708 4 : }
709 0 : void KJones::setApply() {
710 :
711 : // Call parent to do conventional things
712 0 : GJones::setApply();
713 :
714 : // Enforce calWt() = false for delays
715 0 : calWt()=false;
716 :
717 : // Set the ref freqs to something usable
718 0 : KrefFreqs_.resize(nSpw());
719 0 : KrefFreqs_.set(0.0);
720 :
721 0 : }
722 :
723 0 : void KJones::setCallib(const Record& callib,
724 : const MeasurementSet& selms) {
725 :
726 : // Call parent to do conventional things
727 0 : SolvableVisCal::setCallib(callib,selms);
728 :
729 0 : if (calWt())
730 0 : logSink() << " (" << this->typeName() << ": Enforcing calWt()=false for phase/delay-like terms)" << LogIO::POST;
731 :
732 : // Enforce calWt() = false for delays
733 0 : calWt()=false;
734 :
735 : // Extract per-spw ref Freq for phase(delay) calculation
736 : // from the CalTable
737 0 : Int nCTSpw(cpp_->nCTSpw()); // number of spws in _caltable_
738 0 : String ctvers=cpp_->CTCASAvers();
739 0 : if (ctvers==String("Unknown") || // pre-5.3.0-80 (no version recorded in table)
740 0 : ctvers==String("5.3.0-100") || // a few pre-release versions with reverted behavior
741 0 : ctvers==String("5.3.0-101") ||
742 0 : ctvers==String("5.3.0-102") ||
743 0 : ctvers==String("5.3.0-103") ||
744 0 : ctvers==String("5.3.0-104") ||
745 0 : ctvers==String("5.3.0-105") ||
746 0 : ctvers==String("5.3.0-106") ) {
747 0 : KrefFreqs_.assign(cpp_->refFreqIn());
748 0 : if (typeName()!=String("KMBD Jones") &&
749 0 : typeName()!=String("KAntPos Jones") )
750 : logSink() << LogIO::WARN
751 : << " Found pre-5.3.0 CASA K (delay) cal table; using spw REF_FREQUENCY pivot (usually the edge) for phase(freq) calculation."
752 0 : << LogIO::POST;
753 : }
754 : else {
755 : // Extract physical freq
756 : // Nominally, there should be nSpw() (MS) reference frequencies,
757 : // but caltable may have more (or less)
758 0 : KrefFreqs_.resize(max(nSpw(),nCTSpw));
759 0 : KrefFreqs_.set(0.0f);
760 : // Fill KrefFreqs_ with as many as caltable can support (nCTSpw)
761 : // assuming identity with MS spw ids, for now (spwmap applied below)
762 : // (if nCTSpw>nSpw(), maybe spwmap will need more than nSpw() spws)
763 0 : for (Int ispw=0;ispw<nCTSpw;++ispw) {
764 0 : const Vector<Double>& f(cpp_->freqIn(ispw));
765 0 : Int nf=f.nelements();
766 0 : KrefFreqs_[ispw]=f[nf/2]; // center (usually this will be same as [0])
767 : }
768 : }
769 0 : KrefFreqs_/=1.0e9; // In GHz
770 :
771 : // Manage spwmap in interim manner
772 0 : spwMap().assign(callib.asRecord(0).asArrayInt("spwmap"));
773 :
774 : // If more than one callib entry, currently have to assume spwmap uniformity
775 0 : if (spwMap().nelements()>0 && callib.nfields()>3)
776 : logSink() << LogIO::WARN << "Assuming uniform spwmap for frequency pivot point for all callib entries for caltable="
777 0 : << ct_->tableName()
778 0 : << ": spwmap = " << spwMap()
779 0 : << LogIO::POST;
780 :
781 : // Catch spwmap that is too long
782 0 : if (spwMap().nelements()>uInt(nSpw()))
783 0 : throw(AipsError("Specified spwmap has more elements ("+String::toString(spwMap().nelements())+") than the number of spectral windows in the MS ("+String::toString(nSpw())+")."));
784 :
785 : // Catch spwmap indices not available in the caltable
786 0 : if (spwMap().nelements()>0 && max(spwMap())>=nCTSpw)
787 0 : throw(AipsError("Specified spwmap includes calibration spws not available in the caltable ("+Path(ct_->tableName()).baseName()+"); cal spws are all <"+String::toString(nCTSpw)+". "));
788 :
789 : // Re-assign KrefFreq_ according spwmap (if any)
790 0 : if (spwMap().nelements()>0) {
791 0 : Vector<Double> tmpfreqs;
792 0 : tmpfreqs.assign(KrefFreqs_);
793 0 : AlwaysAssert(spwMap().nelements()<=uInt(nSpw()),AipsError); // Should be guaranteed by now
794 0 : for (uInt ispw=0;ispw<spwMap().nelements();++ispw)
795 0 : if (spwMap()(ispw)>-1)
796 0 : KrefFreqs_(ispw)=tmpfreqs(spwMap()(ispw));
797 : }
798 :
799 0 : }
800 :
801 5 : void KJones::setSolve(const Record& solve) {
802 :
803 : // Call parent to do conventional things
804 5 : GJones::setSolve(solve);
805 :
806 : // Trap unspecified refant:
807 5 : if (refant()<0)
808 0 : throw(AipsError("Please specify a good reference antenna (refant) explicitly."));
809 :
810 5 : }
811 :
812 :
813 2 : void KJones::specify(const Record& specify) {
814 :
815 2 : return SolvableVisCal::specify(specify);
816 :
817 : }
818 :
819 12 : void KJones::calcAllJones() {
820 :
821 : //if (prtlev()>6) cout << " VJ::calcAllJones()" << endl;
822 :
823 : // Should handle OK flags in this method, and only
824 : // do Jones calc if OK
825 :
826 24 : Vector<Complex> oneJones;
827 24 : Vector<Bool> oneJOK;
828 24 : Vector<Float> onePar;
829 24 : Vector<Bool> onePOK;
830 :
831 24 : ArrayIterator<Complex> Jiter(currJElem(),1);
832 24 : ArrayIterator<Bool> JOKiter(currJElemOK(),1);
833 24 : ArrayIterator<Float> Piter(currRPar(),1);
834 24 : ArrayIterator<Bool> POKiter(currParOK(),1);
835 :
836 12 : Double phase(0.0);
837 132 : for (Int iant=0; iant<nAnt(); iant++) {
838 :
839 1080 : for (Int ich=0; ich<nChanMat(); ich++) {
840 :
841 960 : oneJones.reference(Jiter.array());
842 960 : oneJOK.reference(JOKiter.array());
843 960 : onePar.reference(Piter.array());
844 960 : onePOK.reference(POKiter.array());
845 :
846 2880 : for (Int ipar=0;ipar<nPar();++ipar) {
847 1920 : if (onePOK(ipar)) {
848 1920 : phase=2.0*C::pi*onePar(ipar)*(currFreq()(ich)-KrefFreqs_(currSpw()));
849 1920 : oneJones(ipar)=Complex(cos(phase),sin(phase));
850 1920 : oneJOK(ipar)=true;
851 : }
852 : }
853 :
854 : // Advance iterators
855 960 : Jiter.next();
856 960 : JOKiter.next();
857 960 : if (freqDepPar()) {
858 0 : Piter.next();
859 0 : POKiter.next();
860 : }
861 :
862 : }
863 : // Step to next antenns's pars if we didn't in channel loop
864 120 : if (!freqDepPar()) {
865 120 : Piter.next();
866 120 : POKiter.next();
867 : }
868 : }
869 12 : }
870 :
871 0 : void KJones::selfSolveOne(VisBuffGroupAcc& vbga) {
872 :
873 : // Forward to MBD solver if more than one VB (more than one spw, probably)
874 0 : if (vbga.nBuf()!=1)
875 : // throw(AipsError("KJones can't process multi-VB vbga."));
876 0 : this->solveOneVBmbd(vbga);
877 :
878 : // otherwise, call the single-VB solver with the first VB in the vbga
879 : else
880 0 : this->solveOneVB(vbga(0));
881 :
882 0 : }
883 :
884 :
885 12 : void KJones::selfSolveOne(SDBList& sdbs) {
886 :
887 : // Forward to MBD solver if more than one VB (more than one spw, probably)
888 12 : if (sdbs.nSDB()!=1)
889 12 : this->solveOneSDBmbd(sdbs);
890 :
891 : // otherwise, call the single-VB solver with the first SDB in the SDBList
892 : else
893 0 : this->solveOneSDB(sdbs(0));
894 :
895 12 : }
896 :
897 :
898 :
899 0 : void KJones::solveOneVBmbd(VisBuffGroupAcc& vbga) {
900 :
901 0 : Int nbuf=vbga.nBuf();
902 :
903 0 : Vector<Int> nch(nbuf,0);
904 0 : Vector<Double> f0(nbuf,0.0);
905 0 : Vector<Double> df(nbuf,0.0);
906 0 : Double flo(1e15),fhi(0.0);
907 :
908 0 : for (Int ibuf=0;ibuf<nbuf;++ibuf) {
909 0 : VisBuffer& vb(vbga(ibuf));
910 0 : Vector<Double>& chf(vb.frequency());
911 0 : nch(ibuf)=vbga(ibuf).nChannel();
912 0 : f0(ibuf)=chf(0)/1.0e9; // GHz
913 0 : df(ibuf)=(chf(1)-chf(0))/1.0e9; // GHz
914 0 : flo=min(flo,f0[ibuf]);
915 0 : fhi=max(fhi,f0[ibuf]+nch[ibuf]*df[ibuf]);
916 : }
917 0 : Double tbw=fhi-flo;
918 : // cout << "tbw = " << tbw << " (" << flo << "-" << fhi << ")" << endl;
919 :
920 0 : Double ptbw=tbw*8; // pad total bw by 8X
921 : // TBD: verifty that all df are factors of tbw
922 :
923 0 : Int nCor=vbga(0).nCorr();
924 :
925 0 : DelayFFT sumfft(f0[0],min(df),ptbw,(nCor>1 ? 2 : 1),vbga(0).numberAnt(),refant(),Complex(0.0));
926 0 : for (Int ibuf=0;ibuf<nbuf;++ibuf) {
927 0 : DelayFFT delfft1(vbga(ibuf),ptbw,refant());
928 0 : delfft1.FFT();
929 0 : delfft1.shift(f0[0]);
930 0 : sumfft.add(delfft1);
931 0 : delfft1.searchPeak();
932 : // cout << ibuf << " "
933 : // << delfft1.delay()(Slice(0,1,1),Slice()) << endl;
934 : }
935 :
936 0 : sumfft.searchPeak();
937 :
938 : // cout << "sum: " << sumfft.delay()(Slice(0,1,1),Slice()) << endl;
939 :
940 : // Keep solutions
941 0 : Matrix<Float> sRP(solveRPar().nonDegenerate(1));
942 0 : sRP=sumfft.delay();
943 0 : Matrix<Bool> sPok(solveParOK().nonDegenerate(1));
944 0 : sPok=(!sumfft.flag());
945 :
946 0 : }
947 :
948 :
949 12 : void KJones::solveOneSDBmbd(SDBList& sdbs) {
950 :
951 12 : Int nbuf=sdbs.nSDB();
952 :
953 24 : Vector<Int> nch(nbuf,0);
954 24 : Vector<Double> f0(nbuf,0.0); // the effective ref frequency for each
955 24 : Vector<Double> df(nbuf,0.0); // negative will mean LSB
956 12 : Double flo(1e15),fhi(0.0); // absolute low and high frequencies (SB-independent)
957 :
958 36 : for (Int ibuf=0;ibuf<nbuf;++ibuf) {
959 24 : SolveDataBuffer& sdb(sdbs(ibuf));
960 48 : Vector<Double> chf(sdb.freqs());
961 24 : nch(ibuf)=sdbs(ibuf).nChannels();
962 24 : f0(ibuf)=chf(0)/1.0e9; // GHz
963 24 : df(ibuf)=(chf(1)-chf(0))/1.0e9; // GHz
964 24 : flo=min(flo, min(chf));
965 24 : fhi=max(fhi, max(chf));
966 : }
967 12 : flo/=1e9; // GHz
968 12 : fhi/=1e9; // GHz
969 :
970 12 : Double sb(0.0);
971 12 : if (allLT(df,0.0)) sb=-1.0; // All LSB
972 12 : if (allGT(df,0.0)) sb=1.0; // All USB
973 12 : if (sb==0.0) // Mixed! (FORBIDDEN!)
974 0 : throw(AipsError("Multi-spw delay solutions cannot currently handle mixed sidebands!"));
975 :
976 : /*
977 : cout << "df=" << df << endl;
978 : cout << "f0=" << f0 << endl;
979 : cout << "nch=" << nch << endl;
980 : */
981 :
982 : // Channel bandwidth extremes
983 12 : Double adfmax=max(abs(df)); // unsigned!
984 12 : Double adfmin=min(abs(df)); // unsigned!
985 12 : Double sdfmin(sb*adfmin); // Signed ( <0.0 for LSB)
986 :
987 : //cout << "adfmax=" << adfmax << endl;
988 :
989 : // Total spanned bandwidth as integral multiple of largest absolute channel bandwidth
990 : // NB: Nominal frequency span may not be integral, eg when spws strictly not on same
991 : // grid and/or with different channel bandwidths
992 : // NB: Assumes max channel bandwidth is an integral multiple of all narrower ones.
993 : // This is generally true, and is more strictly enforced in the addWithDupAndShift
994 : // call below....
995 12 : Double tbw=floor(2.0+(fhi-flo)/adfmax)*adfmax;
996 :
997 : // Pad the total bandwith 8X
998 12 : Double ptbw=tbw*8;
999 : // TBD: verifty that all df are factors of tbw
1000 :
1001 : /*
1002 : cout << "tbw = " << tbw << " (" << flo << "-" << fhi << ")" << " resoln=" << 1.0/tbw
1003 : << "; ptbw = " << ptbw << " resoln=" << 1/ptbw
1004 : << endl;
1005 : */
1006 :
1007 12 : Int nCor=sdbs(0).nCorrelations();
1008 :
1009 24 : DelayFFT sumfft(f0[0],sdfmin,ptbw,(nCor>1 ? 2 : 1),nAnt(),refant(),Complex(0.0));
1010 :
1011 36 : for (Int ibuf=0;ibuf<nbuf;++ibuf) {
1012 : //cout << endl << "ibuf=" << ibuf << "-----------------" << endl;
1013 48 : DelayFFT delfft1(sdbs(ibuf),ptbw,refant(),nAnt());
1014 : //delfft1.state();
1015 24 : delfft1.FFT();
1016 24 : sumfft.addWithDupAndShift(delfft1);
1017 : //delfft1.shift(f0[0]);
1018 : //delfft1.searchPeak();
1019 : //sumfft.add(delfft1);
1020 : }
1021 :
1022 : // cout << endl << "Aggregate---------------" << endl;
1023 : // cout << "sumfft.state()=" << endl;
1024 : // sumfft.state();
1025 12 : sumfft.searchPeak();
1026 :
1027 : // Keep solutions
1028 24 : Matrix<Float> sRP(solveRPar().nonDegenerate(1));
1029 12 : sRP=sumfft.delay();
1030 12 : Matrix<Bool> sPok(solveParOK().nonDegenerate(1));
1031 12 : sPok=(!sumfft.flag());
1032 :
1033 12 : }
1034 :
1035 :
1036 : // Do the FFTs
1037 0 : void KJones::solveOneVB(const VisBuffer& vb) {
1038 :
1039 0 : Int nChan=vb.nChannel();
1040 :
1041 0 : solveRPar()=0.0;
1042 0 : solveParOK()=false;
1043 :
1044 : // cout << "solveRPar().shape() = " << solveRPar().shape() << endl;
1045 : // cout << "vb.nCorr() = " << vb.nCorr() << endl;
1046 : // cout << "vb.corrType() = " << vb.corrType() << endl;
1047 :
1048 : // FFT parallel-hands only
1049 0 : Int nCorr=vb.nCorr();
1050 0 : Int nC= (nCorr>1 ? 2 : 1); // number of parallel hands
1051 0 : Int sC= (nCorr>2 ? 3 : 1); // step by 3 for full pol data
1052 :
1053 : // I/O shapes
1054 0 : Int fact(8);
1055 0 : Int nPadChan=nChan*fact;
1056 :
1057 0 : IPosition ip0=vb.visCube().shape();
1058 0 : IPosition ip1=ip0;
1059 0 : ip1(0)=nC; // the number of correlations to FFT
1060 0 : ip1(1)=nPadChan; // padded channel axis
1061 :
1062 : // I/O slicing
1063 0 : Slicer sl0(Slice(0,nC,sC),Slice(),Slice());
1064 0 : Slicer sl1(Slice(),Slice(nChan*(fact-1)/2,nChan,1),Slice());
1065 :
1066 : // Fill the (padded) transform array
1067 : // TBD: only do ref baselines
1068 0 : Cube<Complex> vpad(ip1);
1069 0 : Cube<Complex> slvis=vb.visCube();
1070 0 : vpad.set(Complex(0.0));
1071 0 : vpad(sl1)=slvis(sl0);
1072 :
1073 : // We will only transform frequency axis of 3D array
1074 0 : Vector<Bool> ax(3,false);
1075 0 : ax(1)=true;
1076 :
1077 : // Do the FFT
1078 0 : ArrayLattice<Complex> c(vpad);
1079 0 : LatticeFFT::cfft(c,ax);
1080 :
1081 : // Find peak in each FFT
1082 0 : Int ipk=0;
1083 0 : Float amax(0.0);
1084 0 : Vector<Float> amp;
1085 :
1086 : // cout << "Time=" << MVTime(refTime()/C::day).string(MVTime::YMD,7)
1087 : // << " Spw=" << currSpw() << ":" << endl;
1088 :
1089 0 : for (Int irow=0;irow<vb.nRow();++irow) {
1090 0 : if (!vb.flagRow()(irow) &&
1091 0 : vb.antenna1()(irow)!=vb.antenna2()(irow) &&
1092 0 : (vb.antenna1()(irow)==refant() ||
1093 0 : vb.antenna2()(irow)==refant()) ) {
1094 :
1095 0 : for (Int icor=0;icor<ip1(0);++icor) {
1096 0 : amp=amplitude(vpad(Slice(icor,1,1),Slice(),Slice(irow,1,1)));
1097 0 : ipk=1;
1098 0 : amax=0;
1099 0 : for (Int ich=1;ich<nPadChan-1;++ich) {
1100 0 : if (amp(ich)>amax) {
1101 0 : ipk=ich;
1102 0 : amax=amp(ich);
1103 : }
1104 : } // ich
1105 :
1106 : /*
1107 : cout << vb.antenna1()(irow) << " " << vb.antenna2()(irow) << " "
1108 : << ntrue(vb.flagCube()(Slice(0,2,3),Slice(),Slice(irow,1,1))) << " "
1109 : << ntrue(vb.flag()(Slice(),Slice(irow,1,1)))
1110 : << endl;
1111 : */
1112 :
1113 : // Derive refined peak (fractional) channel
1114 : // via parabolic interpolation of peak and neighbor channels
1115 :
1116 0 : Vector<Float> amp3(amp(IPosition(1,ipk-1),IPosition(1,ipk+1)));
1117 0 : Float denom(amp3(0)-2.0*amp3(1)+amp3(2));
1118 :
1119 0 : if (amax>0.0 && abs(denom)>0) {
1120 :
1121 0 : Float fipk=Float(ipk)+0.5-(amp3(2)-amp3(1))/denom;
1122 :
1123 : // Handle FFT offset and scale
1124 0 : Float delay=(fipk-Float(nPadChan/2))/Float(nPadChan); // cycles/sample
1125 :
1126 : // Convert to cycles/Hz and then to nsec
1127 0 : Double df=vb.frequency()(1)-vb.frequency()(0);
1128 0 : delay/=df;
1129 0 : delay/=1.0e-9;
1130 :
1131 : // cout << " Antenna ID=";
1132 0 : if (vb.antenna1()(irow)==refant()) {
1133 : // cout << vb.antenna2()(irow)
1134 : // << ", pol=" << icor << " delay(nsec)="<< -delay;
1135 0 : solveRPar()(icor,0,vb.antenna2()(irow))=-delay;
1136 0 : solveParOK()(icor,0,vb.antenna2()(irow))=true;
1137 : }
1138 0 : else if (vb.antenna2()(irow)==refant()) {
1139 : // cout << vb.antenna1()(irow)
1140 : // << ", pol=" << icor << " delay(nsec)="<< delay;
1141 0 : solveRPar()(icor,0,vb.antenna1()(irow))=delay;
1142 0 : solveParOK()(icor,0,vb.antenna1()(irow))=true;
1143 : }
1144 : // cout << " (refant ID=" << refant() << ")" << endl;
1145 :
1146 : /*
1147 : cout << irow << " "
1148 : << vb.antenna1()(irow) << " "
1149 : << vb.antenna2()(irow) << " "
1150 : << icor << " "
1151 : << ipk << " "
1152 : << fipk << " "
1153 : << delay << " "
1154 : << endl;
1155 : */
1156 : } // amax > 0
1157 : /*
1158 : else {
1159 : cout << "No solution found for antenna ID= ";
1160 : if (vb.antenna1()(irow)==refant())
1161 : cout << vb.antenna2()(irow);
1162 : else if (vb.antenna2()(irow)==refant())
1163 : cout << vb.antenna1()(irow);
1164 : cout << " in polarization " << icor << endl;
1165 : }
1166 : */
1167 :
1168 : } // icor
1169 : } // !flagrRow, etc.
1170 :
1171 : } // irow
1172 :
1173 : // Ensure refant has zero delay and is NOT flagged
1174 0 : solveRPar()(Slice(),Slice(),Slice(refant(),1,1)) = 0.0;
1175 0 : solveParOK()(Slice(),Slice(),Slice(refant(),1,1)) = true;
1176 :
1177 : /*
1178 : if (nfalse(solveParOK())>0) {
1179 : cout << "NB: No delay solutions found for antenna IDs: ";
1180 : Int nant=solveParOK().shape()(2);
1181 : for (Int iant=0;iant<nant;++iant)
1182 : if (!(solveParOK()(0,0,iant)))
1183 : cout << iant << " ";
1184 : cout << endl;
1185 : }
1186 : */
1187 :
1188 0 : }
1189 : // Do the FFTs
1190 0 : void KJones::solveOneSDB(SolveDataBuffer& sdb) {
1191 :
1192 0 : Int nChan=sdb.nChannels();
1193 :
1194 0 : solveRPar()=0.0;
1195 0 : solveParOK()=false;
1196 :
1197 : // cout << "solveRPar().shape() = " << solveRPar().shape() << endl;
1198 : // cout << "sdb.nCorrelations() = " << sdb.nCorrelations() << endl;
1199 : // cout << "sdb.correlationTypes() = " << sdb.correlationTypes() << endl;
1200 :
1201 : // FFT parallel-hands only
1202 0 : Int nCorr=sdb.nCorrelations();
1203 0 : Int nC= (nCorr>1 ? 2 : 1); // number of parallel hands
1204 0 : Int sC= (nCorr>2 ? 3 : 1); // step by 3 for full pol data
1205 :
1206 : // I/O shapes
1207 0 : Int fact(8);
1208 0 : Int nPadChan=nChan*fact;
1209 :
1210 0 : IPosition ip0=sdb.visCubeCorrected().shape();
1211 0 : IPosition ip1=ip0;
1212 0 : ip1(0)=nC; // the number of correlations to FFT
1213 0 : ip1(1)=nPadChan; // padded channel axis
1214 :
1215 : // I/O slicing
1216 0 : Slicer sl0(Slice(0,nC,sC),Slice(),Slice());
1217 0 : Slicer sl1(Slice(),Slice(nChan*(fact-1)/2,nChan,1),Slice());
1218 :
1219 : // Fill the (padded) transform array
1220 : // TBD: only do ref baselines
1221 0 : Cube<Complex> vpad(ip1);
1222 0 : Cube<Complex> slvis=sdb.visCubeCorrected();
1223 0 : vpad.set(Complex(0.0));
1224 0 : vpad(sl1)=slvis(sl0);
1225 :
1226 : // We will only transform frequency axis of 3D array
1227 0 : Vector<Bool> ax(3,false);
1228 0 : ax(1)=true;
1229 :
1230 : // Do the FFT
1231 0 : ArrayLattice<Complex> c(vpad);
1232 0 : LatticeFFT::cfft(c,ax);
1233 :
1234 : // Find peak in each FFT
1235 0 : Int ipk=0;
1236 0 : Float amax(0.0);
1237 0 : Vector<Float> amp;
1238 :
1239 : // cout << "Time=" << MVTime(refTime()/C::day).string(MVTime::YMD,7)
1240 : // << " Spw=" << currSpw() << ":" << endl;
1241 :
1242 0 : for (Int irow=0;irow<sdb.nRows();++irow) {
1243 0 : if (!sdb.flagRow()(irow) &&
1244 0 : sdb.antenna1()(irow)!=sdb.antenna2()(irow) &&
1245 0 : (sdb.antenna1()(irow)==refant() ||
1246 0 : sdb.antenna2()(irow)==refant()) ) {
1247 :
1248 0 : for (Int icor=0;icor<ip1(0);++icor) {
1249 0 : amp=amplitude(vpad(Slice(icor,1,1),Slice(),Slice(irow,1,1)));
1250 0 : ipk=1;
1251 0 : amax=0;
1252 0 : for (Int ich=1;ich<nPadChan-1;++ich) {
1253 0 : if (amp(ich)>amax) {
1254 0 : ipk=ich;
1255 0 : amax=amp(ich);
1256 : }
1257 : } // ich
1258 :
1259 : /*
1260 : cout << sdb.antenna1()(irow) << " " << sdb.antenna2()(irow) << " "
1261 : << ntrue(sdb.flagCube()(Slice(0,2,3),Slice(),Slice(irow,1,1))) << " "
1262 : << ntrue(sdb.flag()(Slice(),Slice(irow,1,1)))
1263 : << endl;
1264 : */
1265 :
1266 : // Derive refined peak (fractional) channel
1267 : // via parabolic interpolation of peak and neighbor channels
1268 :
1269 0 : Vector<Float> amp3(amp(IPosition(1,ipk-1),IPosition(1,ipk+1)));
1270 0 : Float denom(amp3(0)-2.0*amp3(1)+amp3(2));
1271 :
1272 0 : if (amax>0.0 && abs(denom)>0) {
1273 :
1274 0 : Float fipk=Float(ipk)+0.5-(amp3(2)-amp3(1))/denom;
1275 :
1276 : // Handle FFT offset and scale
1277 0 : Float delay=(fipk-Float(nPadChan/2))/Float(nPadChan); // cycles/sample
1278 :
1279 : // Convert to cycles/Hz and then to nsec
1280 0 : Double df=sdb.freqs()(1)-sdb.freqs()(0);
1281 0 : delay/=df;
1282 0 : delay/=1.0e-9;
1283 :
1284 : // cout << " Antenna ID=";
1285 0 : if (sdb.antenna1()(irow)==refant()) {
1286 : // cout << vb.antenna2()(irow)
1287 : // << ", pol=" << icor << " delay(nsec)="<< -delay;
1288 0 : solveRPar()(icor,0,sdb.antenna2()(irow))=-delay;
1289 0 : solveParOK()(icor,0,sdb.antenna2()(irow))=true;
1290 : }
1291 0 : else if (sdb.antenna2()(irow)==refant()) {
1292 : // cout << vb.antenna1()(irow)
1293 : // << ", pol=" << icor << " delay(nsec)="<< delay;
1294 0 : solveRPar()(icor,0,sdb.antenna1()(irow))=delay;
1295 0 : solveParOK()(icor,0,sdb.antenna1()(irow))=true;
1296 : }
1297 : // cout << " (refant ID=" << refant() << ")" << endl;
1298 :
1299 : /*
1300 : cout << irow << " "
1301 : << vb.antenna1()(irow) << " "
1302 : << vb.antenna2()(irow) << " "
1303 : << icor << " "
1304 : << ipk << " "
1305 : << fipk << " "
1306 : << delay << " "
1307 : << endl;
1308 : */
1309 : } // amax > 0
1310 : /*
1311 : else {
1312 : cout << "No solution found for antenna ID= ";
1313 : if (vb.antenna1()(irow)==refant())
1314 : cout << vb.antenna2()(irow);
1315 : else if (vb.antenna2()(irow)==refant())
1316 : cout << vb.antenna1()(irow);
1317 : cout << " in polarization " << icor << endl;
1318 : }
1319 : */
1320 :
1321 : } // icor
1322 : } // !flagrRow, etc.
1323 :
1324 : } // irow
1325 :
1326 : // Ensure refant has zero delay and is NOT flagged
1327 0 : solveRPar()(Slice(),Slice(),Slice(refant(),1,1)) = 0.0;
1328 0 : solveParOK()(Slice(),Slice(),Slice(refant(),1,1)) = true;
1329 :
1330 : /*
1331 : if (nfalse(solveParOK())>0) {
1332 : cout << "NB: No delay solutions found for antenna IDs: ";
1333 : Int nant=solveParOK().shape()(2);
1334 : for (Int iant=0;iant<nant;++iant)
1335 : if (!(solveParOK()(0,0,iant)))
1336 : cout << iant << " ";
1337 : cout << endl;
1338 : }
1339 : */
1340 :
1341 0 : }
1342 :
1343 : // **********************************************************
1344 : // KcrossJones Implementations
1345 : //
1346 :
1347 0 : KcrossJones::KcrossJones(VisSet& vs) :
1348 : VisCal(vs), // virtual base
1349 : VisMueller(vs), // virtual base
1350 0 : KJones(vs) // immediate parent
1351 : {
1352 0 : if (prtlev()>2) cout << "Kx::Kx(vs)" << endl;
1353 :
1354 : // Extract per-spw ref Freq for phase(delay) calculation
1355 : // TBD: these should be in the caltable!!
1356 0 : MSSpectralWindow msSpw(vs.spectralWindowTableName());
1357 0 : MSSpWindowColumns msCol(msSpw);
1358 0 : msCol.refFrequency().getColumn(KrefFreqs_,true);
1359 0 : KrefFreqs_/=1.0e9; // in GHz
1360 :
1361 0 : }
1362 :
1363 0 : KcrossJones::KcrossJones(String msname,Int MSnAnt,Int MSnSpw) :
1364 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
1365 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
1366 0 : KJones(msname,MSnAnt,MSnSpw) // immediate parent
1367 : {
1368 0 : if (prtlev()>2) cout << "Kx::Kx(msname,MSnAnt,MSnSpw)" << endl;
1369 :
1370 : // Extract per-spw ref Freq for phase(delay) calculation
1371 : // TBD: these should be in the caltable!!
1372 : /* DEPRECATED, because we get it from ct_?
1373 : MSSpectralWindow msSpw(vs.spectralWindowTableName());
1374 : MSSpWindowColumns msCol(msSpw);
1375 : msCol.refFrequency().getColumn(KrefFreqs_,true);
1376 : KrefFreqs_/=1.0e9; // in GHz
1377 : */
1378 0 : }
1379 :
1380 1 : KcrossJones::KcrossJones(const MSMetaInfoForCal& msmc) :
1381 : VisCal(msmc), // virtual base
1382 : VisMueller(msmc), // virtual base
1383 1 : KJones(msmc) // immediate parent
1384 : {
1385 1 : if (prtlev()>2) cout << "Kx::Kx(msmc)" << endl;
1386 :
1387 : // Extract per-spw ref Freq for phase(delay) calculation
1388 : // TBD: these should be in the caltable!!
1389 : /* DEPRECATED, because we get it from ct_?
1390 : MSSpectralWindow msSpw(vs.spectralWindowTableName());
1391 : MSSpWindowColumns msCol(msSpw);
1392 : msCol.refFrequency().getColumn(KrefFreqs_,True);
1393 : KrefFreqs_/=1.0e9; // in GHz
1394 : */
1395 1 : }
1396 :
1397 0 : KcrossJones::KcrossJones(const Int& nAnt) :
1398 : VisCal(nAnt),
1399 : VisMueller(nAnt),
1400 0 : KJones(nAnt)
1401 : {
1402 0 : if (prtlev()>2) cout << "Kx::Kx(nAnt)" << endl;
1403 0 : }
1404 :
1405 2 : KcrossJones::~KcrossJones() {
1406 1 : if (prtlev()>2) cout << "Kx::~Kx()" << endl;
1407 2 : }
1408 :
1409 0 : void KcrossJones::selfSolveOne(VisBuffGroupAcc& vbga) {
1410 :
1411 : // Trap MBD attempt (NYI)
1412 0 : if (vbga.nBuf()!=1)
1413 0 : throw(AipsError("KcrossJones does not yet support MBD"));
1414 : // this->solveOneVBmbd(vbga);
1415 :
1416 : // otherwise, call the single-VB solver with the first VB in the vbga
1417 : else
1418 0 : this->solveOneVB(vbga(0));
1419 :
1420 0 : }
1421 :
1422 :
1423 4 : void KcrossJones::selfSolveOne(SDBList& sdbs) {
1424 :
1425 : // Do MBD if more than one SolveDataBuffer
1426 4 : if (sdbs.nSDB()>1)
1427 0 : this->solveOneSDBmbd(sdbs);
1428 :
1429 : // otherwise, call the single-SDB solver with the first (and only) SDB in SDBLIST
1430 : else
1431 : //this->solveOneSDBmbd(sdbs);
1432 4 : this->solveOneSDB(sdbs(0));
1433 :
1434 4 : }
1435 :
1436 :
1437 : // Do the FFTs
1438 0 : void KcrossJones::solveOneVB(const VisBuffer& vb) {
1439 :
1440 0 : solveRPar()=0.0;
1441 0 : solveParOK()=false;
1442 :
1443 0 : Int fact(8);
1444 0 : Int nChan=vb.nChannel();
1445 0 : Int nPadChan=nChan*fact;
1446 :
1447 : // Collapse cross-hands over baseline
1448 0 : Vector<Complex> sumvis(nPadChan);
1449 0 : sumvis.set(Complex(0.0));
1450 0 : Vector<Complex> slsumvis(sumvis(Slice(nChan*(fact-1)/2,nChan,1)));
1451 0 : Vector<Float> sumwt(nChan);
1452 0 : sumwt.set(0.0);
1453 0 : for (Int irow=0;irow<vb.nRow();++irow) {
1454 0 : if (!vb.flagRow()(irow) &&
1455 0 : vb.antenna1()(irow)!=vb.antenna2()(irow)) {
1456 :
1457 0 : for (Int ich=0;ich<nChan;++ich) {
1458 :
1459 0 : if (!vb.flag()(ich,irow)) {
1460 : // 1st cross-hand
1461 0 : slsumvis(ich)+=(vb.visCube()(1,ich,irow)*vb.weightMat()(1,irow));
1462 0 : sumwt(ich)+=vb.weightMat()(1,irow);
1463 : // 2nd cross-hand
1464 0 : slsumvis(ich)+=conj(vb.visCube()(2,ich,irow)*vb.weightMat()(2,irow));
1465 0 : sumwt(ich)+=vb.weightMat()(2,irow);
1466 : }
1467 : }
1468 : }
1469 : }
1470 : // Normalize the channelized sum
1471 0 : for (int ich=0;ich<nChan;++ich)
1472 0 : if (sumwt(ich)>0)
1473 0 : slsumvis(ich)/=sumwt(ich);
1474 : else
1475 0 : slsumvis(ich)=Complex(0.0);
1476 :
1477 : /*
1478 : cout << "slsumvis.nelements() = " << slsumvis.nelements() << endl;
1479 : cout << "sumwt = " << sumwt/sumwt(0) << endl;
1480 : cout << "amp = " << amplitude(slsumvis) << endl;
1481 : cout << "phase = " << phase(slsumvis)*180.0/C::pi << endl;
1482 : */
1483 :
1484 : // Do the FFT
1485 0 : ArrayLattice<Complex> c(sumvis);
1486 0 : LatticeFFT::cfft(c,true);
1487 :
1488 : // Find peak in each FFT
1489 0 : Vector<Float> amp=amplitude(sumvis);
1490 :
1491 0 : Int ipk=0;
1492 0 : Float amax(0.0);
1493 0 : for (Int ich=0;ich<nPadChan;++ich) {
1494 0 : if (amp(ich)>amax) {
1495 0 : ipk=ich;
1496 0 : amax=amp(ich);
1497 : }
1498 : } // ich
1499 :
1500 : // Derive refined peak (fractional) channel
1501 : // via parabolic interpolation of peak and neighbor channels
1502 0 : Float fipk=ipk;
1503 : // Interpolate the peak (except at edges!)
1504 0 : if (ipk>0 && ipk<(nPadChan-1)) {
1505 0 : Vector<Float> amp3(amp(IPosition(1,ipk-1),IPosition(1,ipk+1)));
1506 0 : fipk=Float(ipk)+0.5-(amp3(2)-amp3(1))/(amp3(0)-2.0*amp3(1)+amp3(2));
1507 0 : Vector<Float> pha3=phase(sumvis(IPosition(1,ipk-1),IPosition(1,ipk+1)));
1508 : }
1509 :
1510 : // Handle FFT offset and scale
1511 0 : Float delay=(fipk-Float(nPadChan/2))/Float(nPadChan); // cycles/sample
1512 :
1513 : // Convert to cycles/Hz and then to nsec
1514 0 : Double df=vb.frequency()(1)-vb.frequency()(0);
1515 0 : delay/=df;
1516 0 : delay/=1.0e-9;
1517 :
1518 0 : solveRPar()(Slice(0,1,1),Slice(),Slice())=delay;
1519 0 : solveParOK()=true;
1520 :
1521 0 : logSink() << " Time="<< MVTime(refTime()/C::day).string(MVTime::YMD,7)
1522 0 : << " Spw=" << currSpw()
1523 : << " Global cross-hand delay=" << delay << " nsec"
1524 0 : << LogIO::POST;
1525 0 : }
1526 :
1527 : // Do the FFTs
1528 4 : void KcrossJones::solveOneSDB(SolveDataBuffer& sdb) {
1529 :
1530 4 : solveRPar()=0.0;
1531 4 : solveParOK()=false;
1532 :
1533 4 : Int fact(8);
1534 4 : Int nChan=sdb.nChannels();
1535 4 : Int nPadChan=nChan*fact;
1536 :
1537 : // Collapse cross-hands over baseline
1538 8 : Vector<Complex> sumvis(nPadChan);
1539 4 : sumvis.set(Complex(0.0));
1540 8 : Vector<Complex> slsumvis(sumvis(Slice(nChan*(fact-1)/2,nChan,1)));
1541 8 : Vector<Float> sumwt(nChan);
1542 4 : sumwt.set(0.0);
1543 184 : for (Int irow=0;irow<sdb.nRows();++irow) {
1544 360 : if (!sdb.flagRow()(irow) &&
1545 180 : sdb.antenna1()(irow)!=sdb.antenna2()(irow)) {
1546 :
1547 1620 : for (Int ich=0;ich<nChan;++ich) {
1548 1440 : if (!sdb.flagCube()(1,ich,irow)) {
1549 : // 1st cross-hand
1550 1440 : Float& wt(sdb.weightSpectrum()(1,ich,irow));
1551 1440 : slsumvis(ich)+=(sdb.visCubeCorrected()(1,ich,irow)*wt);
1552 1440 : sumwt(ich)+=wt;
1553 : }
1554 1440 : if (!sdb.flagCube()(2,ich,irow)) {
1555 : // 2nd cross-hand
1556 1440 : Float& wt(sdb.weightSpectrum()(2,ich,irow));
1557 1440 : slsumvis(ich)+=conj(sdb.visCubeCorrected()(2,ich,irow)*wt);
1558 1440 : sumwt(ich)+=wt;
1559 : }
1560 : }
1561 : }
1562 : }
1563 : // Normalize the channelized sum
1564 36 : for (int ich=0;ich<nChan;++ich)
1565 32 : if (sumwt(ich)>0)
1566 32 : slsumvis(ich)/=sumwt(ich);
1567 : else
1568 0 : slsumvis(ich)=Complex(0.0);
1569 :
1570 : /*
1571 : cout << "slsumvis.nelements() = " << slsumvis.nelements() << endl;
1572 : cout << "amp = " << amplitude(slsumvis) << endl;
1573 : cout << "phase = " << phase(slsumvis)*180.0/C::pi << endl;
1574 : cout << "sumwt = " << sumwt/sumwt(0) << endl;
1575 : */
1576 :
1577 : // Do the FFT
1578 8 : ArrayLattice<Complex> c(sumvis);
1579 4 : LatticeFFT::cfft(c,true);
1580 :
1581 : // Find peak in each FFT
1582 4 : Vector<Float> amp=amplitude(sumvis);
1583 :
1584 4 : Int ipk=0;
1585 4 : Float amax(0.0);
1586 260 : for (Int ich=0;ich<nPadChan;++ich) {
1587 256 : if (amp(ich)>amax) {
1588 34 : ipk=ich;
1589 34 : amax=amp(ich);
1590 : }
1591 : } // ich
1592 :
1593 : // Derive refined peak (fractional) channel
1594 : // via parabolic interpolation of peak and neighbor channels
1595 4 : Float fipk=ipk;
1596 : // Interpolate the peak (except at edges!)
1597 4 : if (ipk>0 && ipk<(nPadChan-1)) {
1598 9 : Vector<Float> amp3(amp(IPosition(1,ipk-1),IPosition(1,ipk+1)));
1599 3 : fipk=Float(ipk)+0.5-(amp3(2)-amp3(1))/(amp3(0)-2.0*amp3(1)+amp3(2));
1600 6 : Vector<Float> pha3=phase(sumvis(IPosition(1,ipk-1),IPosition(1,ipk+1)));
1601 : }
1602 :
1603 : // Handle FFT offset and scale
1604 4 : Float delay=(fipk-Float(nPadChan/2))/Float(nPadChan); // cycles/sample
1605 :
1606 : // Convert to cycles/Hz and then to nsec
1607 4 : Double df=sdb.freqs()(1)-sdb.freqs()(0);
1608 4 : delay/=df;
1609 4 : delay/=1.0e-9;
1610 :
1611 4 : solveRPar()(Slice(0,1,1),Slice(),Slice())=delay;
1612 4 : solveParOK()=true;
1613 :
1614 8 : logSink() << " Time="<< MVTime(refTime()/C::day).string(MVTime::YMD,7)
1615 4 : << " Spw=" << currSpw()
1616 : << " Global cross-hand delay=" << delay << " nsec"
1617 8 : << LogIO::POST;
1618 4 : }
1619 :
1620 0 : void KcrossJones::solveOneSDBmbd(SDBList& sdbs) {
1621 :
1622 0 : Int nbuf=sdbs.nSDB();
1623 :
1624 0 : Vector<Int> nch(nbuf,0);
1625 0 : Vector<Double> f0(nbuf,0.0);
1626 0 : Vector<Double> df(nbuf,0.0);
1627 0 : Double flo(1e15),fhi(0.0);
1628 :
1629 0 : for (Int ibuf=0;ibuf<nbuf;++ibuf) {
1630 0 : SolveDataBuffer& sdb(sdbs(ibuf));
1631 0 : Vector<Double> chf(sdb.freqs());
1632 0 : nch(ibuf)=sdbs(ibuf).nChannels();
1633 0 : f0(ibuf)=chf(0)/1.0e9; // GHz
1634 0 : df(ibuf)=(chf(1)-chf(0))/1.0e9; // GHz
1635 0 : flo=min(flo,f0[ibuf]);
1636 0 : fhi=max(fhi,f0[ibuf]+nch[ibuf]*df[ibuf]);
1637 : }
1638 0 : Double tbw=fhi-flo;
1639 :
1640 0 : Double ptbw=tbw*8; // pad total bw by 8X
1641 : // TBD: verifty that all df are factors of tbw
1642 :
1643 : /*
1644 : cout << "tbw = " << tbw << " (" << flo << "-" << fhi << ")" << " resoln=" << 1.0/tbw
1645 : << "; ptbw = " << ptbw << " resoln=" << 1/ptbw
1646 : << endl;
1647 : */
1648 :
1649 : // Must always have 4 correlations when doing cross-hand delays
1650 0 : Int nCor=sdbs(0).nCorrelations();
1651 0 : AlwaysAssert(nCor==4,AipsError);
1652 :
1653 0 : CrossDelayFFT sumfft(f0[0],min(df),ptbw);
1654 0 : for (Int ibuf=0;ibuf<nbuf;++ibuf) {
1655 0 : CrossDelayFFT delfft1(sdbs(ibuf),ptbw);
1656 0 : delfft1.FFT();
1657 0 : delfft1.shift(f0[0]);
1658 0 : sumfft.add(delfft1);
1659 0 : delfft1.searchPeak();
1660 : }
1661 :
1662 0 : sumfft.searchPeak();
1663 :
1664 : // Keep solution (same for all antennas in first pol)
1665 0 : solveRPar()(Slice(0,1,1),Slice(),Slice())=sumfft.delay()(0,0);
1666 0 : solveParOK()=true;
1667 :
1668 0 : logSink() << " Time="<< MVTime(refTime()/C::day).string(MVTime::YMD,7)
1669 : //<< " Spw=" << currSpw()
1670 0 : << " Multi-band cross-hand delay=" << sumfft.delay()(0,0) << " nsec"
1671 0 : << LogIO::POST;
1672 :
1673 :
1674 0 : }
1675 :
1676 :
1677 :
1678 : // **********************************************************
1679 : // KMBDJones Implementations
1680 : //
1681 :
1682 0 : KMBDJones::KMBDJones(VisSet& vs) :
1683 : VisCal(vs), // virtual base
1684 : VisMueller(vs), // virtual base
1685 0 : KJones(vs) // immediate parent
1686 : {
1687 0 : if (prtlev()>2) cout << "Kmbd::Kmbd(vs)" << endl;
1688 :
1689 : // For MBD, the ref frequencies are zero
1690 0 : KrefFreqs_.resize(nSpw());
1691 0 : KrefFreqs_.set(0.0);
1692 0 : }
1693 :
1694 0 : KMBDJones::KMBDJones(String msname,Int MSnAnt,Int MSnSpw) :
1695 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
1696 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
1697 0 : KJones(msname,MSnAnt,MSnSpw) // immediate parent
1698 : {
1699 0 : if (prtlev()>2) cout << "Kmbd::Kmbd(msname,MSnAnt,MSnSpw)" << endl;
1700 :
1701 : // For MBD, the ref frequencies are zero
1702 0 : KrefFreqs_.resize(MSnSpw);
1703 0 : KrefFreqs_.set(0.0);
1704 0 : }
1705 :
1706 0 : KMBDJones::KMBDJones(const MSMetaInfoForCal& msmc) :
1707 : VisCal(msmc), // virtual base
1708 : VisMueller(msmc), // virtual base
1709 0 : KJones(msmc) // immediate parent
1710 : {
1711 0 : if (prtlev()>2) cout << "Kmbd::Kmbd(msmc)" << endl;
1712 :
1713 : // For MBD, the ref frequencies are zero
1714 0 : KrefFreqs_.resize(nSpw());
1715 0 : KrefFreqs_.set(0.0);
1716 0 : }
1717 :
1718 0 : KMBDJones::KMBDJones(const Int& nAnt) :
1719 : VisCal(nAnt),
1720 : VisMueller(nAnt),
1721 0 : KJones(nAnt)
1722 : {
1723 :
1724 0 : if (prtlev()>2) cout << "Kmbd::Kmbd(nAnt)" << endl;
1725 : // For MBD, the ref frequencies are zero
1726 : // TBD: these should be in the caltable!!
1727 0 : KrefFreqs_.resize(nSpw());
1728 0 : KrefFreqs_.set(0.0);
1729 :
1730 0 : }
1731 :
1732 0 : KMBDJones::~KMBDJones() {
1733 0 : if (prtlev()>2) cout << "Kmbd::~Kmbd()" << endl;
1734 0 : }
1735 :
1736 :
1737 0 : void KMBDJones::setApply(const Record& apply) {
1738 0 : if (prtlev()>2) cout << "Kmbd::setApply()" << endl;
1739 0 : KJones::setApply(apply);
1740 0 : KrefFreqs_.set(0.0); // MBD is ALWAYS ref'd to zero freq
1741 : logSink() << LogIO::WARN
1742 : << " Found pre-5.3.0 CASA multi-band delay cal table; using zero frequency pivot for phase(freq) calculation."
1743 0 : << LogIO::POST;
1744 :
1745 0 : }
1746 :
1747 :
1748 : // **********************************************************
1749 : // KAntPosJones Implementations
1750 : //
1751 :
1752 0 : KAntPosJones::KAntPosJones(VisSet& vs) :
1753 : VisCal(vs), // virtual base
1754 : VisMueller(vs), // virtual base
1755 : KJones(vs), // immediate parent
1756 : doTrDelCorr_(false),
1757 : userEterm_(0.0),
1758 0 : eterm_(0.0)
1759 : {
1760 0 : if (prtlev()>2) cout << "Kap::Kap(vs)" << endl;
1761 :
1762 0 : epochref_p="UTC";
1763 :
1764 0 : }
1765 :
1766 0 : KAntPosJones::KAntPosJones(String msname,Int MSnAnt,Int MSnSpw) :
1767 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
1768 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
1769 : KJones(msname,MSnAnt,MSnSpw), // immediate parent
1770 : doTrDelCorr_(false),
1771 : userEterm_(0.0),
1772 0 : eterm_(0.0)
1773 : {
1774 0 : if (prtlev()>2) cout << "Kap::Kap(msname,MSnAnt,MSnSpw)" << endl;
1775 :
1776 0 : epochref_p="UTC";
1777 :
1778 0 : }
1779 :
1780 4 : KAntPosJones::KAntPosJones(const MSMetaInfoForCal& msmc) :
1781 : VisCal(msmc), // virtual base
1782 : VisMueller(msmc), // virtual base
1783 : KJones(msmc), // immediate parent
1784 : doTrDelCorr_(false),
1785 : userEterm_(0.0),
1786 4 : eterm_(0.0)
1787 : {
1788 4 : if (prtlev()>2) cout << "Kap::Kap(msmc)" << endl;
1789 :
1790 4 : epochref_p="UTC";
1791 :
1792 4 : }
1793 :
1794 :
1795 0 : KAntPosJones::KAntPosJones(const Int& nAnt) :
1796 : VisCal(nAnt),
1797 : VisMueller(nAnt),
1798 : KJones(nAnt),
1799 : doTrDelCorr_(false),
1800 : userEterm_(0.0),
1801 0 : eterm_(0.0)
1802 : {
1803 0 : if (prtlev()>2) cout << "Kap::Kap(nAnt)" << endl;
1804 0 : }
1805 :
1806 8 : KAntPosJones::~KAntPosJones() {
1807 4 : if (prtlev()>2) cout << "Kap::~Kap()" << endl;
1808 8 : }
1809 :
1810 0 : void KAntPosJones::setApply(const Record& apply) {
1811 :
1812 : // TBD: Handle missing solutions in spw 0?
1813 :
1814 : // Force spwmap to all 0 (antpos is not spw-dep)
1815 : // NB: this is required before calling parents, because
1816 : // SVC::setApply sets up the CTPatchedInterp with spwMap()
1817 0 : logSink() << " (" << this->typeName()
1818 0 : << ": Overriding with spwmap=[0] since " << this->typeName()
1819 : << " is not spw-dependent)"
1820 0 : << LogIO::POST;
1821 0 : spwMap().assign(Vector<Int>(1,0));
1822 :
1823 : // Remove spwmap from record, and pass along to generic code
1824 0 : Record newapply;
1825 0 : newapply=apply;
1826 0 : if (newapply.isDefined("spwmap"))
1827 0 : newapply.removeField("spwmap");
1828 :
1829 : // Call parent to do conventional things
1830 0 : KJones::setApply(newapply);
1831 :
1832 : // Arrange for Trop Del correction, if applicable
1833 : // (force check for CalTable keyword
1834 0 : if (vlaTrDelCorrApplicable(true))
1835 0 : initTrDelCorr();
1836 :
1837 0 : }
1838 :
1839 :
1840 :
1841 0 : void KAntPosJones::setCallib(const Record& callib,
1842 : const MeasurementSet& selms)
1843 : {
1844 :
1845 : // Enforce hardwired spwmap in a new revised callib
1846 0 : Record newcallib;
1847 0 : newcallib.define("calwt",Bool(callib.asBool("calwt")));
1848 0 : newcallib.define("tablename",String(callib.asString("tablename")));
1849 :
1850 : // Loop over separate callib slices, if any
1851 : // (formally, probably not needed, but some pipelines are redundant, e.g., over fields groups)
1852 0 : for (uInt icls=0;icls<callib.nfields();++icls) {
1853 : //cout << icls << " name=" << callib.name(icls) << " type=" << callib.type(icls) << " isRecord=" << (callib.type(icls)==TpRecord) << endl;
1854 0 : if (callib.type(icls)==TpRecord) {
1855 0 : Record thiscls;
1856 0 : String slname(callib.name(icls));
1857 0 : thiscls=callib.asRecord(slname); // copy
1858 0 : thiscls.removeField("spwmap");
1859 0 : thiscls.define("spwmap",Vector<Int>(nSpw(),0));
1860 0 : newcallib.defineRecord(slname,thiscls);
1861 : }
1862 : }
1863 :
1864 : // Call generic to do conventional things
1865 0 : SolvableVisCal::setCallib(newcallib,selms);
1866 :
1867 0 : if (calWt())
1868 0 : logSink() << " (" << this->typeName() << ": Enforcing calWt()=false for phase/delay-like terms)" << LogIO::POST;
1869 :
1870 : // Enforce calWt() = false for delays
1871 0 : calWt()=false;
1872 :
1873 : // Force spwmap to all 0 (antpos is not spw-dep)
1874 : // NB: this is required before calling parents, because
1875 : // SVC::setApply sets up the CTPatchedInterp with spwMap()
1876 0 : logSink() << " (" << this->typeName()
1877 0 : << ": Overriding with spwmap=[0] since " << this->typeName()
1878 : << " is not spw-dependent)"
1879 0 : << LogIO::POST;
1880 0 : spwMap().assign(Vector<Int>(1,0));
1881 :
1882 : // Extract per-spw ref Freq for phase(delay) calculation
1883 : // from the CalTable
1884 0 : KrefFreqs_.assign(cpp_->refFreqIn());
1885 0 : KrefFreqs_/=1.0e9; // in GHz
1886 :
1887 : // Re-assign KrefFreq_ according spwmap (if any)
1888 0 : if (spwMap().nelements()>0) {
1889 0 : Vector<Double> tmpfreqs;
1890 0 : tmpfreqs.assign(KrefFreqs_);
1891 0 : for (uInt ispw=0;ispw<spwMap().nelements();++ispw)
1892 0 : if (spwMap()(ispw)>-1)
1893 0 : KrefFreqs_(ispw)=tmpfreqs(spwMap()(ispw));
1894 : }
1895 :
1896 : // Arrange for Trop Del correction, if applicable
1897 : // (force check for CalTable keyword
1898 0 : if (vlaTrDelCorrApplicable(true))
1899 0 : initTrDelCorr();
1900 :
1901 :
1902 0 : }
1903 :
1904 4 : void KAntPosJones::specify(const Record& specify) {
1905 :
1906 12 : LogMessage message(LogOrigin("KAntPosJones","specify"));
1907 :
1908 8 : Vector<Int> spws;
1909 8 : Vector<Int> antennas;
1910 8 : Vector<Double> parameters;
1911 :
1912 4 : Int Nant(0);
1913 :
1914 : // Handle old VLA rotation, if necessary
1915 4 : Bool doVLARot(false);
1916 8 : Matrix<Double> vlaRot=Rot3D(0,0.0);
1917 4 : if (specify.isDefined("caltype") ) {
1918 8 : String caltype=upcase(specify.asString("caltype"));
1919 4 : if (upcase(caltype).contains("VLA")) {
1920 0 : doVLARot=true;
1921 0 : MPosition vlaCenter;
1922 0 : AlwaysAssert(MeasTable::Observatory(vlaCenter,"VLA"),AipsError);
1923 0 : Double vlalong=vlaCenter.getValue().getLong();
1924 : // vlalong=-107.617722*C::pi/180.0;
1925 0 : cout << "We will rotate specified offsets by VLA longitude = "
1926 0 : << vlalong*180.0/C::pi << endl;
1927 0 : vlaRot=Rot3D(2,vlalong);
1928 : }
1929 : }
1930 :
1931 8 : IPosition ip0(3,0,0,0);
1932 8 : IPosition ip1(3,2,0,0);
1933 :
1934 4 : if (specify.isDefined("antenna")) {
1935 : // TBD: the antennas (in order) identifying the solutions
1936 4 : antennas=specify.asArrayInt("antenna");
1937 : // cout << "antenna indices = " << antennas << endl;
1938 4 : Nant=antennas.nelements();
1939 4 : if (Nant<1) {
1940 : // Use specified values for _all_ antennas implicitly
1941 0 : Nant=1; // For the antenna loop below
1942 0 : ip0(2)=0;
1943 0 : ip1(2)=nAnt()-1;
1944 : }
1945 : else {
1946 : // Point to first antenna
1947 4 : ip0(2)=antennas(0);
1948 4 : ip1(2)=ip0(2);
1949 : }
1950 : }
1951 :
1952 4 : if (specify.isDefined("parameter")) {
1953 : // TBD: the actual cal values
1954 4 : parameters=specify.asArrayDouble("parameter");
1955 : // cout << "parameters = ]" << parameters << "[" << endl;
1956 : }
1957 :
1958 4 : Int npar=parameters.nelements();
1959 :
1960 : // Can't proceed of no parameters were specified
1961 4 : if (npar==0)
1962 0 : throw(AipsError("No antenna position corrections specified!"));
1963 :
1964 : // Must be a multiple of 3
1965 4 : if (npar%3 != 0)
1966 0 : throw(AipsError("For antenna position corrections, 3 parameters per antenna are required."));
1967 :
1968 : // cout << "Shapes = " << parameters.nelements() << " "
1969 : // << Nant*3 << endl;
1970 :
1971 : // cout << "parameters = " << parameters << endl;
1972 :
1973 : // Always _ONLY_ spw=0 for antpos corrections
1974 4 : currSpw()=0;
1975 :
1976 : // Loop over specified antennas
1977 4 : Int ipar(0);
1978 21 : for (Int iant=0;iant<Nant;++iant) {
1979 17 : if (Nant>1)
1980 17 : ip1(2)=ip0(2)=antennas(iant);
1981 :
1982 : // make sure ipar doesn't exceed specified list
1983 17 : ipar=ipar%npar;
1984 :
1985 : // The current 3-vector of position corrections
1986 51 : Vector<Double> apar(parameters(IPosition(1,ipar),IPosition(1,ipar+2)));
1987 :
1988 : // If old VLA, rotate them
1989 17 : if (doVLARot) {
1990 0 : cout << "id = " << antennas(iant) << " " << apar;
1991 0 : apar = product(vlaRot,apar);
1992 0 : cout << "--(rotation VLA to ITRF)-->" << apar << endl;
1993 : }
1994 :
1995 : // Loop over 3 parameters, each antenna
1996 68 : for (Int ipar0=0;ipar0<3;++ipar0) {
1997 51 : ip1(0)=ip0(0)=ipar0;
1998 :
1999 51 : Array<Float> sl(solveAllRPar()(ip0,ip1));
2000 :
2001 : // Acccumulation is addition for ant pos corrections
2002 51 : sl+=Float(apar(ipar0));
2003 51 : ++ipar;
2004 : }
2005 : }
2006 :
2007 : // Store in the memory caltable
2008 : // (currSpw()=0 is the only one we need)
2009 4 : keepNCT();
2010 :
2011 : // Detect if Trop Del correction applicable
2012 4 : if (vlaTrDelCorrApplicable()) {
2013 0 : markCalTableForTrDelCorr();
2014 : }
2015 :
2016 4 : }
2017 :
2018 : // KAntPosJones needs ant-based phase direction and position frame
2019 0 : void KAntPosJones::syncMeta(const VisBuffer& vb) {
2020 :
2021 : // Call parent (sets currTime())
2022 0 : KJones::syncMeta(vb);
2023 :
2024 0 : phasedir_p=vb.msColumns().field().phaseDirMeas(currField());
2025 0 : antpos0_p=vb.msColumns().antenna().positionMeas()(0);
2026 :
2027 0 : if (doTrDelCorr_)
2028 : // capture azel info
2029 0 : azel_.reference(vb.azel(currTime()));
2030 :
2031 0 : }
2032 :
2033 : // KAntPosJones needs ant-based phase direction and position frame
2034 0 : void KAntPosJones::syncMeta2(const vi::VisBuffer2& vb) {
2035 :
2036 : // Call parent (sets currTime())
2037 0 : KJones::syncMeta2(vb);
2038 :
2039 0 : phasedir_p=vb.subtableColumns().field().phaseDirMeas(currField());
2040 0 : antpos0_p=vb.subtableColumns().antenna().positionMeas()(0);
2041 :
2042 :
2043 0 : if (doTrDelCorr_)
2044 : // capture azel info
2045 0 : azel_.reference(vb.azel(currTime()));
2046 :
2047 0 : }
2048 :
2049 :
2050 0 : void KAntPosJones::calcAllJones() {
2051 :
2052 0 : if (prtlev()>6) cout << " Kap::calcAllJones()" << endl;
2053 :
2054 : // The relevant timestamp
2055 0 : MEpoch epoch(Quantity(currTime(),"s"));
2056 0 : epoch.setRefString(epochref_p);
2057 :
2058 : // The frame in which we convert our MBaseline from earth to sky and to uvw
2059 0 : MeasFrame mframe(antpos0_p,epoch,phasedir_p);
2060 :
2061 : // template MBaseline, that will be used in calculations below
2062 0 : MBaseline::Ref mbearthref(MBaseline::ITRF,mframe);
2063 0 : MBaseline mb;
2064 0 : MVBaseline mvb;
2065 0 : mb.set(mvb,mbearthref);
2066 :
2067 : // A converter that takes the MBaseline from earth to sky frame
2068 0 : MBaseline::Ref mbskyref(MBaseline::fromDirType(MDirection::castType(phasedir_p.getRef().getType())));
2069 0 : MBaseline::Convert mbcverter(mb,mbskyref);
2070 :
2071 0 : Double phase(0.0);
2072 0 : for (Int iant=0; iant<nAnt(); iant++) {
2073 :
2074 0 : Vector<Float> rpars(currRPar().xyPlane(iant).column(0));
2075 0 : Vector<Double> dpars(rpars.nelements());
2076 0 : convertArray(dpars,rpars);
2077 :
2078 : // We need the w offset (in direction of source) implied
2079 : // by the antenna position correction
2080 0 : Double dw(0.0);
2081 :
2082 : // Only do complicated calculation if there
2083 : // is a non-zero ant pos error
2084 0 : if (max(abs(rpars))>0.0) {
2085 :
2086 : // The current antenna's error as an MBaseline (earth frame)
2087 0 : mvb=MVBaseline(dpars);
2088 0 : mb.set(mvb,mbearthref);
2089 :
2090 : // Convert to sky frame
2091 0 : MBaseline mbdir = mbcverter(mb);
2092 :
2093 : // Get implied uvw
2094 0 : MVuvw uvw(mbdir.getValue(),phasedir_p.getValue());
2095 :
2096 : // dw is third element
2097 0 : dw=uvw.getVector()(2); // in m
2098 :
2099 : }
2100 : //if (iant==26) cout << iant << " " << dw << " -> ";
2101 :
2102 : // Add on the Tropo Delay correction
2103 0 : if (doTrDelCorr_)
2104 0 : dw+=calcTrDelError(iant); // still in m
2105 : //if (iant==26) cout << dw;
2106 :
2107 0 : if (abs(dw)>0.0) {
2108 :
2109 : // To delay units
2110 0 : dw/=C::c; // to sec
2111 0 : dw*=1.0e9; // to nsec
2112 :
2113 : //if (iant==26) cout << " -> " << dw << endl;
2114 :
2115 : // Form the complex corrections per chan (freq)
2116 0 : for (Int ich=0; ich<nChanMat(); ++ich) {
2117 :
2118 : // NB: currFreq() is in GHz
2119 0 : phase=2.0*C::pi*dw*currFreq()(ich);
2120 0 : currJElem()(0,ich,iant)=Complex(cos(phase),sin(phase));
2121 0 : currJElemOK()(0,ich,iant)=true;
2122 :
2123 : }
2124 : }
2125 : else {
2126 : // No correction
2127 0 : currJElem().xyPlane(iant)=Complex(1.0);
2128 0 : currJElemOK().xyPlane(iant)=true;
2129 : }
2130 :
2131 : } // iant
2132 :
2133 :
2134 0 : }
2135 :
2136 4 : bool KAntPosJones::vlaTrDelCorrApplicable(bool checkCalTable) {
2137 :
2138 : // Nominally OFF
2139 4 : doTrDelCorr_=false;
2140 :
2141 4 : Int nobs(msmc().msmd().nObservations());
2142 : // MUST be VLA
2143 8 : Vector<String> obsnames(msmc().msmd().getObservatoryNames());
2144 7 : for (Int iobs=0;iobs<nobs;++iobs)
2145 4 : if (!obsnames(iobs).contains("VLA")) return false;
2146 :
2147 : // Reach here only if all obs are VLA
2148 :
2149 : // Parse boundary dates
2150 3 : if (MJDlim_.nelements()==0) {
2151 3 : MJDlim_.resize(2);
2152 6 : QuantumHolder qh;
2153 3 : String er;
2154 3 : qh.fromString(er,MJD0);
2155 3 : MJDlim_[0]=qh.asQuantity().getValue();
2156 3 : qh.fromString(er,MJD1);
2157 3 : MJDlim_[1]=qh.asQuantity().getValue();
2158 : }
2159 :
2160 : // Are we in affected date range?
2161 3 : double iMJD=msmc().msmd().getTimeRangesOfObservations()[0].first.getValue().get(); // days
2162 4 : if (iMJD > MJDlim_[0] &&
2163 1 : iMJD < MJDlim_[1]) {
2164 : // TURN IT ON!
2165 0 : doTrDelCorr_=true;
2166 : logSink() << "NB: This EVLA dataset appears to fall within the period" << endl
2167 : << " of semester 16B during which the online tropospheric" << endl
2168 0 : << " delay model was mis-applied." << LogIO::WARN << LogIO::POST;
2169 : }
2170 :
2171 :
2172 : // Check table for user-specified scale, if requested
2173 : // (setApply context)
2174 3 : if (checkCalTable) {
2175 0 : const TableRecord& tr(ct_->keywordSet());
2176 0 : if (tr.isDefined("VLATrDelCorr")) {
2177 0 : userEterm_ =tr.asDouble("VLATrDelCorr");
2178 0 : if (userEterm_==0.0) {
2179 : // keyword value says turn it off!
2180 0 : if (doTrDelCorr_)
2181 : logSink() << "Found VLATrDelCorr=0.0 in the antpos caltable; turning trop delay correction OFF."
2182 0 : << LogIO::WARN << LogIO::POST;
2183 0 : doTrDelCorr_=false;
2184 : }
2185 : else {
2186 0 : if (!doTrDelCorr_)
2187 : // user (via caltable) is insisting, even if we are out of date range
2188 : logSink() << "Found VLATrDelCorr keyword in the antpos caltable; turning trop delay correction ON."
2189 0 : << LogIO::WARN << LogIO::POST;
2190 0 : doTrDelCorr_=true;
2191 : }
2192 :
2193 : }
2194 : else {
2195 : // Keyword not present in caltable == TURN IT OFF
2196 0 : doTrDelCorr_=false;
2197 : logSink() << "No VLATrDelCorr keyword in the antpos caltable; turning trop delay correction OFF."
2198 0 : << LogIO::POST;
2199 : }
2200 : }
2201 :
2202 3 : if (doTrDelCorr_)
2203 : logSink() << "A correction for the online tropospheric delay model error WILL BE APPLIED!"
2204 0 : << LogIO::WARN << LogIO::POST;
2205 :
2206 3 : return doTrDelCorr_;
2207 : }
2208 :
2209 0 : void KAntPosJones::markCalTableForTrDelCorr() {
2210 :
2211 : // Only do this if turned on!
2212 0 : AlwaysAssert(doTrDelCorr_,AipsError);
2213 :
2214 : logSink() << "Marking antpos caltable to turn ON the trop delay correction."
2215 0 : << LogIO::WARN << LogIO::POST;
2216 :
2217 : // Add a Table keyword to signal the correction
2218 0 : TableRecord& tr(ct_->rwKeywordSet());
2219 0 : userEterm_=1.0; // non-zero is meaningful, 1.0 is nominal
2220 0 : tr.define("VLATrDelCorr",userEterm_);
2221 :
2222 0 : }
2223 :
2224 :
2225 0 : void KAntPosJones::initTrDelCorr() {
2226 :
2227 : // Must have turned TrDelCorr on!
2228 0 : AlwaysAssert(doTrDelCorr_,AipsError);
2229 :
2230 : // correction scale factor
2231 0 : eterm_=-1.0e-15; // s/m (nominal)
2232 :
2233 : // Apply user-supplied factor from keyword
2234 0 : if (userEterm_!=1.0) {
2235 0 : logSink() << "Applying user-supplied scalefactor="+String::toString<Double>(userEterm_)+" to the trop delay correction."
2236 0 : << LogIO::WARN << LogIO::POST;
2237 0 : eterm_*=userEterm_;
2238 : }
2239 :
2240 0 : logSink() << "Tropospheric delay error correction coefficient="+String::toString<Double>(eterm_/1e-12)+" (ps/m)"
2241 0 : << LogIO::WARN << LogIO::POST;
2242 :
2243 0 : losDist_.resize(nAnt()); // distance from center
2244 0 : armAz_.resize(nAnt());
2245 :
2246 : // Local VLA geometry
2247 0 : Vector<String> sta(msmc().msmd().getAntennaStations());
2248 0 : Vector<Double> vlapos(msmc().msmd().getObservatoryPosition(0).get("m").getValue("m"));
2249 0 : Vector<MPosition> mpos(msmc().msmd().getAntennaPositions());
2250 :
2251 : // Process each antenna
2252 0 : for (int iant=0;iant<nAnt();++iant) {
2253 :
2254 0 : String pre("");
2255 0 : if (sta(iant).startsWith("EVLA:")) pre="EVLA:";
2256 0 : if (sta(iant).startsWith("VLA:")) pre="VLA:";
2257 :
2258 : // Set VLA arm Az (deg)
2259 0 : if (sta(iant).contains(pre+"N")) {
2260 0 : armAz_(iant) = 355.0;
2261 : }
2262 0 : else if (sta(iant).contains(pre+"E")) {
2263 0 : armAz_(iant) = 115.0;
2264 : }
2265 0 : else if (sta(iant).contains(pre+"W")) {
2266 0 : armAz_(iant) = 236.0;
2267 : }
2268 :
2269 : // Set distance from center (m)
2270 0 : Vector<Double> ipos(mpos(iant).get("m").getValue("m"));
2271 0 : losDist_(iant)=sqrt(square(ipos(0)-vlapos(0))+
2272 0 : square(ipos(1)-vlapos(1))+
2273 0 : square(ipos(2)-vlapos(2))); // m from center
2274 :
2275 : }
2276 :
2277 : // Handle units
2278 0 : armAz_*=(C::pi/180.0); // to rad
2279 0 : losDist_*=eterm_; // to s (light travel time)
2280 0 : losDist_*=C::c; // to m (light travel distance)
2281 0 : }
2282 :
2283 :
2284 :
2285 0 : Double KAntPosJones::calcTrDelError(Int iant) {
2286 :
2287 : // Time-dep part of the calculation
2288 0 : Double dgeo(1.0);
2289 0 : Vector<Double> iazel(azel_(iant).getAngle().getValue());
2290 0 : Double az=iazel(0);
2291 0 : Double el=iazel(1);
2292 0 : dgeo*=cos(az-armAz_(iant));
2293 0 : dgeo/=tan(el);
2294 0 : dgeo/=sin(el);
2295 :
2296 : // if (iant==26) cout << az*180/C::pi << " " << el*180/C::pi << " " << dgeo << endl;
2297 :
2298 0 : return dgeo*losDist_(iant);
2299 :
2300 : }
2301 :
2302 :
2303 :
2304 :
2305 :
2306 :
2307 : } //# NAMESPACE CASA - END
|
