Line data Source code
1 : //# LJJones.cc: Implementation of LJJOnes VisCal type
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 :
27 : #include <synthesis/MeasurementComponents/LJJones.h>
28 :
29 : #include <msvis/MSVis/VisBuffer.h>
30 : #include <msvis/MSVis/VisBuffAccumulator.h>
31 : #include <casacore/ms/MeasurementSets/MSColumns.h>
32 : #include <synthesis/MeasurementEquations/VisEquation.h>
33 : #include <synthesis/TransformMachines/Utils.h>
34 : #include <synthesis/MeasurementComponents/SteepestDescentSolver.h>
35 : #include <casacore/casa/Quanta/Quantum.h>
36 : #include <casacore/casa/Quanta/QuantumHolder.h>
37 :
38 : #include <casacore/tables/TaQL/ExprNode.h>
39 :
40 : #include <casacore/casa/Arrays/ArrayMath.h>
41 : #include <casacore/casa/BasicSL/String.h>
42 : #include <casacore/casa/Utilities/Assert.h>
43 : #include <casacore/casa/Exceptions/Error.h>
44 : #include <casacore/casa/OS/Memory.h>
45 : #include <casacore/casa/System/Aipsrc.h>
46 :
47 : #include <sstream>
48 :
49 : #include <casacore/casa/Logging/LogMessage.h>
50 : #include <casacore/casa/Logging/LogSink.h>
51 :
52 : using namespace casacore;
53 : namespace casa { //# NAMESPACE CASA - BEGIN
54 :
55 : #define NEED_UNDERSCORES
56 : #if defined(NEED_UNDERSCORES)
57 : #define leakyantso leakyantso_
58 : #endif
59 :
60 : extern "C" {
61 : void leakyantso(Complex *corgain,
62 : Float *corwt,
63 : Int* nant,
64 : Complex *antgain, Complex *leakage,
65 : Int *mode, Float *oresid, Float *nresid,
66 : Float *scanwt, Int* refant, Int* stat, Int *ByPass);
67 : };
68 : // **********************************************************
69 : // LJJones
70 : //
71 0 : LJJones::LJJones(VisSet& vs):
72 : VisCal(vs),
73 : VisMueller(vs),
74 : SolvableVisJones(vs),
75 : pointPar_(),
76 : // ms_p(0),
77 : vs_p(&vs),
78 : maxTimePerSolution(0), minTimePerSolution(10000000), avgTimePerSolution(0),
79 0 : timer(), polMap_p(), tolerance_p(1e-12), gain_p(0.1), niter_p(500), modelImageName_p("")
80 : {
81 0 : if (prtlev()>2) cout << "LJJ::LJJ(vs)" << endl;
82 0 : pbwp_p = NULL;
83 : // String msName = vs.msName();
84 : // ms_p = new MeasurementSet(msName);
85 0 : setParType(VisCalEnum::REAL);
86 0 : byPass_p=0;
87 0 : }
88 0 : LJJones::LJJones(VisSet& vs, MeasurementSet& ) :
89 : VisCal(vs),
90 : VisMueller(vs),
91 : SolvableVisJones(vs),
92 : pointPar_(),
93 : // ms_p(&ms),
94 : vs_p(&vs),
95 : maxTimePerSolution(0), minTimePerSolution(10000000), avgTimePerSolution(0),
96 0 : timer(), polMap_p(), tolerance_p(1e-12), gain_p(0.1), niter_p(500),modelImageName_p("")
97 : {
98 0 : if (prtlev()>2) cout << "LJJ::LJJ(vs)" << endl;
99 0 : pbwp_p = NULL;
100 : // String msName = vs.msName();
101 : // ms_p = new MeasurementSet(msName);
102 0 : setParType(VisCalEnum::REAL);
103 0 : byPass_p=0;
104 0 : }
105 : //
106 : //-----------------------------------------------------------------------
107 : //
108 0 : LJJones::~LJJones()
109 : {
110 0 : if (prtlev()>2) cout << "LJJ::~LJJ()" << endl;
111 0 : }
112 : //
113 : //-----------------------------------------------------------------------
114 : //
115 0 : void LJJones::makeComplexGrid(TempImage<Complex>& Grid,
116 : PagedImage<Float>& ModelImage,
117 : VisBuffer& vb)
118 : {
119 0 : Vector<Int> whichStokes(0);
120 : CoordinateSystem cimageCoord =
121 : StokesImageUtil::CStokesCoord(//cimageShape,
122 : ModelImage.coordinates(),
123 : whichStokes,
124 0 : StokesImageUtil::CIRCULAR);
125 :
126 0 : Grid.resize(IPosition(ModelImage.ndim(),
127 0 : ModelImage.shape()(0),
128 0 : ModelImage.shape()(1),
129 0 : ModelImage.shape()(2),
130 0 : ModelImage.shape()(3)));
131 :
132 0 : Grid.setCoordinateInfo(cimageCoord);
133 :
134 0 : Grid.setMiscInfo(ModelImage.miscInfo());
135 0 : StokesImageUtil::From(Grid,ModelImage);
136 :
137 0 : if(vb.polFrame()==MSIter::Linear)
138 0 : StokesImageUtil::changeCStokesRep(Grid,StokesImageUtil::LINEAR);
139 0 : else StokesImageUtil::changeCStokesRep(Grid,StokesImageUtil::CIRCULAR);
140 0 : }
141 : //
142 : //-----------------------------------------------------------------------
143 : //
144 0 : void LJJones::setModel(const String& modelImageName)
145 : {
146 0 : ROVisIter& vi(vs_p->iter());
147 0 : VisBuffer vb(vi);
148 :
149 0 : PagedImage<Float> modelImage(modelImageName);
150 0 : makeComplexGrid(targetVisModel_,modelImage,vb);
151 0 : vi.originChunks();
152 0 : vi.origin();
153 0 : pbwp_p->initializeToVis(targetVisModel_,vb);
154 0 : modelImageName_p = modelImageName;
155 : // polMap_p = pbwp_p->getPolMap();
156 : // cout << "Pol Map = " << polMap_p << endl;
157 0 : }
158 : //
159 : //-----------------------------------------------------------------------
160 : //
161 0 : void LJJones::setSolve(const Record& solve)
162 : {
163 : // Int nFacets=1;
164 0 : Long cachesize=200000000;
165 0 : Float paInc=1.0; // 1 deg.
166 0 : String cfCacheDirName("tmpCFCache.dir");
167 : //Bool doPBCorr=false;// applyPointingOffsets=true;
168 0 : if (solve.isDefined("cfcache"))
169 0 : cfCacheDirName=solve.asString("cfcache");
170 0 : if (solve.isDefined("painc"))
171 0 : paInc=solve.asDouble("painc");
172 : //if (solve.isDefined("pbcorr"))
173 : // doPBCorr=solve.asBool("pbcorr");
174 :
175 0 : if (pbwp_p) delete pbwp_p;
176 :
177 0 : pbwp_p = new GridFT(cachesize,16);
178 : // pbwp_p = new nPBWProjectFT(//*ms_p,
179 : // nFacets,
180 : // cachesize,
181 : // cfCacheDirName,
182 : // applyPointingOffsets,
183 : // doPBCorr, //Bool do PB correction before prediction
184 : // 16, //Int tilesize=16
185 : // paInc //Float paSteps=1.0
186 : // );
187 0 : casacore::Quantity patol(paInc,"deg");
188 0 : logSink() << LogOrigin("LJJones","setSolve")
189 : << "Parallactic Angle tolerance set to " << patol.getValue("deg") << " deg"
190 0 : << LogIO::POST;
191 : // pbwp_p->setPAIncrement(patol);
192 : // pbwp_p->setLJJones(this);
193 :
194 : // pbwp_p->setPAIncrement(paInc);
195 : //
196 : // What the HELL does the following correspond to? It's not a syntax error!
197 : // azOff(IPosition(1,nAnt())), elOff(IPosition(1,nAnt()));
198 :
199 : // azOff.resize(IPosition(1,nAnt()));
200 : // elOff.resize(IPosition(1,nAnt()));
201 : // azOff = elOff = 0.0;
202 :
203 0 : pointPar_.resize(nPar(),1,nAnt());
204 0 : pointPar_ = 0;
205 :
206 : // calTableName_="test";
207 : // SolvableVisCal::setSolve(solve);
208 :
209 : // if (solve.isDefined("t"))
210 : // interval()=solve.asFloat("t");
211 0 : if (solve.isDefined("solint"))
212 0 : solint()=solve.asString("solint");
213 :
214 0 : QuantumHolder qhsolint;
215 0 : String error;
216 0 : Quantity qsolint;
217 0 : qhsolint.fromString(error,solint());
218 0 : if (error.length()!=0)
219 0 : throw(AipsError("LJJones::setsolve(): Unrecognized units for solint."));
220 0 : qsolint=qhsolint.asQuantumDouble();
221 0 : if (qsolint.isConform("s"))
222 0 : interval()=qsolint.get("s").getValue();
223 : else
224 : {
225 : // assume seconds
226 0 : interval()=qsolint.getValue();
227 0 : solint()=solint()+"s";
228 : }
229 :
230 0 : if (solve.isDefined("preavg"))
231 0 : preavg()=solve.asFloat("preavg");
232 :
233 0 : if (solve.isDefined("refant")) {
234 0 : refantlist().resize();
235 0 : refantlist()=solve.asArrayInt("refant");
236 : }
237 :
238 0 : if (solve.isDefined("phaseonly"))
239 0 : if (solve.asBool("phaseonly"))
240 0 : apmode()="phaseonly";
241 :
242 0 : if (solve.isDefined("table"))
243 0 : calTableName()=solve.asString("table");
244 :
245 0 : if (solve.isDefined("append"))
246 0 : append()=solve.asBool("append");
247 :
248 : // TBD: Warn if table exists (and append=F)!
249 :
250 : // If normalizable & preavg<0, use inteval for preavg
251 : // (or handle this per type, e.g. D)
252 : // TBD: make a nice log message concerning preavg
253 : // TBD: make this work better with solnboundary par
254 : //
255 0 : if (preavg()<0.0) {
256 0 : if (interval()>0.0)
257 : // use interval
258 0 : preavg()=interval();
259 : else
260 : // scan-based, so max out preavg() to get full-chunk time-average
261 0 : preavg()=DBL_MAX;
262 : }
263 : // This is the solve context
264 0 : setSolved(true);
265 0 : setApplied(false);
266 : // SolvableVisCal::setSolve(solve);
267 0 : cs_ = new CalSet<Complex>(nSpw());
268 0 : };
269 : //
270 : //-----------------------------------------------------------------------
271 : //
272 0 : void LJJones::guessPar(VisBuffer& /*vb*/)
273 : {
274 0 : pointPar_=0;
275 : // solveRPar() = 0;
276 0 : }
277 : //
278 : //-----------------------------------------------------------------------
279 : //
280 0 : Cube<Complex>& LJJones::loadPar()
281 : {
282 0 : return pointPar_;
283 : }
284 : //
285 : //-----------------------------------------------------------------------
286 : //
287 : // Specialized setapply extracts image info ("extracts image info"?)
288 : //
289 0 : void LJJones::setApply(const Record& apply)
290 : {
291 : // Call generic
292 : //
293 : // Yet again, this assumes a complex CalSet!
294 : // SolvableVisCal::setApply(applypar);
295 0 : if (prtlev()>2) cout << "LJJones::setApply(apply)" << endl;
296 :
297 : // Call VisCal version for generic stuff
298 : // Sets the value returned by currSpw().
299 : // Resizes currCPar or currRPar to nPar() x nChanPar() x nElem()
300 : // Resizes currParOK() to nChanPar() x nElem()
301 : // Set currParOK() = true
302 : //
303 0 : VisCal::setApply(apply);
304 :
305 : // Collect Cal table parameters
306 0 : if (apply.isDefined("table")) {
307 0 : calTableName()=apply.asString("table");
308 0 : verifyCalTable(calTableName());
309 : }
310 :
311 0 : if (apply.isDefined("select"))
312 0 : calTableSelect()=apply.asString("select");
313 :
314 : // Does this belong here?
315 0 : if (apply.isDefined("append"))
316 0 : append()=apply.asBool("append");
317 :
318 : // Collect interpolation parameters
319 0 : if (apply.isDefined("interp"))
320 0 : tInterpType()=apply.asString("interp");
321 :
322 : // TBD: move spw to VisCal version?
323 0 : if (apply.isDefined("spwmap")) {
324 0 : Vector<Int> spwmap(apply.asArrayInt("spwmap"));
325 0 : spwMap()(IPosition(1,0),IPosition(1,spwmap.nelements()-1))=spwmap;
326 : }
327 :
328 : // TBD: move interval to VisCal version?
329 0 : if (apply.isDefined("t"))
330 0 : interval()=apply.asFloat("t");
331 :
332 : // This is apply context
333 0 : setApplied(true);
334 0 : setSolved(false);
335 :
336 : // TBD: "Arranging to apply...."
337 :
338 :
339 : // Create CalSet, from table
340 0 : switch(parType())
341 : {
342 0 : case VisCalEnum::COMPLEX:
343 : {
344 0 : cs_ = new CalSet<Complex>(calTableName(),calTableSelect(),nSpw(),nPar(),nElem());
345 : // rcs().initCalTableDesc(typeName(),parType());
346 0 : break;
347 : }
348 0 : default:
349 0 : throw(AipsError("Unsupported parameter type found in LJJones::setapply(). Bailing out."));
350 : }
351 :
352 0 : }
353 : //
354 : //-----------------------------------------------------------------------
355 : //
356 0 : void LJJones::applyCal(VisBuffer& vb, Cube<Complex>& Mout)
357 : {
358 : //
359 : // Inflate model data in VB, Mout references it
360 : // (In this type, model data is always re-calc'd from scratch)
361 : //
362 0 : vb.modelVisCube(true);
363 0 : Mout.reference(vb.modelVisCube());
364 0 : }
365 : //
366 : //-----------------------------------------------------------------------
367 : //
368 0 : void LJJones::differentiate(VisBuffer& /*vb*/,
369 : Cube<Complex>& /*Mout*/,
370 : Array<Complex>& /*dMout*/,
371 : Matrix<Bool>& /*Mflg*/)
372 : {
373 : /*
374 : Int nCorr(2); // TBD
375 : //
376 : // Load the offsets from the internal LJJones storage
377 : // (ultimately change the data structures to not need these copies)
378 : //
379 : IPosition ndx(1);
380 : Cube<Float> pointingOffsets(nPar(),1,nAnt());
381 : for(ndx(0)=0;ndx(0)<nAnt();ndx(0)++)
382 : for(Int j=0;j<nPar();j++)
383 : {
384 : // Use solveRPar()(nPar,0,Ant)
385 : pointingOffsets(j,0,ndx(0)) = solveRPar()(j,0,ndx(0));//pointPar_(0,0,ndx(0));
386 : }
387 : //
388 : // Size differentiated model
389 : //
390 : dMout.resize(IPosition(5,nCorr,nPar(),1,vb.nRow(),2));
391 : //
392 : // Model vis shape must match visibility
393 : //
394 : vb.modelVisCube(false);
395 : Mout.reference(vb.modelVisCube());
396 :
397 : //
398 : // Compute the corrupted model and the derivatives.
399 : //
400 :
401 : // Cube<Complex> dVAz(Mout.shape()), dVEl(Mout.shape());
402 : // pbwp_p->nget(vb, azOff, elOff, Mout,dVAz, dVEl,0,1);
403 :
404 : VisBuffer dAZVB(vb), dELVB(vb);
405 : Cube<Complex> dVAz, dVEl;
406 : dAZVB.modelVisCube().resize(vb.modelVisCube().shape());
407 : dELVB.modelVisCube().resize(vb.modelVisCube().shape());
408 : dVAz.reference(dAZVB.modelVisCube());
409 : dVEl.reference(dELVB.modelVisCube());
410 :
411 : pbwp_p->get(vb, dAZVB, dELVB, pointingOffsets);
412 :
413 : // for(Int i=0;i<vb.modelVisCube().shape()(2);i++)
414 : // {
415 : // cout << i
416 : // << " " << vb.modelVisCube()(0,0,i)
417 : // << " " << vb.modelVisCube()(1,0,i)
418 : // << " " << vb.visCube()(0,0,i)
419 : // << " " << vb.visCube()(1,0,i)
420 : // << " " << vb.flag()(0,i)
421 : // << " " << vb.flag()(1,i)
422 : // << " " << vb.antenna1()(i)
423 : // << " " << vb.antenna2()(i)
424 : // << " " << vb.flagRow()(i)
425 :
426 : // << endl;
427 : // }
428 :
429 : //
430 : // For now, copy the derivatives to the required data structure.
431 : //
432 : for(Int j=0;j<nCorr;j++)
433 : for(Int i=0;i<vb.nRow();i++)
434 : {
435 : dMout(IPosition(5,j,0,0,i,0))=dVAz(j,0,i);
436 : dMout(IPosition(5,j,1,0,i,0))=dVEl(j,0,i);
437 : //
438 : // Not sure if the following is what's needed by the solver
439 : //
440 : dMout(IPosition(5,j,0,0,i,1))=conj(dVAz(j,0,i));
441 : dMout(IPosition(5,j,1,0,i,1))=conj(dVEl(j,0,i));
442 : dMout(IPosition(5,j,0,0,i,1))=(dVAz(j,0,i));
443 : dMout(IPosition(5,j,1,0,i,1))=(dVEl(j,0,i));
444 : }
445 :
446 : Mflg.reference(vb.flag());
447 : */
448 0 : }
449 : //
450 : //-----------------------------------------------------------------------
451 : //
452 0 : void LJJones::differentiate(VisBuffer& /*vb*/,VisBuffer& /*dAZVB*/,VisBuffer& /*dELVB*/,
453 : Matrix<Bool>& /*Mflg*/)
454 : {
455 : /*
456 : //
457 : // Load the offsets from the internal LJJones storage
458 : // (ultimately change the data structures to not need these copies)
459 : //
460 : IPosition ndx(1);
461 : Cube<Float> pointingOffsets(nPar(),1,nAnt());
462 : for(ndx(0)=0;ndx(0)<nAnt();ndx(0)++)
463 : for(Int j=0;j<nPar();j++)
464 : {
465 : pointingOffsets(j,0,ndx(0)) = solveRPar()(0,0,ndx(0));//pointPar_(0,0,ndx(0));
466 : }
467 : //
468 : // Model vis shape must match visibility
469 : //
470 : vb.modelVisCube(false);
471 : //
472 : // Compute the corrupted model and the derivatives.
473 : //
474 : dAZVB = vb;
475 : dELVB = vb;
476 : Cube<Complex> dVAz, dVEl;
477 : dAZVB.modelVisCube().resize(vb.modelVisCube().shape());
478 : dELVB.modelVisCube().resize(vb.modelVisCube().shape());
479 : dVAz.reference(dAZVB.modelVisCube());
480 : dVEl.reference(dELVB.modelVisCube());
481 :
482 : pbwp_p->get(vb, dAZVB, dELVB, pointingOffsets);
483 :
484 : // cout << pointingOffsets << endl;
485 : // for(Int i=0;i<vb.modelVisCube().shape()(2);i++)
486 : // {
487 : // cout << "Model: " << i
488 : // << " " << abs(vb.modelVisCube()(0,0,i))
489 : // << " " << arg(vb.modelVisCube()(0,0,i))*57.295
490 : // // << " " << vb.modelVisCube()(1,0,i)
491 : // << " " << abs(vb.visCube()(0,0,i))
492 : // << " " << arg(vb.visCube()(0,0,i))*57.295
493 : // // << " " << vb.visCube()(1,0,i)
494 : // << " " << vb.flag()(0,i)
495 : // // << " " << vb.flag()(1,i)
496 : // << " " << vb.antenna1()(i)
497 : // << " " << vb.antenna2()(i)
498 : // << " " << vb.uvw()(i)(0)
499 : // << " " << vb.uvw()(i)(1)
500 : // << " " << vb.uvw()(i)(2)
501 : // << " " << vb.flagRow()(i)
502 : // << " " << vb.flagCube()(0,0,i)
503 : // // << " " << vb.flagCube()(1,0,i)
504 : // << " " << vb.weight()(i)
505 : // << endl;
506 : // }
507 : // exit(0);
508 : Mflg.reference(vb.flag());
509 : */
510 0 : }
511 :
512 0 : void LJJones::keep(const Int& slot)
513 : {
514 0 : if (prtlev()>4) cout << " SVC::keep(i)" << endl;
515 :
516 0 : if (slot<cs().nTime(currSpw()))
517 : {
518 0 : cs().fieldId(currSpw())(slot)=currField();
519 0 : cs().time(currSpw())(slot)=refTime();
520 : //
521 : // Only stop-start diff matters
522 : // TBD: change CalSet to use only the interval
523 : // TBD: change VisBuffAcc to calculate exposure properly
524 : //
525 0 : cs().startTime(currSpw())(slot)=0.0;
526 0 : cs().stopTime(currSpw())(slot)=interval();
527 : //
528 : // For now, just make these non-zero:
529 : //
530 0 : cs().iFit(currSpw()).column(slot)=1.0;
531 0 : cs().iFitwt(currSpw()).column(slot)=1.0;
532 0 : cs().fit(currSpw())(slot)=1.0;
533 0 : cs().fitwt(currSpw())(slot)=1.0;
534 :
535 0 : IPosition blc4(4,0, focusChan(),0, slot);
536 0 : IPosition trc4(4,nPar()-1,focusChan(),nElem()-1,slot);
537 0 : cs().par(currSpw())(blc4,trc4).nonDegenerate(3) = solvePar();
538 :
539 0 : IPosition blc3(3,focusChan(),0, slot);
540 0 : IPosition trc3(3,focusChan(),nElem()-1,slot);
541 0 : cs().parOK(currSpw())(blc4,trc4).nonDegenerate(3)= solveParOK();
542 0 : cs().parErr(currSpw())(blc4,trc4).nonDegenerate(3)= solveParErr();
543 0 : cs().parSNR(currSpw())(blc4,trc4).nonDegenerate(3)= solveParSNR();
544 0 : cs().solutionOK(currSpw())(slot) = anyEQ(solveParOK(),true);
545 :
546 : }
547 : else
548 0 : throw(AipsError("SVJ::keep: Attempt to store solution in non-existent CalSet slot"));
549 0 : }
550 : //
551 : //-----------------------------------------------------------------------
552 : //
553 0 : void LJJones::setSolve()
554 : {
555 0 : if (prtlev()>2) cout << "LJJ::setSolve()" << endl;
556 :
557 0 : interval()=10.0;
558 0 : refant()=-1;
559 0 : apmode()="<none>";
560 0 : calTableName()="<none>";
561 :
562 : // This is the solve context
563 0 : setSolved(true);
564 0 : setApplied(false);
565 :
566 : // Create a pristine CalSet
567 : // TBD: move this to inflate()?
568 0 : cs_ = new CalSet<Complex>(nSpw());
569 0 : }
570 : //
571 : //-----------------------------------------------------------------------
572 : //
573 0 : void LJJones::inflate(const Vector<Int>& nChan,const Vector<Int>& startChan,
574 : const Vector<Int>& nSlot)
575 : {
576 0 : if (prtlev()>3) cout << " LJJ::inflate(,,)" << endl;
577 : //
578 : // Size up the CalSet
579 : //
580 0 : cs().resize(nPar(),nChan,nElem(),nSlot);
581 0 : cs().setStartChan(startChan);
582 0 : }
583 : //
584 : //-----------------------------------------------------------------------
585 : //
586 0 : void LJJones::initSolve(VisSet& vs)
587 : {
588 0 : if (prtlev()>2) cout << "LJJ::initSolve(vs)" << endl;
589 :
590 : // Determine solving channelization according to VisSet & type
591 0 : setSolveChannelization(vs);
592 :
593 : // Nominal spwMap in solve is identity
594 0 : spwMap().resize(vs.numberSpw());
595 0 : indgen(spwMap());
596 :
597 : // Inflate the CalSet according to VisSet
598 0 : SolvableVisCal::inflate(vs);
599 :
600 : //
601 0 : cs().initCalTableDesc(typeName(),parType());
602 :
603 : // Size the solvePar arrays
604 0 : initSolvePar();
605 0 : }
606 : //
607 : //-----------------------------------------------------------------------
608 : //
609 0 : void LJJones::initSolvePar()
610 : {
611 0 : if (prtlev()>3) cout << " LJJ::initSolvePar()" << endl;
612 :
613 0 : for (Int ispw=0;ispw<nSpw();++ispw)
614 : {
615 0 : currSpw()=ispw;
616 :
617 : // Int nchan = cs().nChan(spwMap()(ispw));
618 0 : solvePar().resize(nPar(),1,nAnt());
619 0 : solveParOK().resize(nPar(),1,nAnt());
620 0 : solveParErr().resize(nPar(),1,nAnt());
621 :
622 0 : solvePar()=(0.0);
623 0 : solveParOK()=true;
624 :
625 0 : solveParSNR().resize(nPar(),1,nAnt());
626 0 : solveParSNR()=0.0;
627 : }
628 :
629 : // currSpw()=0;
630 0 : }
631 : //
632 : //-----------------------------------------------------------------------
633 : //
634 0 : void LJJones::store()
635 : {
636 0 : if (prtlev()>3) cout << " LJJ::store()" << endl;
637 :
638 0 : if (append())
639 0 : logSink() << "Appending solutions to table: " << calTableName()
640 0 : << LogIO::POST;
641 : else
642 0 : logSink() << "Writing solutions to table: " << calTableName()
643 0 : << LogIO::POST;
644 :
645 0 : cs().store(calTableName(),typeName(),append());
646 0 : }
647 : //
648 : //-----------------------------------------------------------------------
649 : //
650 0 : void LJJones::store(const String& table,const Bool& append)
651 : {
652 0 : if (prtlev()>3) cout << " LJJ::store(table,append)" << endl;
653 :
654 : // Override tablename
655 0 : calTableName()=table;
656 0 : SolvableVisCal::append()=append;
657 :
658 : // Call conventional store
659 0 : store();
660 0 : }
661 : //
662 : //-----------------------------------------------------------------------
663 : //
664 0 : Bool LJJones::verifyForSolve(VisBuffer& vb)
665 : {
666 : // cout << "verifyForSolve..." << endl;
667 :
668 0 : Int nAntForSolveFinal(-1);
669 0 : Int nAntForSolve(0);
670 :
671 : // We will count baselines and weights per ant
672 : // and set solveParOK accordingly
673 0 : Vector<Int> blperant(nAnt(),0);
674 0 : Vector<Double> wtperant(nAnt(),0.0);
675 0 : Vector<Bool> antOK(nAnt(),false);
676 :
677 :
678 0 : while (nAntForSolve!=nAntForSolveFinal)
679 : {
680 0 : nAntForSolveFinal=nAntForSolve;
681 0 : nAntForSolve=0;
682 : // TBD: optimize indexing with pointers in the following
683 0 : blperant=0;
684 0 : wtperant=0.0;
685 0 : for (Int irow=0;irow<vb.nRow();++irow)
686 : {
687 0 : Int a1=vb.antenna1()(irow);
688 0 : Int a2=vb.antenna2()(irow);
689 0 : if (!vb.flagRow()(irow) && a1!=a2)
690 : {
691 0 : if (!vb.flag()(focusChan(),irow))
692 : {
693 0 : blperant(a1)+=1;
694 0 : blperant(a2)+=1;
695 :
696 0 : wtperant(a1)+=Double(sum(vb.weightMat().column(irow)));
697 0 : wtperant(a2)+=Double(sum(vb.weightMat().column(irow)));
698 : }
699 : }
700 : }
701 :
702 0 : antOK=false;
703 0 : for (Int iant=0;iant<nAnt();++iant)
704 : {
705 0 : if (blperant(iant)>3 &&
706 0 : wtperant(iant)>0.0)
707 : {
708 : // This antenna is good, keep it
709 0 : nAntForSolve+=1;
710 0 : antOK(iant)=true;
711 : }
712 : else
713 : {
714 : // This antenna under-represented; flag it
715 0 : vb.flagRow()(vb.antenna1()==iant)=true;
716 0 : vb.flagRow()(vb.antenna2()==iant)=true;
717 : }
718 : }
719 : // cout << "blperant = " << blperant << endl;
720 : // cout << "wtperant = " << wtperant << endl;
721 : // cout << "nAntForSolve = " << nAntForSolve << " " << antOK << endl;
722 : }
723 : // We've converged on the correct good antenna count
724 0 : nAntForSolveFinal=nAntForSolve;
725 :
726 : // Set a priori solution flags
727 0 : solveParOK() = false;
728 0 : for (Int iant=0;iant<nAnt();++iant)
729 0 : if (antOK(iant))
730 : // This ant ok
731 0 : solveParOK().xyPlane(iant) = true;
732 : else
733 : // This ant not ok, set soln to zero
734 0 : solveRPar().xyPlane(iant)=0.0;
735 :
736 : // cout << "antOK = " << antOK << endl;
737 : // cout << "solveParOK() = " << solveParOK() << endl;
738 : // cout << "amp(solvePar()) = " << amplitude(solvePar()) << endl;
739 :
740 0 : return (nAntForSolve>3);
741 : }
742 : //
743 : //-----------------------------------------------------------------------
744 : //
745 0 : void LJJones::postSolveMassage(const VisBuffer& vb)
746 : {
747 0 : Array<Float> sol;
748 0 : Double PA = getPA(vb);
749 : Float dL, dM;
750 0 : IPosition ndx(3,0,0,0);
751 :
752 0 : for(ndx(2)=0;ndx(2)<nAnt();ndx(2)++)
753 : {
754 0 : ndx(0)=0; dL = solveRPar()(ndx);
755 0 : ndx(0)=1; dM = solveRPar()(ndx);
756 0 : ndx(0)=0;
757 0 : solveRPar()(ndx) = dL*cos(PA) - dM*sin(PA);
758 0 : ndx(0)=1;
759 0 : solveRPar()(ndx) = dL*sin(PA) - dM*cos(PA);
760 : }
761 0 : };
762 : //
763 : //-----------------------------------------------------------------------
764 : //
765 0 : void LJJones::printActivity(const Int slotNo, const Int fieldId,
766 : const Int spw, const Int nSolutions)
767 : {
768 : Int nSlots, nMesg;
769 :
770 0 : nSlots = cs().nTime(spw);
771 :
772 0 : Double timeTaken = timer.all();
773 0 : if (maxTimePerSolution < timeTaken) maxTimePerSolution = timeTaken;
774 0 : if (minTimePerSolution > timeTaken) minTimePerSolution = timeTaken;
775 0 : avgTimePerSolution += timeTaken;
776 0 : Double avgT = avgTimePerSolution/(nSolutions>0?nSolutions:1);
777 : //
778 : // Boost the no. of messages printed if the next message, based on
779 : // the average time per solution, is going to appear after a time
780 : // longer than my (SB) patience would permit! The limit of
781 : // patience is set to 10 min.
782 : //
783 0 : Float boost = avgT*printFraction(nSlots)*nSlots/(10*60.0);
784 0 : boost = (boost < 1.0)? 1.0:boost;
785 0 : nMesg = (Int)(nSlots*printFraction(nSlots)/boost);
786 0 : nMesg = (nMesg<1?1:nMesg);
787 :
788 0 : Int tmp=abs(nSlots-slotNo); Bool print;
789 0 : print = false;
790 0 : if ((slotNo == 0) || (slotNo == nSlots-1)) print=true;
791 0 : else if ((tmp > 0 ) && ((slotNo+1)%nMesg ==0)) print=true;
792 0 : else print=false;
793 :
794 0 : if (print)
795 : {
796 0 : Int f = (Int)(100*(slotNo+1)/(nSlots>0?nSlots:1));
797 : logSink()<< LogIO::NORMAL
798 : << "Spw=" << spw << " slot=" << slotNo << "/" << nSlots
799 : << " field=" << fieldId << ". Done " << f << "%"
800 : << " Time taken per solution (max/min/avg): "
801 : << maxTimePerSolution << "/"
802 0 : << (minTimePerSolution<0?1:minTimePerSolution) << "/"
803 0 : << avgT << " sec" << LogIO::POST;
804 : }
805 0 : }
806 : //
807 : //-----------------------------------------------------------------------
808 : //
809 0 : void LJJones::selfGatherAndSolve(VisSet& vs, VisEquation& ve)
810 : {
811 : //
812 : // Inform logger/history
813 : //
814 0 : logSink() << LogOrigin("LJJones", "selfGatherAndSolve") << "Solving for " << typeName()
815 0 : << LogIO::POST;
816 : //
817 : // Arrange for iteration over data - set up the VisIter and the VisBuffer
818 : //
819 0 : Block<Int> columns;
820 0 : if (interval()==0.0)
821 : {
822 0 : columns.resize(5);
823 0 : columns[0]=MS::ARRAY_ID;
824 0 : columns[1]=MS::SCAN_NUMBER;
825 0 : columns[2]=MS::FIELD_ID;
826 0 : columns[3]=MS::DATA_DESC_ID;
827 0 : columns[4]=MS::TIME;
828 : }
829 : else
830 : {
831 0 : columns.resize(4);
832 0 : columns[0]=MS::ARRAY_ID;
833 0 : columns[1]=MS::FIELD_ID;
834 0 : columns[2]=MS::DATA_DESC_ID;
835 0 : columns[3]=MS::TIME;
836 : }
837 0 : vs.resetVisIter(columns,interval());
838 0 : VisIter& vi(vs.iter());
839 0 : VisBuffer vb(vi);
840 :
841 : //
842 : // Make an initial guess for the solutions
843 : //
844 0 : guessPar(vb);
845 0 : initSolve(vs);
846 0 : Vector<Int> islot(vs.numberSpw(),0);
847 0 : Array<Complex> visArray,antgain,leakage;
848 0 : Array<Float> visWts;
849 0 : Int refantNdx=refant()+1, ByPass=0,stat,nant,mode=2; // mode=0 ==> Phase-only, 1 ==> Amp-only, 2 ==> AP
850 0 : Int iPol, outPol, spw, nGood(0);
851 : Float oresid, nresid,scanwt;
852 0 : Bool ppPolFound=false, keepSolutions=false;
853 0 : antgain.set(1.0);
854 0 : leakage.set(0.0);
855 0 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk())
856 : {
857 0 : VisBuffAccumulator avgXij(nAnt(),preavg(),false);
858 0 : spw=vi.spectralWindow();
859 0 : currSpw()=spw;
860 0 : Int n=0;
861 0 : for (vi.origin(); vi.more(); vi++)
862 0 : if (syncSolveMeta(vb,vi.fieldId()))
863 : {
864 : // cerr << "VB: " << vb.spectralWindow() << endl;
865 0 : ve.collapse(vb);
866 0 : if (!freqDepPar()) vb.freqAveCubes();
867 0 : if (modelImageName_p != "")
868 : {
869 0 : pbwp_p->get(vb);
870 0 : vb.visCube() /= vb.modelVisCube();
871 : }
872 0 : avgXij.accumulate(vb);
873 0 : n++;
874 : }
875 0 : avgXij.finalizeAverage();
876 0 : VisBuffer& avgVB(avgXij.aveVisBuff());
877 0 : Int nPol=avgVB.corrType().shape()[0];
878 : // cerr << spwMap() << " " << avgVB.spectralWindow() << endl;
879 : Bool vbOk;
880 0 : if (n==0) vbOk=false;
881 0 : else vbOk = syncSolveMeta(avgVB,vi.fieldId());
882 :
883 0 : if (vbOk)
884 : {
885 : //Bool totalGoodSol(false);
886 0 : antgain.set(1.0); leakage.set(0.0);
887 0 : for (Int ich=nChanPar()-1;ich>-1;--ich)
888 : {
889 0 : keepSolutions=false;
890 0 : for(iPol=0;iPol<nPol;iPol++)
891 : {
892 0 : ppPolFound=false;
893 0 : if (avgVB.corrType()[iPol] == Stokes::RR) {outPol=0; ppPolFound=true;}
894 0 : else if (avgVB.corrType()[iPol] == Stokes::LL) {outPol=2; ppPolFound=true;}
895 0 : if (ppPolFound)
896 : {
897 0 : focusChan()=ich;
898 0 : timer.mark();
899 0 : reformVisibilities(avgVB,ich,iPol,visArray, visWts);
900 :
901 : Bool del;
902 0 : nant = nAnt();
903 0 : antgain.resize(IPosition(1,nant),true);
904 0 : leakage.resize(IPosition(1,nant),true);
905 0 : Complex *visArrayPtr=visArray.getStorage(del),
906 0 : *antgainPtr = antgain.getStorage(del),
907 0 : *leakagePtr = leakage.getStorage(del);
908 0 : Float *visWtsPtr=visWts.getStorage(del);
909 0 : ByPass=0;
910 0 : leakyantso(visArrayPtr, visWtsPtr, &nant, antgainPtr, leakagePtr, &mode,
911 : &oresid, &nresid, &scanwt, &refantNdx, &stat, &ByPass);
912 0 : for (IPosition ndx(1,0); ndx(0)<nant;ndx(0)++)
913 0 : if (!avgVB.msColumns().antenna().flagRow()(ndx(0)))
914 0 : solvePar()(outPol+0,0,ndx(0)) = antgain(ndx);
915 0 : for (IPosition ndx(1,0); ndx(0)<nant;ndx(0)++)
916 0 : if (!avgVB.msColumns().antenna().flagRow()(ndx(0)))
917 0 : solvePar()(outPol+1,0,ndx(0)) = leakage(ndx);
918 : }
919 0 : keepSolutions = (stat>0) || keepSolutions;
920 : } // All Pols
921 0 : if (keepSolutions)
922 : {
923 0 : keep(islot(spw));
924 0 : printActivity(islot(spw),vi.fieldId(),spw,nGood);
925 : }
926 : } // All Chans
927 0 : if (keepSolutions) nGood++;
928 0 : islot(spw)++;
929 : } // vbOk
930 : } // chunks
931 : logSink() << " Found " << nGood << " good "
932 0 : << typeName() << " solutions."
933 0 : << LogIO::POST;
934 0 : store();
935 0 : };
936 : //
937 : //-----------------------------------------------------------------------
938 : //
939 0 : void LJJones::diffResiduals(VisIter& /*vi*/, VisEquation& /*ve*/,
940 : VisBuffer& /*residuals*/,
941 : VisBuffer& /*dAZVB*/,
942 : VisBuffer& /*dELVB*/,
943 : Matrix<Bool>& /*Mflg*/)
944 : {
945 : /*
946 : VisBuffAccumulator resAvgr(nAnt(),preavg(),false),
947 : dRes1Avgr(nAnt(), preavg(), false),
948 : dRes2Avgr(nAnt(), preavg(), false);
949 : VisBuffer vb(vi);
950 : IPosition shp;
951 : //
952 : // Load the offsets from the internal LJJones storage
953 : // (ultimately change the data structures to not need these copies)
954 : //
955 : // Cube<Float> pointingOffsets(2,1,nAnt());
956 : // Int nchan=1;
957 : // Cube<Float> pointingOffsets(nPar(),nchan,nAnt());
958 : // for(Int ant=0;ant<nAnt();ant++)
959 : // for(Int par=0;par<nPar();par++)
960 : // for(Int chan=0;chan<nchan;chan++)
961 : // pointingOffsets(par,chan,ant) = solveRPar()(par,chan,ant);//pointPar_(0,0,ndx(0));
962 :
963 : // cout << "LJJ: " << pointingOffsets << endl;
964 : //
965 : // Model vis shape must match visibility
966 : //
967 : residuals.modelVisCube(false);
968 :
969 : residuals.modelVisCube() = Complex(0,0);
970 : dAZVB.modelVisCube() = dELVB.modelVisCube() = Complex(0,0);
971 : static Int j=0;
972 :
973 : // cout << "LJJ::diff: " << residuals.modelVisCube().shape() << " "
974 : // << vb.visCube().shape() << " "
975 : // << vb.modelVisCube().shape() << " "
976 : // << dAZVB.modelVisCube().shape() << " "
977 : // << dELVB.modelVisCube().shape() << endl;
978 :
979 : for (vi.origin(); vi.more(); vi++)
980 : {
981 : // ve.collapse(vb);
982 : dAZVB = dELVB = vb;
983 : shp = vb.modelVisCube().shape();
984 :
985 : // Use the target VBs as temp. storage as well
986 : // (reduce the max. mem. footprint)
987 : dAZVB.modelVisCube().resize(shp);
988 : dELVB.modelVisCube().resize(shp);
989 : vb.modelVisCube() = dAZVB.modelVisCube() = dELVB.modelVisCube() = Complex(0,0);
990 :
991 : // pbwp_p->get(vb, dAZVB, dELVB, pointingOffsets);
992 : pbwp_p->get(vb, dAZVB, dELVB, solveRPar());
993 : // pbwp_p->get(vb);
994 :
995 : vb.modelVisCube() -= vb.visCube(); // Residual = VModel - VObs
996 :
997 : resAvgr.accumulate(vb);
998 : dRes1Avgr.accumulate(dAZVB);
999 : dRes2Avgr.accumulate(dELVB);
1000 : j++;
1001 : }
1002 :
1003 :
1004 : resAvgr.finalizeAverage();
1005 : dRes1Avgr.finalizeAverage();
1006 : dRes2Avgr.finalizeAverage();
1007 : //
1008 : // First copy the internals of the averaged VisBuffers (i.e, Time,
1009 : // UVW, Weights, etc.)
1010 : //
1011 : residuals = resAvgr.aveVisBuff();
1012 : dAZVB = dRes1Avgr.aveVisBuff();
1013 : dELVB = dRes2Avgr.aveVisBuff();
1014 : //
1015 : // Now resize the modelVisCube() of the target VisBuffers (Not
1016 : // resizing the LHS in LHS=RHS of CASA Arrays must be leading to
1017 : // extra code of this type all over the place)
1018 : //
1019 : // shp = resAvgr.aveVisBuff().modelVisCube().shape();
1020 : //
1021 : // The data cubes...
1022 : //
1023 : residuals.modelVisCube().resize(resAvgr.aveVisBuff().modelVisCube().shape());
1024 : dAZVB.modelVisCube().resize(dRes1Avgr.aveVisBuff().modelVisCube().shape());
1025 : dELVB.modelVisCube().resize(dRes2Avgr.aveVisBuff().modelVisCube().shape());
1026 : // The flag cubes..
1027 : residuals.flagCube().resize(resAvgr.aveVisBuff().flagCube().shape());
1028 : dAZVB.flagCube().resize(dRes1Avgr.aveVisBuff().flagCube().shape());
1029 : dELVB.flagCube().resize(dRes2Avgr.aveVisBuff().flagCube().shape());
1030 : // The flags...
1031 : residuals.flag().resize(resAvgr.aveVisBuff().flag().shape());
1032 : dAZVB.flag().resize(dRes1Avgr.aveVisBuff().flag().shape());
1033 : dELVB.flag().resize(dRes2Avgr.aveVisBuff().flag().shape());
1034 : // The row flags....
1035 : residuals.flagRow().resize(resAvgr.aveVisBuff().flagRow().shape());
1036 : dAZVB.flagRow().resize(dRes1Avgr.aveVisBuff().flagRow().shape());
1037 : dELVB.flagRow().resize(dRes2Avgr.aveVisBuff().flagRow().shape());
1038 :
1039 : residuals.weight().resize(resAvgr.aveVisBuff().weight().shape());
1040 : //
1041 : // Now copy the modelVisCube() from the averaged VisBuffers to the
1042 : // target VisBuffers().
1043 : //
1044 : // The data cubes...
1045 : residuals.modelVisCube() = resAvgr.aveVisBuff().modelVisCube();
1046 : dAZVB.modelVisCube() = dRes1Avgr.aveVisBuff().modelVisCube();
1047 : dELVB.modelVisCube() = dRes2Avgr.aveVisBuff().modelVisCube();
1048 :
1049 : // The flag cubes...
1050 : residuals.flagCube() = resAvgr.aveVisBuff().flagCube();
1051 : dAZVB.flagCube() = dRes1Avgr.aveVisBuff().flagCube();
1052 : dELVB.flagCube() = dRes2Avgr.aveVisBuff().flagCube();
1053 : // The flags...
1054 : residuals.flag() = resAvgr.aveVisBuff().flag();
1055 : dAZVB.flag() = dRes1Avgr.aveVisBuff().flag();
1056 : dELVB.flag() = dRes2Avgr.aveVisBuff().flag();
1057 : // The row flags...
1058 : residuals.flagRow() = resAvgr.aveVisBuff().flagRow();
1059 : dAZVB.flagRow() = dRes1Avgr.aveVisBuff().flagRow();
1060 : dELVB.flagRow() = dRes2Avgr.aveVisBuff().flagRow();
1061 :
1062 : residuals.weight() = resAvgr.aveVisBuff().weight();
1063 : //
1064 : // Average the residuals and the derivates in frequency.
1065 : //
1066 : residuals.freqAveCubes();
1067 : dAZVB.freqAveCubes();
1068 : dELVB.freqAveCubes();
1069 : Mflg.reference(residuals.flag());
1070 : */
1071 0 : }
1072 :
1073 : //
1074 : // Quick-n-dirty implementation - to be replaced by use of CalInterp
1075 : // class if-and-when that's ready for use
1076 : //
1077 0 : void LJJones::nearest(const Double thisTime, Array<Complex>& vals)
1078 : {
1079 0 : Array<Complex> par = getOffsets(currSpw());
1080 0 : Vector<Double> time = getTime(currSpw());
1081 : // Int nant=nAnt();
1082 0 : uInt nTimes = time.nelements(), slot=0;
1083 0 : Double dT=abs(time[0]-thisTime);
1084 0 : IPosition shp=par.shape();
1085 :
1086 0 : for(uInt i=0;i<nTimes;i++)
1087 0 : if (abs(time[i]-thisTime) < dT)
1088 : {
1089 0 : dT = abs(time[i]-thisTime);
1090 0 : slot=i;
1091 : }
1092 :
1093 0 : if (slot >= nTimes) throw(AipsError("LJJones::nearest(): Internal problem - "
1094 0 : "nearest slot is out of range"));
1095 0 : Array<Complex> tmp=par(IPosition(4,0,0,0,slot), IPosition(4,shp[0]-1,shp[1]-1,shp[2]-1,slot));
1096 0 : vals.resize();
1097 0 : vals = tmp;
1098 0 : }
1099 :
1100 0 : void LJJones::printRPar()
1101 : {
1102 0 : Int n=solveRPar().shape()(2);
1103 0 : for(Int i=0;i<n;i++)
1104 0 : cout << solveRPar()(0,0,i) << " "
1105 0 : << solveRPar()(1,0,i) << " "
1106 0 : << solveRPar()(2,0,i) << " "
1107 0 : << solveRPar()(3,0,i)
1108 0 : << endl;
1109 0 : }
1110 :
1111 0 : void LJJones::reformVisibilities(const VisBuffer& vb, const Int whichChan,
1112 : const Int whichPol,
1113 : Array<Complex>& visArray,
1114 : Array<Float>& visWts)
1115 : {
1116 0 : Int nant = nAnt(), nRows = vb.nRow();
1117 0 : Bool useModelVis = (modelImageName_p != "");
1118 0 : IPosition shp(2,nant,nant);
1119 0 : visArray.resize(shp);
1120 0 : visWts.resize(shp);
1121 0 : visArray.set(0.0);
1122 0 : visWts.set(0.0);
1123 0 : for(Int r=0; r<nRows; r++)
1124 : {
1125 0 : Int a1=vb.antenna1()(r), a2=vb.antenna2()(r);
1126 0 : Complex vis=vb.visCube()(whichPol,whichChan,r),
1127 0 : mvis=vb.modelVisCube()(whichPol,whichChan,r);
1128 0 : Float wts = vb.weight()(r);
1129 0 : if (!useModelVis) mvis=1.0;
1130 0 : mvis=1.0;
1131 0 : if (vb.flagRow()(r) || vb.flagCube()(whichPol, whichChan, r))
1132 : {
1133 0 : vis = 0.0;
1134 0 : wts = 0.0;
1135 : }
1136 0 : if ((a1 != a2) && (wts > 0.0))
1137 : {
1138 0 : shp(0)=a1; shp(1)=a2;
1139 0 : visArray(shp) = vis/mvis;
1140 0 : visWts(shp)=wts;
1141 0 : shp(0)=a2; shp(1)=a1;
1142 0 : visArray(shp) = conj(vis/mvis);
1143 0 : visWts(shp)=wts;
1144 : // cout << a1 << " " << a2 << " " << vis << " " << mvis << endl;
1145 : }
1146 : }
1147 : // cout << "----------------------------------------------" << endl;
1148 : // for(shp(0)=0;shp(0)<nant;shp(0)++)
1149 : // {
1150 : // for(shp(1)=0;shp(1)<nant;shp(1)++)
1151 : // cout << abs(visArray(shp)) << " ";
1152 : // cout << endl;
1153 : // }
1154 : // cout << "----------------------------------------------" << endl;
1155 : // for(shp(0)=0;shp(0)<nant;shp(0)++)
1156 : // {
1157 : // for(shp(1)=0;shp(1)<nant;shp(1)++)
1158 : // cout << (visWts(shp)) << " ";
1159 : // cout << endl;
1160 : // }
1161 : // cout << "----------------------------------------------" << endl;
1162 0 : }
1163 :
1164 0 : Int LJJones::getRefAnt(const Array<Float>& visWts)
1165 : {
1166 0 : IPosition shp(visWts.shape());
1167 : Float sumWts;
1168 0 : Int nx=shp[0], ny=shp[1];
1169 0 : for (shp[0]=0;shp[0]<nx; shp[0]++)
1170 : {
1171 0 : sumWts=0;
1172 0 : for (shp[1]=0;shp[1]<ny; shp[1]++)
1173 0 : if (visWts(shp) > 0) sumWts += visWts(shp);
1174 0 : if (sumWts > 0) {return shp[0];}
1175 : }
1176 0 : return -1;
1177 : }
1178 0 : void LJJones::getAvgVB(VisIter& vi, VisEquation& ve, VisBuffer& vb)
1179 : {
1180 0 : VisBuffAccumulator avgXij(nAnt(),preavg(),false);
1181 0 : VisBuffer accumVB(vi);
1182 0 : Int n=0;
1183 :
1184 0 : for (vi.origin(); vi.more(); vi++)
1185 : {
1186 0 : ve.collapse(accumVB);
1187 0 : if (modelImageName_p != "")
1188 : {
1189 0 : pbwp_p->get(accumVB);
1190 0 : accumVB.visCube() /= accumVB.modelVisCube();
1191 : }
1192 0 : avgXij.accumulate(accumVB);
1193 0 : n++;
1194 : }
1195 0 : cerr << "No. of time-stamps average = " << n << endl;
1196 0 : avgXij.finalizeAverage();
1197 0 : vb = avgXij.aveVisBuff();
1198 0 : }
1199 : } //# NAMESPACE CASA - END
1200 :
1201 :
1202 : // cout << pointingOffsets << endl;
1203 : // for(Int i=0;i<vb.modelVisCube().shape()(2);i++)
1204 : // {
1205 : // cout << "Model: " << i
1206 : // << " " << abs(vb.modelVisCube()(0,0,i))
1207 : // << " " << arg(vb.modelVisCube()(0,0,i))*57.295
1208 : // // << " " << vb.modelVisCube()(1,0,i)
1209 : // << " " << abs(vb.visCube()(0,0,i))
1210 : // << " " << arg(vb.visCube()(0,0,i))*57.295
1211 : // // << " " << vb.visCube()(1,0,i)
1212 : // << " " << vb.flag()(0,i)
1213 : // // << " " << vb.flag()(1,i)
1214 : // << " " << vb.antenna1()(i)
1215 : // << " " << vb.antenna2()(i)
1216 : // << " " << vb.uvw()(i)(0)
1217 : // << " " << vb.uvw()(i)(1)
1218 : // << " " << vb.uvw()(i)(2)
1219 : // << " " << vb.flagRow()(i)
1220 : // << " " << vb.flagCube()(0,0,i)
1221 : // // << " " << vb.flagCube()(1,0,i)
1222 : // << " " << vb.weight()(i)
1223 : // << endl;
1224 : // }
1225 : // exit(0);
|
