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