Line data Source code
1 : //# StandardVisCal.cc: Implementation of Standard VisCal types
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/EJones.h>
28 :
29 : #include <synthesis/MeasurementComponents/MSMetaInfoForCal.h>
30 :
31 : #include <msvis/MSVis/VisBuffer.h>
32 : #include <msvis/MSVis/VisBuffAccumulator.h>
33 : #include <casacore/ms/MeasurementSets/MSColumns.h>
34 : #include <synthesis/MeasurementEquations/VisEquation.h>
35 :
36 : #include <casacore/tables/Tables/Table.h>
37 : #include <casacore/tables/Tables/TableIter.h>
38 : #include <casacore/tables/Tables/TableUtil.h>
39 : #include <casacore/tables/TaQL/ExprNode.h>
40 :
41 : #include <casacore/casa/Arrays/ArrayMath.h>
42 : #include <casacore/casa/BasicSL/String.h>
43 : #include <casacore/casa/Utilities/Assert.h>
44 : #include <casacore/casa/Exceptions/Error.h>
45 : #include <casacore/casa/OS/Memory.h>
46 : #include <casacore/casa/System/Aipsrc.h>
47 : #include <casacore/scimath/Functionals/ScalarSampledFunctional.h>
48 : #include <casacore/scimath/Functionals/Interpolate1D.h>
49 : #include <casacore/scimath/Mathematics/Combinatorics.h>
50 :
51 : #include <sstream>
52 :
53 : #include <casacore/casa/Logging/LogMessage.h>
54 : #include <casacore/casa/Logging/LogSink.h>
55 :
56 : using namespace casacore;
57 : namespace casa { //# NAMESPACE CASA - BEGIN
58 :
59 :
60 : // **********************************************************
61 : // EGainCurve
62 : //
63 :
64 0 : EGainCurve::EGainCurve(VisSet& vs) :
65 : VisCal(vs),
66 : VisMueller(vs),
67 : SolvableVisJones(vs),
68 : za_(),
69 0 : eff_(nSpw(),1.0)
70 : {
71 0 : if (prtlev()>2) cout << "EGainCurve::EGainCurve(vs)" << endl;
72 0 : }
73 :
74 0 : EGainCurve::EGainCurve(String msname,Int MSnAnt,Int MSnSpw) :
75 : VisCal(msname,MSnAnt,MSnSpw),
76 : VisMueller(msname,MSnAnt,MSnSpw),
77 : SolvableVisJones(msname,MSnAnt,MSnSpw),
78 : za_(),
79 0 : eff_(nSpw(),1.0)
80 : {
81 0 : if (prtlev()>2) cout << "EGainCurve::EGainCurve(msname,MSnAnt,MSnSpw)" << endl;
82 0 : }
83 :
84 3 : EGainCurve::EGainCurve(const MSMetaInfoForCal& msmc) :
85 : VisCal(msmc),
86 : VisMueller(msmc),
87 : SolvableVisJones(msmc),
88 : za_(),
89 3 : eff_(nSpw(),1.0)
90 : {
91 3 : if (prtlev()>2) cout << "EGainCurve::EGainCurve(msmc)" << endl;
92 3 : }
93 :
94 :
95 :
96 4 : EGainCurve::~EGainCurve() {
97 3 : if (prtlev()>2) cout << "EGainCurve::~EGainCurve()" << endl;
98 4 : }
99 :
100 1 : void EGainCurve::setApply(const Record& applypar) {
101 :
102 2 : LogMessage message;
103 :
104 : // Extract user's table name
105 2 : String usertab("");
106 1 : if (applypar.isDefined("table")) {
107 1 : usertab=applypar.asString("table");
108 : }
109 :
110 1 : if (usertab=="") {
111 :
112 0 : { ostringstream o;
113 0 : o<< "Invoking gaincurve application without a caltable (e.g., " << endl
114 0 : << " via gaincurve=T in calibration tasks) will soon be disabled." << endl
115 0 : << " Please begin using gencal to generate a gaincurve caltable, " << endl
116 0 : << " and then supply that table in the standard manner." << endl;
117 0 : message.message(o);
118 0 : message.priority(LogMessage::WARN);
119 0 : logSink().post(message);
120 : }
121 :
122 0 : String tempname="temporary.gaincurve";
123 :
124 : // Generate automatically via specify mechanism
125 0 : Record specpar;
126 0 : specpar.define("caltable",tempname);
127 0 : specpar.define("caltype","gc");
128 0 : setSpecify(specpar);
129 0 : specify(specpar);
130 0 : storeNCT();
131 0 : delete ct_; ct_=NULL; // so we can form it in the standard way...
132 :
133 : // Load the caltable for standard apply
134 0 : Record newapplypar=applypar;
135 0 : newapplypar.define("table",tempname);
136 0 : SolvableVisCal::setApply(newapplypar);
137 :
138 : // Delete the temporary gaincurve disk table (in-mem copy still ok)
139 0 : if (calTableName()==tempname &&
140 0 : TableUtil::canDeleteTable(calTableName()) ) {
141 0 : TableUtil::deleteTable(calTableName());
142 : }
143 :
144 : // Revise name that will appear in log messages, etc.
145 0 : calTableName()="<"+tempname+">";
146 :
147 : }
148 : else
149 : // Standard apply via table
150 1 : SolvableVisCal::setApply(applypar);
151 :
152 : // Resize za()
153 1 : za().resize(nAnt());
154 :
155 1 : }
156 :
157 0 : void EGainCurve::setCallib(const Record& applypar,const MeasurementSet& selms) {
158 :
159 0 : LogMessage message;
160 :
161 : // Standard setCallib
162 0 : SolvableVisCal::setCallib(applypar,selms);
163 :
164 : // Resize za()
165 0 : za().resize(nAnt());
166 :
167 0 : }
168 :
169 1 : void EGainCurve::setSpecify(const Record& specify) {
170 :
171 : // Get some information from the MS to help us find
172 : // the right gain curves
173 :
174 2 : MeasurementSet ms(msName());
175 2 : MSColumns mscol(ms);
176 :
177 : // The antenna names
178 1 : const MSAntennaColumns& antcol(mscol.antenna());
179 1 : antnames_ = antcol.name().getColumn();
180 :
181 : // Observation info
182 1 : const MSObservationColumns& obscol(mscol.observation());
183 :
184 2 : String telescope(obscol.telescopeName()(0));
185 :
186 : // Parse infile
187 1 : if (specify.isDefined("infile"))
188 1 : gainCurveSrc_=specify.asString("infile");
189 :
190 : // Use external file, if specified
191 1 : if (gainCurveSrc_!="") {
192 0 : if ( !Table::isReadable(gainCurveSrc_) )
193 0 : throw(AipsError("Specified gain curve table: "+gainCurveSrc_+" is unreadable or absent."));
194 :
195 : // Specified table exists, so we'll try to use it
196 : logSink() << "Using user-specified gaincurve table: "
197 0 : << gainCurveSrc_
198 0 : << LogIO::POST;
199 :
200 : }
201 : // If VLA, use standard file
202 1 : else if (telescope.contains("VLA")) {
203 0 : gainCurveSrc_=Aipsrc::aipsRoot() + "/data/nrao/VLA/GainCurves";
204 0 : if ( !Table::isReadable(gainCurveSrc_) )
205 0 : throw(AipsError("Standard VLA gain curve table "+gainCurveSrc_+" is unreadable or absent."));
206 :
207 : // VLA gaincurve exists, so we'll try to use it
208 : logSink() << "Using standard VLA gaincurve table from the data repository: "
209 0 : << gainCurveSrc_
210 0 : << LogIO::POST;
211 :
212 : // Strip any ea/va baloney from the MS antenna names so
213 : // they match the integers (as strings) in the GainCurves table
214 0 : for (uInt iant=0;iant<antnames_.nelements();++iant) {
215 0 : antnames_(iant)=antnames_(iant).from(RXint);
216 0 : if (antnames_(iant).find('0')==0)
217 0 : antnames_(iant)=antnames_(iant).after('0');
218 : }
219 : }
220 : // Unspecified and not VLA: gain curve unsupported...
221 : else {
222 1 : throw(AipsError("Automatic gain curve not supported for "+telescope));
223 : }
224 :
225 0 : Vector<Double> timerange(obscol.timeRange()(0));
226 0 : obstime_ = timerange(0);
227 :
228 0 : const MSSpWindowColumns& spwcol(mscol.spectralWindow());
229 0 : spwfreqs_.resize(nSpw());
230 0 : spwfreqs_.set(0.0);
231 0 : spwbands_.resize(nSpw());
232 0 : spwbands_.set(String("?"));
233 0 : for (Int ispw=0;ispw<nSpw();++ispw) {
234 0 : Vector<Double> chanfreqs=spwcol.chanFreq()(ispw);
235 0 : spwfreqs_(ispw)=chanfreqs(chanfreqs.nelements()/2);
236 0 : String bname=spwcol.name()(ispw);
237 0 : if (bname.contains("EVLA"))
238 0 : spwbands_(ispw)=String(bname.before("#")).after("EVLA_");
239 : }
240 : // cout << "spwfreqs_ = " << spwfreqs_ << endl;
241 : // cout << "spwbands_ = " << spwbands_ << endl;
242 :
243 : // Neither applying nor solving in specify context
244 0 : setSolved(false);
245 0 : setApplied(false);
246 :
247 : // Collect Cal table parameters
248 0 : if (specify.isDefined("caltable")) {
249 0 : calTableName()=specify.asString("caltable");
250 :
251 0 : if (Table::isReadable(calTableName()))
252 : logSink() << "FYI: We are going to overwrite an existing CalTable: "
253 0 : << calTableName()
254 0 : << LogIO::POST;
255 : }
256 :
257 : // Create a new caltable to fill
258 0 : createMemCalTable();
259 :
260 : // Setup shape of solveAllRPar
261 0 : initSolvePar();
262 :
263 : /* From AIPS TYAPL (2012Sep14):
264 :
265 : DATA TF / 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8,
266 : * 1.9, 2.0, 2.0, 2.3, 2.7, 3.0, 3.5, 3.7, 4.0, 4.0,
267 : * 5.0, 6.0, 7.0, 8.0, 8.0, 12.0, 12.0, 13.0, 14.0, 15.0,
268 : * 16.0, 17.0, 18.0, 19.0, 24.0, 26.0, 26.5, 28.0, 33.0, 38.0,
269 : * 40.0, 43.0, 48.0/
270 : C Rick Perley efficiencies
271 : DATA TE /0.45, 0.48, 0.48, 0.45, 0.46, 0.45, 0.43, 0.44, 0.44,
272 : * 0.49, 0.48, 0.52, 0.52, 0.51, 0.53, 0.55, 0.53, 0.54, 0.55,
273 : * 0.54, 0.56, 0.62, 0.64, 0.60, 0.60, 0.65, 0.65, 0.62, 0.58,
274 : * 0.59, 0.60, 0.60, 0.57, 0.52, 0.48, 0.50, 0.49, 0.42, 0.35,
275 : * 0.29, 0.28, 0.26/
276 : */
277 :
278 :
279 0 : Double Efffrq[] = {1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8,
280 : 1.9, 2.0, 2.0+1e-9, 2.3, 2.7, 3.0, 3.5, 3.7, 4.0, 4.0+1e-9,
281 : 5.0, 6.0, 7.0, 8.0, 8.0+1e-9, 12.0, 12.0+1e-9, 13.0, 14.0, 15.0,
282 : 16.0, 17.0, 18.0, 19.0, 24.0, 26.0, 26.5, 28.0, 33.0, 38.0,
283 : 40.0, 43.0, 48.0};
284 0 : Double Eff[] = {0.45, 0.48, 0.48, 0.45, 0.46, 0.45, 0.43, 0.44, 0.44,
285 : 0.49, 0.48, 0.52, 0.52, 0.51, 0.53, 0.55, 0.53, 0.54, 0.55,
286 : 0.54, 0.56, 0.62, 0.64, 0.60, 0.60, 0.65, 0.65, 0.62, 0.58,
287 : 0.59, 0.60, 0.60, 0.57, 0.52, 0.48, 0.50, 0.49, 0.42, 0.35,
288 : 0.29, 0.28, 0.26};
289 :
290 0 : Vector<Double> f,ef;
291 :
292 0 : f.takeStorage(IPosition(1,42),Efffrq);
293 0 : ef.takeStorage(IPosition(1,42),Eff);
294 :
295 : // Fractional -> K/Jy (for 25m)
296 0 : ef/=Double(5.622);
297 :
298 : // convert to voltagey units
299 0 : ef=sqrt(ef);
300 :
301 0 : ScalarSampledFunctional<Double> fx(f);
302 0 : ScalarSampledFunctional<Double> fy(ef);
303 0 : Interpolate1D<Double,Double> eff(fx,fy);
304 0 : eff.setMethod(Interpolate1D<Double,Double>::linear);
305 0 : for (Int ispw=0;ispw<nSpw();++ispw) {
306 0 : Double fghz=spwfreqs_(ispw)/1.e9;
307 0 : eff_(ispw)=eff(fghz);
308 : }
309 :
310 : // cout << "spwfreqs_ = " << spwfreqs_ << endl;
311 : // cout << "eff_ = " << eff_ << endl;
312 :
313 :
314 0 : }
315 :
316 0 : void EGainCurve::specify(const Record& specify) {
317 :
318 0 : LogMessage message;
319 :
320 0 : Bool doeff(false);
321 0 : Bool dogc(false);
322 0 : if (specify.isDefined("caltype")) {
323 0 : String caltype=specify.asString("caltype");
324 : //cout << "caltype=" << caltype << endl;
325 0 : doeff = (caltype.contains("eff"));
326 0 : dogc = (caltype.contains("gc"));
327 : }
328 :
329 : // Open raw gain curve source table
330 0 : Table rawgaintab(gainCurveSrc_);
331 :
332 0 : logSink() << "Using " << Path(gainCurveSrc_).absoluteName()
333 0 : << " nrow=" << rawgaintab.nrow() << LogIO::POST;
334 :
335 : // Select on Time
336 0 : Table timegaintab = rawgaintab(rawgaintab.col("BTIME") <= obstime_
337 0 : && rawgaintab.col("ETIME") > obstime_);
338 :
339 : // ...for each spw:
340 0 : Vector<Bool> spwFound_(nSpw(),false);
341 0 : spwFound_.set(false); // will set true when we find them
342 :
343 0 : for (Int ispw=0; ispw<nSpw(); ispw++) {
344 :
345 0 : currSpw()=ispw;
346 :
347 :
348 0 : if (dogc) {
349 :
350 : // Select on freq:
351 0 : Table freqgaintab=timegaintab(timegaintab.col("BANDNAME")==spwbands_(ispw));
352 :
353 : // If we fail to select anything by bandname, try to select by freq
354 : // (this can get wrong answer near band edges if center freq "out-of-band")
355 0 : if (freqgaintab.nrow()==0)
356 0 : freqgaintab = timegaintab(timegaintab.col("BFREQ") <= spwfreqs_(ispw)
357 0 : && timegaintab.col("EFREQ") > spwfreqs_(ispw));
358 :
359 0 : if (freqgaintab.nrow() > 0) {
360 :
361 : /*
362 : { logSink() // << "EGainCurve::specify:"
363 : << " Found the following gain curve coefficient data" << endl
364 : << " for spectral window = " << ispw << " (band=" << spwbands_(ispw) << ", center freq="
365 : << spwfreqs_(ispw) << "):" << LogIO::POST;
366 : }
367 : */
368 : // Find nominal gain curve
369 : // Nominal antenna is named "0" (this is a VLA convention)
370 0 : Matrix<Float> nomgain(4,2,0.0);
371 : {
372 0 : Table nomgaintab = freqgaintab(freqgaintab.col("ANTENNA")=="0");
373 0 : if (nomgaintab.nrow() > 0) {
374 0 : ArrayColumn<Float> gain(nomgaintab,"GAIN");
375 0 : nomgain=gain(0);
376 : } else {
377 : // nominal gain is unity
378 0 : nomgain(0,0)=1.0;
379 0 : nomgain(0,1)=1.0;
380 : }
381 : }
382 :
383 0 : solveAllParOK()=true;
384 :
385 0 : ArrayIterator<Float> piter(solveAllRPar(),1);
386 :
387 0 : for (Int iant=0; iant<nAnt(); iant++,piter.next()) {
388 :
389 : // Select antenna by name
390 0 : Table antgaintab = freqgaintab(freqgaintab.col("ANTENNA")==antnames_(iant));
391 0 : if (antgaintab.nrow() > 0) {
392 0 : ArrayColumn<Float> gain(antgaintab,"GAIN");
393 0 : piter.array().nonDegenerate().reform(gain(0).shape())=gain(0);
394 : } else {
395 0 : piter.array().nonDegenerate().reform(nomgain.shape())=nomgain;
396 : }
397 :
398 : /*
399 : {
400 : logSink() << " Antenna=" << antnames_(iant) << ": "
401 : << piter.array().nonDegenerate() << LogIO::POST;
402 : }
403 : */
404 : }
405 :
406 0 : spwFound_(currSpw())=true;
407 :
408 : }
409 : else {
410 : logSink() << "Could not find gain curve data for Spw="
411 0 : << ispw << " (reffreq=" << spwfreqs_(ispw)/1.0e9 << " GHz) "
412 0 : << "Using flat unit gaincurve." << LogIO::POST;
413 : // We used to punt here
414 : //throw(AipsError(o.str()));
415 :
416 : // Use unity
417 0 : solveAllParOK()=true;
418 0 : solveAllRPar().set(0.0);
419 0 : solveAllRPar()(Slice(0,1,1),Slice(),Slice()).set(1.0);
420 0 : solveAllRPar()(Slice(4,1,1),Slice(),Slice()).set(1.0);
421 :
422 : }
423 : }
424 : else {
425 : // Use unity, flat
426 0 : solveAllParOK()=true;
427 0 : solveAllRPar().set(0.0);
428 0 : solveAllRPar()(Slice(0,1,1),Slice(),Slice()).set(1.0);
429 0 : solveAllRPar()(Slice(4,1,1),Slice(),Slice()).set(1.0);
430 0 : spwFound_(currSpw())=true;
431 : }
432 :
433 : // Scale by efficiency factor, if requested
434 0 : if (doeff) {
435 0 : solveAllRPar()*=Float(eff_(ispw));
436 : }
437 :
438 : // Record in the memory caltable
439 0 : keepNCT();
440 :
441 : } // ispw
442 :
443 0 : if (allEQ(spwFound_,false))
444 0 : throw(AipsError("Found no gaincurve data for any spw."));
445 :
446 :
447 : // Reset currSpw()
448 0 : currSpw()=0;
449 :
450 : // Resize za()
451 0 : za().resize(nAnt());
452 :
453 0 : }
454 :
455 :
456 0 : void EGainCurve::guessPar(VisBuffer&) {
457 :
458 0 : throw(AipsError("Spurious attempt to guess EGainCurve for solving!"));
459 :
460 : }
461 :
462 : // EGainCurve needs to zenith angle (old VB version)
463 0 : void EGainCurve::syncMeta(const VisBuffer& vb) {
464 :
465 : // Call parent (sets currTime())
466 0 : SolvableVisJones::syncMeta(vb);
467 :
468 : // Current time's zenith angle... (in _degrees_)
469 0 : za().resize(nAnt());
470 0 : Vector<MDirection> antazel(vb.azel(currTime()));
471 0 : Double* a=za().data();
472 0 : for (Int iant=0;iant<nAnt();++iant,++a)
473 0 : (*a)=90.0 - antazel(iant).getAngle().getValue()(1)*180.0/C::pi;
474 :
475 0 : }
476 :
477 : // EGainCurve needs to zenith angle (VB2 version)
478 628 : void EGainCurve::syncMeta2(const vi::VisBuffer2& vb) {
479 :
480 : // Call parent (sets currTime())
481 628 : SolvableVisJones::syncMeta2(vb);
482 :
483 : // Current time's zenith angle...(in _degrees_)
484 628 : za().resize(nAnt());
485 1256 : Vector<MDirection> antazel(vb.azel(currTime()));
486 628 : Double* a=za().data();
487 6908 : for (Int iant=0;iant<nAnt();++iant,++a)
488 6280 : (*a)=90.0 - antazel(iant).getAngle().getValue()(1)*180.0/C::pi;
489 :
490 628 : }
491 :
492 :
493 :
494 628 : void EGainCurve::calcPar() {
495 :
496 628 : if (prtlev()>6) cout << " EGainCurve::calcPar()" << endl;
497 :
498 : // Could do the following in setApply, so it only happens once?
499 628 : if (ci_ || cpp_)
500 628 : SolvableVisCal::calcPar();
501 : else
502 0 : throw(AipsError("Problem in EGainCurve::calcPar()"));
503 :
504 : // Pars now valid, matrices not yet
505 628 : validateP();
506 628 : invalidateJ();
507 :
508 628 : }
509 :
510 0 : void EGainCurve::calcAllJones() {
511 :
512 0 : if (prtlev()>6) cout << " EGainCurve::calcAllJones()" << endl;
513 :
514 : // Nominally no gain curve effect
515 0 : currJElem()=Complex(1.0);
516 0 : currJElemOK()=false;
517 :
518 : /*
519 : cout << "currSpw() = " << currSpw() << endl;
520 : cout << " spwMap() = " << spwMap() << endl;
521 : cout << " currRPar().shape() = " << currRPar().shape() << endl;
522 : cout << " currRPar() = " << currRPar() << endl;
523 : */
524 :
525 0 : Complex* J=currJElem().data();
526 0 : Bool* JOk=currJElemOK().data();
527 0 : Float* c=currRPar().data();
528 0 : Double* a=za().data();
529 :
530 : Double loss, ang;
531 0 : for (Int iant=0; iant<nAnt(); ++iant,++a)
532 0 : for (Int ipol=0;ipol<2;++ipol,++J,++JOk) {
533 0 : loss=Double(*c);
534 0 : ++c;
535 0 : ang=1.0;
536 0 : for (Int ipar=1;ipar<4;++ipar,++c) {
537 0 : ang*=(*a);
538 0 : loss+=((*c)*ang);
539 : }
540 0 : (*J) = Complex(loss);
541 0 : (*JOk) = true;
542 : }
543 :
544 0 : }
545 :
546 :
547 0 : EPowerCurve::EPowerCurve(VisSet& vs) :
548 : VisCal(vs),
549 : VisMueller(vs),
550 : EGainCurve(vs),
551 0 : gainCurveTabName_("")
552 : {
553 0 : if (prtlev()>2) cout << "EPowerCurve::EPowerCurve(vs)" << endl;
554 0 : }
555 :
556 0 : EPowerCurve::EPowerCurve(String msname,Int MSnAnt,Int MSnSpw) :
557 : VisCal(msname,MSnAnt,MSnSpw),
558 : VisMueller(msname,MSnAnt,MSnSpw),
559 : EGainCurve(msname,MSnAnt,MSnSpw),
560 0 : gainCurveTabName_("")
561 : {
562 0 : if (prtlev()>2) cout << "EPowerCurve::EPowerCurve(msname,MSnAnt,MSnSpw)" << endl;
563 :
564 : // Temporary MS to get some info
565 0 : MeasurementSet ms(msname);
566 :
567 : // The relevant subtable names (there must be a better way...)
568 0 : gainCurveTabName_ = ms.rwKeywordSet().asTable("GAIN_CURVE").tableName();
569 0 : }
570 :
571 2 : EPowerCurve::EPowerCurve(const MSMetaInfoForCal& msmc) :
572 : VisCal(msmc),
573 : VisMueller(msmc),
574 : EGainCurve(msmc),
575 2 : gainCurveTabName_("")
576 : {
577 2 : if (prtlev()>2) cout << "EPowerCurve::EPowerCurve(msmc)" << endl;
578 :
579 : // Temporary MS to get some info
580 2 : MeasurementSet ms(this->msmc().msname());
581 :
582 : // The relevant subtable names (there must be a better way...)
583 2 : gainCurveTabName_ = ms.rwKeywordSet().asTable("GAIN_CURVE").tableName();
584 2 : }
585 :
586 4 : EPowerCurve::~EPowerCurve() {
587 2 : if (prtlev()>2) cout << "EPowerCurve::~EPowerCurve()" << endl;
588 4 : }
589 :
590 1 : void EPowerCurve::setSpecify(const Record& specify) {
591 :
592 : // Neither applying nor solving in specify context
593 1 : setSolved(false);
594 1 : setApplied(false);
595 :
596 : // Collect Cal table parameters
597 1 : if (specify.isDefined("caltable")) {
598 1 : calTableName()=specify.asString("caltable");
599 :
600 1 : if (Table::isReadable(calTableName()))
601 : logSink() << "FYI: We are going to overwrite an existing CalTable: "
602 0 : << calTableName()
603 0 : << LogIO::POST;
604 : }
605 :
606 : // Create a new caltable to fill
607 1 : createMemCalTable();
608 :
609 : // Setup shape of solveAllRPar
610 1 : initSolvePar();
611 1 : }
612 :
613 1 : void EPowerCurve::specify(const Record& specify) {
614 :
615 : // Escape if GAIN_CURVE table absent
616 1 : if (!Table::isReadable(gainCurveTabName()))
617 0 : throw(AipsError("The GAIN_CURVE subtable is not present in the specified MS. Gain curves unavailable."));
618 :
619 : // Construct matrix for EL to ZA conversion of polynomials.
620 2 : Matrix<Float> m_el(nPar(), nPar(), 0.0);
621 9 : for (Int i = 0; i < nPar()/2; i++) {
622 72 : for (Int j = 0; j < nPar()/2; j++) {
623 64 : if (i > j)
624 28 : continue;
625 36 : m_el(i, j) = Combinatorics::choose(j, i) *
626 36 : pow(-1.0, j) * pow(-90.0, (j - i));
627 36 : m_el(nPar()/2 + i, nPar()/2 + j) = m_el(i, j);
628 : }
629 : }
630 :
631 2 : Table gainCurveTab(gainCurveTabName(),Table::Old);
632 :
633 3 : for (Int ispw=0; ispw<nSpw(); ispw++) {
634 6 : Table itab = gainCurveTab(gainCurveTab.col("SPECTRAL_WINDOW_ID")==ispw);
635 :
636 4 : ScalarColumn<Double> timeCol(itab, "TIME");
637 4 : ScalarColumn<Double> intervalCol(itab, "INTERVAL");
638 4 : ScalarColumn<Int> antCol(itab,"ANTENNA_ID");
639 4 : ScalarColumn<Int> spwCol(itab,"SPECTRAL_WINDOW_ID");
640 4 : ScalarColumn<String> typeCol(itab,"TYPE");
641 4 : ScalarColumn<Int> numPolyCol(itab, "NUM_POLY");
642 4 : ArrayColumn<Float> gainCol(itab, "GAIN");
643 2 : ArrayColumn<Float> sensCol(itab, "SENSITIVITY");
644 :
645 22 : for (Int irow=0; irow<itab.nrow(); irow++) {
646 20 : Int iant=antCol(irow);
647 20 : currSpw()=ispw;
648 :
649 40 : Matrix<Float> m;
650 20 : if (typeCol(irow) == "POWER(ZA)" || typeCol(irow) == "VOLTAGE(ZA)")
651 20 : m = Matrix<Float>::identity(nPar());
652 0 : else if (typeCol(irow) == "POWER(EL)" || typeCol(irow) == "VOLTAGE(EL)")
653 0 : m = m_el;
654 : else
655 0 : throw(AipsError("The " + typeCol(irow) + "gain curve type is not supported. Gain curves unavailable."));
656 :
657 : // Initialize solveAllParOK, etc.
658 20 : solveAllParOK()=true; // Assume all ok
659 20 : solveAllParErr()=0.0; // what should we use here?
660 20 : solveAllParSNR()=1.0;
661 :
662 20 : Double begin = timeCol(irow) - intervalCol(irow) / 2;
663 20 : Double end = timeCol(irow) + intervalCol(irow) / 2;
664 :
665 : // Warn if we need to truncate the polynomial?
666 20 : Int npoly = min(numPolyCol(irow), nPar()/2);
667 :
668 40 : Vector<Float> gain(nPar(), 0.0);
669 80 : for (Int i = 0; i < npoly; i++) {
670 60 : gain(i) = gainCol(irow)(IPosition(2, 0, i));
671 60 : gain(nPar()/2 + i) = gainCol(irow)(IPosition(2, 1, i));
672 : }
673 :
674 : // Convert to ZA polynomial
675 20 : gain = product(m, gain);
676 :
677 : // Square voltage to get power.
678 20 : if (typeCol(irow).startsWith("VOLTAGE")) {
679 0 : Vector<Float> v(nPar(), 0.0);
680 0 : for (Int i = 0; i < nPar()/2; i++) {
681 0 : for (Int j = 0; j < nPar()/2; j++) {
682 0 : if (i + j < nPar()/2)
683 0 : v(i + j) += gain(i) * gain(j);
684 : }
685 : }
686 0 : gain = v;
687 : }
688 :
689 180 : for (Int i = 0; i < nPar()/2; i++) {
690 160 : gain(i) *= sensCol(irow)(IPosition(1, 0));
691 160 : gain(nPar()/2 + i) *= sensCol(irow)(IPosition(1, 1));
692 : }
693 :
694 20 : ct_->addRow(1);
695 :
696 20 : CTMainColumns ncmc(*ct_);
697 :
698 : // We are adding to the most-recently added row
699 20 : Int row=ct_->nrow()-1;
700 :
701 : // Meta-info
702 20 : ncmc.time().put(row, (begin + end) / 2);
703 20 : ncmc.fieldId().put(row, currField());
704 20 : ncmc.spwId().put(row, currSpw());
705 20 : ncmc.scanNo().put(row, currScan());
706 20 : ncmc.obsId().put(row, currObs());
707 20 : ncmc.interval().put(row, (end - begin) / 2);
708 20 : ncmc.antenna1().put(row, iant);
709 20 : ncmc.antenna2().put(row, -1);
710 :
711 : // Params
712 20 : ncmc.fparam().put(row,gain);
713 :
714 : // Stats
715 20 : ncmc.paramerr().put(row,solveAllParErr().xyPlane(iant));
716 20 : ncmc.snr().put(row,solveAllParErr().xyPlane(iant));
717 20 : ncmc.flag().put(row,!solveAllParOK().xyPlane(iant));
718 : }
719 :
720 : // This spw now has some solutions in it
721 2 : spwOK_(currSpw())=True;
722 : }
723 1 : }
724 :
725 628 : void EPowerCurve::calcAllJones() {
726 :
727 628 : if (prtlev()>6) cout << " EPowerCurve::calcAllJones()" << endl;
728 :
729 : // Nominally no gain curve effect
730 628 : currJElem()=Complex(1.0);
731 628 : currJElemOK()=false;
732 :
733 : /*
734 : cout << "currSpw() = " << currSpw() << endl;
735 : cout << " spwMap() = " << spwMap() << endl;
736 : cout << " currRPar().shape() = " << currRPar().shape() << endl;
737 : cout << " currRPar() = " << currRPar() << endl;
738 : */
739 :
740 628 : Complex* J=currJElem().data();
741 628 : Bool* JOk=currJElemOK().data();
742 628 : Float* c=currRPar().data();
743 628 : Double* a=za().data();
744 :
745 : Double loss, ang;
746 6908 : for (Int iant=0; iant<nAnt(); ++iant,++a)
747 18840 : for (Int ipol=0;ipol<2;++ipol,++J,++JOk) {
748 12560 : loss=Double(*c);
749 12560 : ++c;
750 12560 : ang=1.0;
751 100480 : for (Int ipar=1;ipar<nPar()/2;++ipar,++c) {
752 87920 : ang*=(*a);
753 87920 : loss+=((*c)*ang);
754 : }
755 12560 : if (loss >= 0) {
756 12560 : (*J) = Complex(sqrt(loss));
757 12560 : (*JOk) = true;
758 : } else {
759 0 : (*J) = 0.0;
760 0 : (*JOk) = false;
761 : }
762 : }
763 628 : }
764 :
765 : } //# NAMESPACE CASA - END
|
