Line data Source code
1 : //# EVLASwPow.cc: Implementation of EVLA Switched Power Calibration
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: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 :
27 : #include <synthesis/MeasurementComponents/EVLASwPow.h>
28 :
29 : #include <synthesis/MeasurementComponents/MSMetaInfoForCal.h>
30 :
31 : #include <casacore/ms/MeasurementSets/MSColumns.h>
32 : #include <casacore/casa/BasicMath/Math.h>
33 : #include <casacore/tables/Tables/Table.h>
34 : #include <casacore/tables/Tables/TableIter.h>
35 : #include <synthesis/CalTables/CTGlobals.h>
36 :
37 : #include <casacore/casa/BasicSL/String.h>
38 : #include <casacore/casa/Utilities/Assert.h>
39 : #include <casacore/casa/Utilities/GenSort.h>
40 : #include <casacore/casa/Exceptions/Error.h>
41 : #include <casacore/casa/System/Aipsrc.h>
42 : #include <casacore/casa/System/ProgressMeter.h>
43 :
44 : #include <sstream>
45 :
46 : #include <casacore/casa/Logging/LogMessage.h>
47 : #include <casacore/casa/Logging/LogSink.h>
48 :
49 :
50 :
51 : using namespace casacore;
52 : namespace casa { //# NAMESPACE CASA - BEGIN
53 :
54 :
55 : // **********************************************************
56 : // EVLASwPow Implementations
57 : //
58 :
59 1 : EVLASwPow::SPType EVLASwPow::sptype(const String name) {
60 1 : String utype=upcase(name);
61 1 : if (utype.contains("SWP/RQ"))
62 0 : return EVLASwPow::SWPOVERRQ;
63 1 : if (utype.contains("SWPOW"))
64 0 : return EVLASwPow::SWPOW;
65 1 : if (utype.contains("EVLAGAIN"))
66 0 : return EVLASwPow::SWPOW;
67 1 : if (utype.contains("RQ"))
68 0 : return EVLASwPow::RQ;
69 1 : if (utype.contains("SWPWTS"))
70 1 : return EVLASwPow::SWPWTS;
71 :
72 : // Only get here if name unrecognized
73 0 : throw(AipsError(name+" is not among recognized EVLA Switched Power types ('swpow','evlagain','rq','swp/rq', 'swpwts')"));
74 :
75 : // Should never reach here, but this is accurate (and avoids compiler warning)
76 : return EVLASwPow::NONE;
77 :
78 1 : }
79 :
80 2 : String EVLASwPow::sptype(EVLASwPow::SPType sptype) {
81 2 : switch (sptype) {
82 0 : case EVLASwPow::SWPOW: {
83 0 : return String("swpow");
84 : break;
85 : }
86 0 : case EVLASwPow::RQ: {
87 0 : return String("rq");
88 : break;
89 : }
90 0 : case EVLASwPow::SWPOVERRQ: {
91 0 : return String("swpow/rq");
92 : break;
93 : }
94 2 : case EVLASwPow::SWPWTS: {
95 2 : return String("swpwts");
96 : break;
97 : }
98 0 : case EVLASwPow::NONE:
99 : default: {
100 0 : throw(AipsError("Unrecognized EVLA Switched Power type"));
101 : }
102 : }
103 : // Should never reach here, but this is accurate (and avoids compiler warning)
104 : return String("None");
105 : }
106 :
107 0 : EVLASwPow::EVLASwPow(VisSet& vs) :
108 : VisCal(vs), // virtual base
109 : VisMueller(vs), // virtual base
110 : GJones(vs), // immediate parent
111 0 : sysPowTabName_(vs.syspowerTableName()),
112 0 : calDevTabName_(vs.caldeviceTableName()),
113 0 : correff_(Float(0.932)), // EVLA-specific net corr efficiency (4bit)
114 0 : frgrotscale_(Float(1.176)), // EVLA-specific fringe rotation mean _scale_
115 : //nyquist_(1.0),
116 0 : effChBW_()
117 :
118 : {
119 0 : if (prtlev()>2) cout << "EVLASwPow::EVLASwPow(vs)" << endl;
120 :
121 0 : nChanParList().set(1);
122 0 : startChanList().set(0);
123 :
124 : // Get spw total bandwidths
125 0 : const MSSpWindowColumns& spwcols = vs.iter().msColumns().spectralWindow();
126 0 : effChBW_.resize(nSpw());
127 0 : for (Int ispw=0;ispw<nSpw();++ispw)
128 0 : effChBW_(ispw)=Vector<Double>(spwcols.effectiveBW()(0))(0);
129 :
130 0 : }
131 :
132 0 : EVLASwPow::EVLASwPow(String msname,Int MSnAnt,Int MSnSpw) :
133 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
134 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
135 : GJones(msname,MSnAnt,MSnSpw), // immediate parent
136 0 : sysPowTabName_(""),
137 0 : calDevTabName_(""),
138 0 : correff_(Float(0.932)), // EVLA-specific net corr efficiency (4bit)
139 0 : frgrotscale_(Float(1.176)), // EVLA-specific fringe rotation mean _scale_
140 : //nyquist_(1.0),
141 0 : effChBW_()
142 :
143 : {
144 0 : if (prtlev()>2) cout << "EVLASwPow::EVLASwPow(msname,MSnAnt,MSnSpw)" << endl;
145 :
146 0 : nChanParList().set(1);
147 0 : startChanList().set(0);
148 :
149 : // Temporary MS to get some info
150 0 : MeasurementSet ms(msname);
151 :
152 : // The relevant subtable names (there must be a better way...)
153 0 : sysPowTabName_ = ms.rwKeywordSet().asTable("SYSPOWER").tableName();
154 0 : calDevTabName_ = ms.rwKeywordSet().asTable("CALDEVICE").tableName();
155 :
156 : // Get spw total bandwidths
157 0 : MSColumns mscol(ms);
158 0 : const MSSpWindowColumns& spwcols = mscol.spectralWindow();
159 0 : effChBW_.resize(nSpw());
160 0 : for (Int ispw=0;ispw<nSpw();++ispw)
161 0 : effChBW_(ispw)=Vector<Double>(spwcols.effectiveBW()(0))(0);
162 :
163 0 : }
164 :
165 1 : EVLASwPow::EVLASwPow(const MSMetaInfoForCal& msmc) :
166 : VisCal(msmc), // virtual base
167 : VisMueller(msmc), // virtual base
168 : GJones(msmc), // immediate parent
169 1 : sysPowTabName_(""),
170 1 : calDevTabName_(""),
171 1 : correff_(Float(0.932)), // EVLA-specific net corr efficiency (4bit)
172 1 : frgrotscale_(Float(1.176)), // EVLA-specific fringe rotation mean _scale_
173 : //nyquist_(1.0),
174 2 : effChBW_()
175 :
176 : {
177 1 : if (prtlev()>2) cout << "EVLASwPow::EVLASwPow(msmc)" << endl;
178 :
179 1 : nChanParList().set(1);
180 1 : startChanList().set(0);
181 :
182 : // Temporary MS to get some info
183 1 : MeasurementSet ms(this->msmc().msname());
184 :
185 : // The relevant subtable names (there must be a better way...)
186 1 : sysPowTabName_ = ms.rwKeywordSet().asTable("SYSPOWER").tableName();
187 1 : calDevTabName_ = ms.rwKeywordSet().asTable("CALDEVICE").tableName();
188 :
189 : // Get spw total bandwidths
190 1 : MSColumns mscol(ms);
191 1 : const MSSpWindowColumns& spwcols = mscol.spectralWindow();
192 1 : effChBW_.resize(nSpw());
193 5 : for (Int ispw=0;ispw<nSpw();++ispw)
194 4 : effChBW_(ispw)=Vector<Double>(spwcols.effectiveBW()(0))(0);
195 :
196 1 : }
197 :
198 0 : EVLASwPow::EVLASwPow(const Int& nAnt) :
199 : VisCal(nAnt),
200 : VisMueller(nAnt),
201 0 : GJones(nAnt)
202 : {
203 0 : if (prtlev()>2) cout << "EVLASwPow::EVLASwPow(nAnt)" << endl;
204 :
205 0 : throw(AipsError("Cannot use EVLASwPow with generic ctor."));
206 :
207 0 : }
208 :
209 2 : EVLASwPow::~EVLASwPow() {
210 1 : if (prtlev()>2) cout << "EVLASwPow::~EVLASwPow()" << endl;
211 2 : }
212 :
213 1 : void EVLASwPow::setSpecify(const Record& specify) {
214 :
215 2 : LogMessage message(LogOrigin("EVLASwPow","setSpecify"));
216 :
217 : // Escape if SYSPOWER or CALDEVICE tables absent
218 1 : if (!Table::isReadable(sysPowTabName_))
219 0 : throw(AipsError("The SYSPOWER subtable is not present in the specified MS."));
220 1 : if (!Table::isReadable(calDevTabName_))
221 0 : throw(AipsError("The CALDEVICE subtable is not present in the specified MS."));
222 :
223 : // Not actually applying or solving
224 1 : setSolved(false);
225 1 : setApplied(false);
226 :
227 : // Collect Cal table parameters
228 1 : if (specify.isDefined("caltable")) {
229 1 : calTableName()=specify.asString("caltable");
230 :
231 1 : if (Table::isReadable(calTableName()))
232 : logSink() << "FYI: We are going to overwrite an existing CalTable: "
233 0 : << calTableName()
234 0 : << LogIO::POST;
235 : }
236 :
237 : // we are creating a table from scratch
238 1 : logSink() << "Creating " << typeName()
239 : // << " table from CANNED VALUES."
240 : << " table from MS SYSPOWER/CALDEVICE subtables."
241 1 : << LogIO::POST;
242 :
243 : // Create a caltable to fill up
244 1 : createMemCalTable();
245 :
246 : // Init the shapes of solveAllRPar, etc.
247 1 : initSolvePar();
248 :
249 1 : }
250 :
251 1 : void EVLASwPow::specify(const Record& specify) {
252 :
253 : // Escape if SYSPOWER or CALDEVICE tables absent
254 1 : if (!Table::isReadable(sysPowTabName_))
255 0 : throw(AipsError("The SYSPOWER subtable is not present in the specified MS."));
256 1 : if (!Table::isReadable(calDevTabName_))
257 0 : throw(AipsError("The CALDEVICE subtable is not present in the specified MS."));
258 :
259 : // Fill the Tcals first
260 1 : fillTcals();
261 :
262 : // Discern which kind of calibration to calculate
263 1 : SPType swptype(EVLASwPow::SWPOW);
264 1 : if (specify.isDefined("caltype")) {
265 1 : String ctype=specify.asString("caltype");
266 1 : swptype=sptype(ctype);
267 : //cout << "caltype=" << ctype << " " << swptype << " " << sptype(swptype) << endl;
268 1 : }
269 :
270 1 : logSink() << "Filling switched-power (" << sptype(swptype) << ") data from the SYSPOWER table." << LogIO::POST;
271 :
272 : // The net digital factor for antenna-based (voltage) gain
273 1 : Float dig=sqrt(correff_*frgrotscale_);
274 :
275 : // Keep a count of the number of Tsys found per spw/ant
276 1 : Matrix<Int> goodcount(nSpw(),nElem(),0), badcount(nSpw(),nElem(),0);
277 :
278 1 : Block<String> columns(2);
279 1 : columns[0] = "TIME";
280 1 : columns[1] = "SPECTRAL_WINDOW_ID";
281 1 : Table sysPowTab(sysPowTabName_,Table::Old);
282 1 : TableIterator sysPowIter(sysPowTab,columns);
283 :
284 : // Count iterations
285 1 : Int niter(0);
286 4907 : while (!sysPowIter.pastEnd()) {
287 4906 : ++niter;
288 4906 : sysPowIter.next();
289 : }
290 1 : sysPowIter.reset();
291 :
292 1 : logSink() << "Found " << niter << " TIME/SPW switched-power samples in SYSPOWER table" << LogIO::POST;
293 :
294 : // Iterate
295 : // Vectors for referencing slices of working info
296 1 : Vector<Float> currpsum,currpdif,currrq,currtcal,gain,tsys;
297 :
298 : // Emit progress meter reports (at least until we improve performance)
299 1 : cerr << "Switched-Power ("+sptype(swptype)+") calculation: 0";
300 2 : ProgressMeter pm(0.,niter , "", "", "", "", true, niter/100);
301 :
302 1 : Int iter(0);
303 4907 : while (!sysPowIter.pastEnd()) {
304 :
305 : // Update the progress meter
306 4906 : pm.update(iter);
307 :
308 4906 : Table itab(sysPowIter.table());
309 :
310 4906 : ScalarColumn<Int> spwCol(itab,"SPECTRAL_WINDOW_ID");
311 4906 : ScalarColumn<Double> timeCol(itab,"TIME");
312 4906 : ScalarColumn<Double> intervalCol(itab,"INTERVAL");
313 4906 : ScalarColumn<Int> antCol(itab,"ANTENNA_ID");
314 4906 : ArrayColumn<Float> swsumCol(itab,"SWITCHED_SUM");
315 4906 : ArrayColumn<Float> swdiffCol(itab,"SWITCHED_DIFF");
316 4906 : ArrayColumn<Float> rqCol(itab,"REQUANTIZER_GAIN");
317 :
318 4906 : Int ispw=spwCol(0);
319 :
320 4906 : Double timestamp=timeCol(0);
321 : // Double interval=intervalCol(0);
322 :
323 4906 : Vector<Int> ants;
324 4906 : antCol.getColumn(ants);
325 :
326 4906 : Matrix<Float> psum,pdif,rq;
327 4906 : swsumCol.getColumn(psum);
328 4906 : swdiffCol.getColumn(pdif);
329 4906 : rqCol.getColumn(rq);
330 4906 : IPosition psumshape(psum.shape());
331 4906 : IPosition pdifshape(pdif.shape());
332 :
333 4906 : AlwaysAssert(psumshape.isEqual(pdifshape),AipsError);
334 :
335 : // the number of receptors
336 4906 : Int nrec=psumshape(0);
337 :
338 : // Now prepare to record pars in the caltable
339 4906 : currSpw()=ispw;
340 4906 : refTime()=timestamp;
341 4906 : currField()=msmc().fieldIdAtTime(timestamp);
342 4906 : currScan()=msmc().scanNumberAtTime(timestamp);
343 :
344 : // Initialize solveAllRPar, etc.
345 4906 : solveAllRPar()=1.0;
346 4906 : solveAllParOK()=false;
347 4906 : solveAllParErr()=0.0; // what should we use here? ~1/bandwidth?
348 4906 : solveAllParSNR()=1.0;
349 :
350 4906 : IPosition blc(3,0,0,0),trc(3,2*nrec-1,0,0),stp(3,nrec,1,1);
351 24527 : for (uInt iant=0;iant<ants.nelements();++iant) {
352 19621 : Int thisant=ants(iant);
353 :
354 : // Reference this ant's info
355 19621 : currpsum.reference(psum.column(iant));
356 19621 : currpdif.reference(pdif.column(iant));
357 19621 : currrq.reference(rq.column(iant));
358 19621 : currtcal.reference(tcals_.xyPlane(ispw).column(thisant));
359 :
360 : // If any of the required values are goofy, we'll skip this antenna
361 : Bool good;
362 19621 : switch (swptype) {
363 0 : case EVLASwPow::SWPOW: {
364 0 : good=(allGT(currtcal,FLT_EPSILON) &&
365 0 : allGT(currpdif,FLT_EPSILON) &&
366 0 : allGT(currpsum,FLT_EPSILON));
367 : }
368 0 : case EVLASwPow::SWPOVERRQ:{
369 0 : good=(allGT(currtcal,FLT_EPSILON) &&
370 0 : allGT(currpdif,FLT_EPSILON) &&
371 0 : allGT(currpsum,FLT_EPSILON) &&
372 0 : allGT(currrq,FLT_EPSILON));
373 0 : break;
374 : }
375 0 : case EVLASwPow::RQ: {
376 0 : good=allGT(currrq,FLT_EPSILON);
377 0 : break;
378 : }
379 19621 : case EVLASwPow::SWPWTS: {
380 39242 : good=(allGT(currpsum,FLT_EPSILON) &&
381 19621 : allGT(currrq,FLT_EPSILON));
382 19621 : break;
383 : }
384 0 : default: {
385 0 : throw(AipsError("Unrecognized EVLA Switched Power type"));
386 : break;
387 : }
388 : }
389 19621 : blc(2)=trc(2)=thisant; // the MS ant index (not loop index)
390 :
391 19621 : blc(0)=0;
392 19621 : gain.reference(solveAllRPar()(blc,trc,stp).nonDegenerate(1)); // 'gain'
393 19621 : blc(0)=1;
394 19621 : tsys.reference(solveAllRPar()(blc,trc,stp).nonDegenerate(1)); // 'tsys'
395 :
396 19621 : if (!good) {
397 :
398 : // ensure transparent values
399 0 : gain=1.0;
400 0 : tsys=1.0;
401 0 : solveAllParOK().xyPlane(thisant)=false;
402 :
403 : // Increment bad counter
404 0 : ++badcount(ispw,thisant);
405 : }
406 : else {
407 :
408 : // Calculate "gain" and "tsys" for different modes
409 : // NB: gain includes correction for digital effects (loss and scale)
410 : // NB: No digital stuff in Tsys! net dig losses included OTF
411 : // in syncWtScale
412 :
413 19621 : switch (swptype) {
414 0 : case EVLASwPow::SWPOW: {
415 0 : gain=sqrt(currpdif/currtcal);
416 0 : gain*=dig; // scale by net digital factor
417 0 : tsys=(currtcal*currpsum/currpdif/2.0); // 'tsys'
418 0 : break;
419 : }
420 0 : case EVLASwPow::RQ: {
421 0 : gain=currrq; // RQ gain only!
422 0 : gain*=dig; // scale by net digital factor
423 0 : tsys=1.0;
424 0 : tsys/=square(currrq); // vis scale by rq req wt scale by 1/rq**2
425 0 : break;
426 : }
427 0 : case EVLASwPow::SWPOVERRQ: {
428 0 : gain=sqrt(currpdif/currtcal); // ordinary sw power gain
429 0 : gain/=currrq; // remove rq effect
430 0 : gain*=dig; // scale by net digital factor
431 0 : tsys=(currtcal*currpsum/currpdif/2.0); // 'tsys'
432 0 : break;
433 : }
434 19621 : case EVLASwPow::SWPWTS:{
435 19621 : gain=(currrq*dig); // include digital stuff, as for RQ
436 19621 : tsys = currpsum/2.0;
437 19621 : tsys/=square(currrq); // includ 1/rq**2, as for RQ
438 19621 : break;
439 : }
440 0 : default: {
441 0 : throw(AipsError("Unrecognized EVLA Switched Power type"));
442 : break;
443 : }
444 : }
445 19621 : solveAllParOK().xyPlane(thisant)=true;
446 :
447 : // Increment good counter
448 19621 : ++goodcount(ispw,thisant);
449 :
450 : }
451 :
452 : }
453 :
454 : // Record in the memory caltable
455 4906 : keepNCT();
456 :
457 4906 : sysPowIter.next();
458 4906 : ++iter;
459 :
460 4906 : }
461 :
462 : // Set scan and fieldid info
463 : // assignCTScanField(*ct_,msName());
464 :
465 1 : logSink() << "GOOD gain counts per spw for antenna Ids 0-"<<nElem()-1<<":" << LogIO::POST;
466 5 : for (Int ispw=0;ispw<nSpw();++ispw) {
467 4 : Vector<Int> goodcountspw(goodcount.row(ispw));
468 4 : if (sum(goodcountspw)>0)
469 : logSink() << " Spw " << ispw << ": " << goodcountspw
470 : << " (" << sum(goodcountspw) << ")"
471 4 : << LogIO::POST;
472 : else
473 0 : logSink() << "Spw " << ispw << ": NONE." << LogIO::POST;
474 4 : }
475 :
476 1 : logSink() << "BAD gain counts per spw for antenna Ids 0-"<<nElem()-1<<":" << LogIO::POST;
477 5 : for (Int ispw=0;ispw<nSpw();++ispw) {
478 4 : Vector<Int> badcountspw(badcount.row(ispw));
479 4 : if (sum(badcountspw)>0)
480 : logSink() << " Spw " << ispw << ": " << badcountspw
481 : << " (" << sum(badcountspw) << ")"
482 0 : << LogIO::POST;
483 4 : }
484 :
485 1 : logSink() << "BAD gain count FRACTION per spw for antenna Ids 0-"<<nElem()-1<<":" << LogIO::POST;
486 5 : for (Int ispw=0;ispw<nSpw();++ispw) {
487 4 : Vector<Float> badcountspw(badcount.row(ispw).shape());
488 4 : Vector<Float> goodcountspw(goodcount.row(ispw).shape());
489 4 : convertArray(badcountspw,badcount.row(ispw));
490 4 : convertArray(goodcountspw,goodcount.row(ispw));
491 4 : if (sum(badcountspw)>0.0f) {
492 0 : Vector<Float> fracbad=badcountspw/(badcountspw+goodcountspw);
493 0 : fracbad=floor(1000.0f*fracbad)/1000.0f;
494 0 : Float fracbadsum=sum(badcountspw)/(sum(badcountspw)+sum(goodcountspw));
495 0 : fracbadsum=floor(1000.0f*fracbadsum)/1000.0f;
496 : logSink() << " Spw " << ispw << ": " << fracbad
497 : << " (" << fracbadsum << ")"
498 0 : << LogIO::POST;
499 0 : }
500 4 : }
501 :
502 :
503 1 : }
504 :
505 1 : void EVLASwPow::fillTcals() {
506 :
507 1 : logSink() << "Filling Tcals from the CALDEVICE table." << LogIO::POST;
508 :
509 1 : Block<String> columns(2);
510 1 : columns[0] = "SPECTRAL_WINDOW_ID";
511 1 : columns[1] = "ANTENNA_ID";
512 1 : Table calDevTab(calDevTabName_,Table::Old);
513 1 : TableIterator calDevIter(calDevTab,columns);
514 :
515 1 : tcals_.resize(2,nElem(),nSpw());
516 1 : tcals_.set(-1.0f);
517 :
518 : // Iterate over CALDEVICE table
519 1 : Int iter(0);
520 1 : Vector<Int> islot(nSpw(),0);
521 :
522 1 : Bool ignoreSolar(false);
523 17 : while (!calDevIter.pastEnd()) {
524 :
525 16 : Table itab(calDevIter.table());
526 :
527 16 : ScalarColumn<Int> spwCol(itab,"SPECTRAL_WINDOW_ID");
528 16 : ScalarColumn<Double> timeCol(itab,"TIME");
529 16 : ScalarColumn<Double> intervalCol(itab,"INTERVAL");
530 16 : ScalarColumn<Int> antCol(itab,"ANTENNA_ID");
531 16 : ScalarColumn<Int> numLoadCol(itab,"NUM_CAL_LOAD");
532 :
533 16 : ArrayColumn<Float> noiseCalCol(itab,"NOISE_CAL");
534 :
535 16 : Int ispw=spwCol(0);
536 16 : Int iant=antCol(0);
537 16 : Int nTcal=noiseCalCol(0).nelements();
538 :
539 16 : if (nTcal==1) {
540 : // One value (not clear this was ever the case)
541 0 : Vector<Float> thisTcal=noiseCalCol(0);
542 0 : AlwaysAssert(thisTcal.nelements()==1,AipsError);
543 0 : tcals_.xyPlane(ispw).column(iant)=thisTcal(0);
544 0 : }
545 16 : else if (nTcal==2) {
546 : // Pre-solar Tcal support: two values
547 16 : Vector<Float> thisTcal=noiseCalCol(0);
548 16 : AlwaysAssert(thisTcal.nelements()==2,AipsError);
549 16 : tcals_.xyPlane(ispw).column(iant)=thisTcal;
550 16 : }
551 0 : else if (nTcal==4) {
552 0 : ignoreSolar=true; // we will ignore solar Tcals
553 0 : Matrix<Float> thisTcalMat=noiseCalCol(0);
554 0 : AlwaysAssert(thisTcalMat.shape()==IPosition(2,2,2),AipsError);
555 0 : Int tcalset(0); // first pair, for now, which is ordinary non-solar Tcals
556 0 : tcals_.xyPlane(ispw).column(iant)=thisTcalMat(Slice(tcalset,1,1),Slice());
557 0 : }
558 :
559 : // Increment the iterator
560 16 : ++iter;
561 16 : calDevIter.next();
562 16 : }
563 :
564 : // Report if we ignored solar Tcals
565 1 : if (ignoreSolar)
566 0 : logSink() << "Ignoring the SOLAR Tcals, which seem to be present in CALDEVICE." << LogIO::POST;
567 :
568 :
569 1 : }
570 :
571 0 : void EVLASwPow::calcAllJones() {
572 :
573 : // 0th and 2nd pars are the gains
574 : // currJElem()=currRPar()(Slice(0,2,2),Slice(),Slice()); // NEWCALTABLE
575 0 : convertArray(currJElem(),currRPar()(Slice(0,2,2),Slice(),Slice()));
576 0 : currJElemOK()=currParOK()(Slice(0,2,2),Slice(),Slice());
577 :
578 0 : }
579 :
580 0 : void EVLASwPow::syncWtScale() {
581 :
582 0 : Int nPolWt=currRPar().shape()(0)/2;
583 0 : currWtScale().resize(nPolWt,1,nAnt());
584 :
585 0 : Cube<Float> Tsys(currRPar()(Slice(1,2,2),Slice(),Slice()));
586 0 : Tsys(Tsys<FLT_MIN)=1.0; // OK?
587 :
588 0 : currWtScale() = 1.0f/Tsys;
589 0 : currWtScale()*=correff_; // reduce by correlator efficiency (per ant)
590 :
591 : // cout << "Tsys = " << Tsys << endl;
592 : // cout << "currWtScale() = " << currWtScale() << endl;
593 :
594 0 : }
595 :
596 : /*
597 : void EVLASwPow::updateWt(Vector<Float>& wt,const Int& a1,const Int& a2) {
598 :
599 : Vector<Float> ws1(currWtScale().column(a1));
600 : Vector<Float> ws2(currWtScale().column(a2));
601 :
602 : if (a1==0 && a2==1) {
603 : cout << "wt = " << wt << endl;
604 : cout << "ws1 = " << ws1 << endl;
605 : cout << "ws2 = " << ws2 << endl;
606 : }
607 :
608 : VisJones::updateWt(wt,a1,a2);
609 :
610 : if (a1==0 && a2==1) {
611 : cout << "wt = " << wt << endl;
612 : }
613 : }
614 : */
615 :
616 : } //# NAMESPACE CASA - END
|
