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