Line data Source code
1 : //# VPManager.cc: Implementation of VPManager.h
2 : //# Copyright (C) 1996-2016
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //# Copyright by ESO (in the framework of the ALMA collaboration)
5 : //#
6 : //# This program is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU General Public License as published by the Free
8 : //# Software Foundation; either version 2 of the License, or (at your option)
9 : //# any later version.
10 : //#
11 : //# This program is distributed in the hope that it will be useful, but WITHOUT
12 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
14 : //# more details.
15 : //#
16 : //# You should have received a copy of the GNU General Public License along
17 : //# with this program; if not, write to the Free Software Foundation, Inc.,
18 : //# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning AIPS++ should be addressed as follows:
21 : //# Internet email: aips2-request@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : //#
27 : //# $Id$
28 :
29 :
30 : #include <casacore/casa/Containers/Record.h>
31 : #include <casacore/casa/Exceptions/Error.h>
32 : #include <iostream>
33 : #include <casacore/tables/Tables/ScalarColumn.h>
34 : #include <casacore/tables/Tables/ScaColDesc.h>
35 : #include <casacore/tables/Tables/ScaRecordColDesc.h>
36 : #include <casacore/tables/Tables/Table.h>
37 : #include <casacore/tables/Tables/SetupNewTab.h>
38 : #include <casacore/tables/TaQL/TableParse.h>
39 : #include <casacore/tables/Tables/TableDesc.h>
40 : #include <casacore/casa/Quanta/QuantumHolder.h>
41 : #include <casacore/casa/Quanta/Quantum.h>
42 : #include <casacore/measures/Measures/MDirection.h>
43 : #include <casacore/measures/Measures/MeasureHolder.h>
44 : #include <synthesis/MeasurementEquations/VPManager.h>
45 : #include <synthesis/TransformMachines/PBMathInterface.h>
46 : #include <synthesis/TransformMachines/PBMath.h>
47 : #include <synthesis/TransformMachines/SynthesisError.h>
48 : #include <synthesis/TransformMachines/ALMACalcIlluminationConvFunc.h>
49 : #include <casacore/casa/Logging.h>
50 : #include <casacore/casa/Logging/LogIO.h>
51 : #include <casacore/casa/Logging/LogSink.h>
52 : #include <casacore/casa/Logging/LogMessage.h>
53 : #include <casacore/casa/OS/Directory.h>
54 : #include <casacore/images/Images/PagedImage.h>
55 : #include <unistd.h>
56 :
57 : using namespace casacore;
58 : namespace casa { //# NAMESPACE CASA - BEGIN
59 :
60 : VPManager* VPManager::instance_p = 0;
61 : std::recursive_mutex VPManager::mutex_p; // to permit calls in same thread
62 :
63 1 : VPManager::VPManager(Bool verbose):
64 : vplist_p(),
65 1 : aR_p()
66 : {
67 1 : reset(verbose);
68 1 : }
69 :
70 :
71 257 : VPManager* VPManager::Instance(){
72 257 : if(instance_p==0){
73 2 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
74 1 : if(instance_p==0){
75 1 : instance_p = new VPManager();
76 : }
77 : }
78 257 : return instance_p;
79 : }
80 :
81 214 : void VPManager::reset(Bool verbose)
82 : {
83 :
84 428 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
85 :
86 214 : vplist_p = Record();
87 214 : vplistdefaults_p.clear();
88 214 : aR_p.init();
89 :
90 428 : LogIO os;
91 214 : os << LogOrigin("VPManager", "reset");
92 :
93 428 : String telName;
94 8988 : for(Int pbtype = static_cast<Int>(PBMath::DEFAULT) + 1;
95 8988 : pbtype < static_cast<Int>(PBMath::NONE); ++pbtype){
96 8774 : PBMath::nameCommonPB(static_cast<PBMath::CommonPB>(pbtype), telName);
97 8774 : if ( vplistdefaults_p.find(telName) != vplistdefaults_p.end( ) ) vplistdefaults_p[telName] = -1;
98 8774 : else vplistdefaults_p.insert(std::pair<casacore::String, casacore::Int >(telName,-1));
99 : }
100 :
101 : // check for available AntennaResponses tables in the Observatories table
102 428 : Vector<String> obsName = MeasTable::Observatories();
103 11556 : for(uInt i=0; i<obsName.size(); i++){
104 :
105 11342 : String telName = obsName(i);
106 :
107 11342 : String antRespPath;
108 11342 : if(!MeasTable::AntennaResponsesPath(antRespPath, telName)) {
109 : // unknown observatory
110 11342 : continue;
111 : }
112 : else{ // remember the corresponding telescopes as special vplist entries
113 0 : if(!aR_p.init(antRespPath)){
114 0 : if(verbose){
115 : os << LogIO::WARN
116 : << "Invalid path defined in Observatories table for \"" << telName << "\":" << endl
117 : << antRespPath << endl
118 0 : << LogIO::POST;
119 : }
120 : }
121 : else{
122 : // init successful
123 0 : Record rec;
124 0 : rec.define("name", "REFERENCE");
125 0 : rec.define("isVP", PBMathInterface::NONE);
126 0 : rec.define("telescope", telName);
127 0 : rec.define("antresppath", antRespPath);
128 :
129 0 : if(verbose){
130 : os << LogIO::NORMAL
131 : << "Will use " << telName << " antenna responses from table "
132 0 : << antRespPath << LogIO::POST;
133 : }
134 :
135 0 : auto vpiter = vplistdefaults_p.find(telName);
136 0 : if( vpiter != vplistdefaults_p.end( ) ){ // there is already a vp for this telName
137 0 : Int ifield = vpiter->second;
138 0 : if(ifield>=0){
139 0 : Record rrec = vplist_p.rwSubRecord(ifield);
140 0 : rrec.define("dopb", false);
141 : }
142 0 : vplistdefaults_p.erase(vpiter);
143 : }
144 0 : if ( vplistdefaults_p.find(telName) != vplistdefaults_p.end( ) ) vplistdefaults_p[telName] = vplist_p.nfields();
145 0 : else vplistdefaults_p.insert(std::pair<casacore::String, casacore::Int >(telName,vplist_p.nfields()));
146 0 : rec.define("dopb", true);
147 :
148 0 : vplist_p.defineRecord(vplist_p.nfields(), rec);
149 : }
150 : }
151 : }
152 214 : if(verbose){
153 0 : os << LogIO::NORMAL << "VPManager initialized." << LogIO::POST;
154 : }
155 :
156 214 : }
157 :
158 :
159 7 : Bool VPManager::saveastable(const String& tablename){
160 :
161 14 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
162 :
163 21 : TableDesc td("vptable", "1", TableDesc::Scratch);
164 7 : td.addColumn(ScalarColumnDesc<String>("telescope"));
165 7 : td.addColumn(ScalarColumnDesc<Int>("antenna"));
166 7 : td.addColumn(ScalarRecordColumnDesc("pbdescription"));
167 14 : SetupNewTable aTab(tablename, td, Table::New);
168 :
169 14 : Table tb(aTab, Table::Plain, vplist_p.nfields());
170 14 : ScalarColumn<String> telcol(tb, "telescope");
171 14 : ScalarColumn<Int> antcol(tb, "antenna");
172 14 : ScalarColumn<TableRecord> pbcol(tb, "pbdescription");
173 25 : for (uInt k=0; k < vplist_p.nfields(); ++k){
174 36 : TableRecord antRec(vplist_p.asRecord(k));
175 36 : String tel=antRec.asString("telescope");
176 18 : telcol.put(k, tel);
177 18 : antcol.put(k,k);
178 18 : pbcol.put(k,antRec);
179 : }
180 :
181 : // create subtable for the vplistdefaults
182 21 : TableDesc td2("vplistdefaultstable", "1", TableDesc::Scratch);
183 7 : td2.addColumn(ScalarColumnDesc<String>("tel_and_anttype"));
184 7 : td2.addColumn(ScalarColumnDesc<Int>("vplistnum"));
185 14 : SetupNewTable defaultsSetup(tablename+"/VPLIST_DEFAULTS", td2, Table::New);
186 :
187 14 : Table tb2(defaultsSetup, Table::Plain, vplistdefaults_p.size( ));
188 14 : ScalarColumn<String> telcol2(tb2, "tel_and_anttype");
189 7 : ScalarColumn<Int> listnumcol(tb2, "vplistnum");
190 7 : uInt keycount = 0;
191 294 : for ( auto iter = vplistdefaults_p.begin( ); iter != vplistdefaults_p.end( ); ++iter, ++keycount ){
192 287 : telcol2.put(keycount, iter->first);
193 287 : listnumcol.put(keycount, iter->second);
194 : }
195 :
196 7 : tb.rwKeywordSet().defineTable("VPLIST_DEFAULTS", tb2);
197 :
198 7 : tb2.flush();
199 7 : tb.flush();
200 14 : return true;
201 : }
202 :
203 :
204 14 : Bool VPManager::loadfromtable(const String& tablename){
205 :
206 28 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
207 :
208 42 : LogIO os(LogOrigin("vpmanager", "loadfromtable"));
209 :
210 28 : Table tb(tablename);
211 28 : ScalarColumn<TableRecord> pbcol(tb, "pbdescription");
212 :
213 28 : Table tb2;
214 :
215 14 : if (tb.keywordSet().isDefined("VPLIST_DEFAULTS")) {
216 14 : tb2 = tb.keywordSet().asTable("VPLIST_DEFAULTS");
217 : }
218 : else{
219 : os << "Format error: table " << tablename
220 0 : << " does not contain a VPLIST_DEFAULTS subtable." << LogIO::POST;
221 0 : return false;
222 : }
223 :
224 28 : ScalarColumn<String> telcol2(tb2, "tel_and_anttype");
225 28 : ScalarColumn<Int> listnumcol(tb2, "vplistnum");
226 :
227 28 : Record tempvplist;
228 28 : std::map<String, Int> tempvplistdefaults;
229 :
230 47 : for (uInt k=0; k < tb.nrow(); k++){
231 33 : tempvplist.defineRecord(k, Record(pbcol(k)));
232 : }
233 588 : for (uInt k=0; k < tb2.nrow(); k++){
234 574 : Int vplistnum = listnumcol(k);
235 574 : if((vplistnum < -1) || (vplistnum>=(Int)tempvplist.nfields())){
236 : os << "Error: invalid vplist number " << vplistnum
237 : << " in subtable VPLIST_DEFAULTS of table "
238 : << tablename << endl
239 0 : << "Valid value range is -1 to " << (Int)tempvplist.nfields()-1
240 0 : << LogIO::POST;
241 0 : return false;
242 : }
243 574 : if ( tempvplistdefaults.find(telcol2(k)) != tempvplistdefaults.end( ) ) tempvplistdefaults[telcol2(k)] = vplistnum;
244 574 : else tempvplistdefaults.insert(std::pair<casacore::String, casacore::Int >(telcol2(k), vplistnum));
245 : }
246 :
247 : // overwrite existing information
248 14 : vplist_p = tempvplist;
249 14 : vplistdefaults_p = tempvplistdefaults;
250 :
251 : os << "Loaded " << tb.nrow() << " VP definitions and " << tb2.nrow()
252 14 : << " VP default settings from table " << tablename << LogIO:: POST;
253 :
254 14 : return true;
255 :
256 : }
257 :
258 :
259 0 : Bool VPManager::summarizevps(const Bool verbose) {
260 :
261 0 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
262 :
263 0 : LogIO os(LogOrigin("vpmanager", "summarizevps"));
264 :
265 : os << LogIO::NORMAL << "Voltage patterns internally defined in CASA (* = global default for this telescope):\n"
266 : << " Telescope: Class"
267 0 : << LogIO::POST;
268 0 : String telName;
269 0 : String pbClassName;
270 0 : for(Int pbtype = static_cast<Int>(PBMath::DEFAULT) + 1;
271 0 : pbtype < static_cast<Int>(PBMath::NONE); ++pbtype){
272 0 : PBMath::nameCommonPB(static_cast<PBMath::CommonPB>(pbtype), telName);
273 0 : auto vpiter = vplistdefaults_p.find(telName);
274 0 : if ( vpiter != vplistdefaults_p.end( ) && vpiter->second == -1 ) {
275 0 : os << LogIO::NORMAL << " * ";
276 : }
277 : else{
278 0 : os << LogIO::NORMAL << " ";
279 : }
280 0 : os << telName;
281 : try{
282 0 : PBMath(static_cast<PBMath::CommonPB>(pbtype)).namePBClass(pbClassName);
283 0 : os << ": " << pbClassName << LogIO::POST;
284 : }
285 0 : catch(AipsError){
286 0 : os << ": not available" << LogIO::POST;
287 : }
288 : }
289 :
290 0 : os << LogIO::NORMAL << "\nExternally defined voltage patterns (* = global default for this telescope):" << LogIO::POST;
291 0 : if (vplist_p.nfields() > 0) {
292 0 : os << "VP# Tel VP Type " << LogIO::POST;
293 0 : for (uInt i=0; i < vplist_p.nfields(); ++i){
294 0 : TableRecord antRec(vplist_p.asRecord(i));
295 0 : String telName = antRec.asString("telescope");
296 0 : auto vpiter = vplistdefaults_p.find(telName);
297 0 : if ( vpiter != vplistdefaults_p.end( ) && ((Int)i == vpiter->second) ) {
298 0 : os << i << " * ";
299 : }
300 : else{
301 0 : os << i << " ";
302 : }
303 0 : os << String(telName+ " ").resize(11);
304 0 : os << antRec.asString("name");
305 :
306 0 : if(antRec.asString("name")=="REFERENCE"){
307 0 : os << ": " << antRec.asString("antresppath");
308 : }
309 : else{
310 : // antenna types
311 0 : uInt counter=0;
312 0 : os << " (used for antenna types ";
313 0 : for ( auto iter = vplistdefaults_p.begin( ); iter != vplistdefaults_p.end( ); ++iter ) {
314 0 : String aDesc = iter->first;
315 0 : if(telName == telFromAntDesc(aDesc)
316 0 : && ((Int)i == iter->second)
317 : ){
318 0 : if(counter>0){
319 0 : os << ", ";
320 : }
321 0 : if(antTypeFromAntDesc(aDesc).empty()){
322 0 : os << "any";
323 : }
324 : else{
325 0 : os << "\"" << antTypeFromAntDesc(aDesc) << "\"";
326 : }
327 0 : counter++;
328 : }
329 : }
330 0 : os << ")";
331 : }
332 :
333 0 : os << LogIO::POST;
334 0 : if (verbose) {
335 0 : ostringstream oss;
336 0 : antRec.print(oss);
337 0 : os << oss.str() << LogIO::POST;
338 : }
339 : }
340 : } else {
341 0 : os << "\tNone" << LogIO::POST;
342 : }
343 :
344 0 : return true;
345 : }
346 :
347 0 : Bool VPManager::setcannedpb(const String& tel,
348 : const String& other,
349 : const Bool dopb,
350 : const String& commonpb,
351 : const Bool dosquint,
352 : const Quantity& paincrement,
353 : const Bool usesymmetricbeam,
354 : Record& rec){
355 :
356 0 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
357 :
358 0 : rec = Record();
359 0 : rec.define("name", "COMMONPB");
360 0 : rec.define("isVP", PBMathInterface::COMMONPB);
361 0 : if (tel=="OTHER") {
362 0 : rec.define("telescope", other);
363 : } else {
364 0 : rec.define("telescope", tel);
365 : }
366 0 : rec.define("dopb", dopb);
367 0 : rec.define("commonpb", commonpb);
368 0 : rec.define("dosquint", dosquint);
369 0 : String error;
370 0 : Record tempholder;
371 0 : QuantumHolder(paincrement).toRecord(error, tempholder);
372 0 : rec.defineRecord("paincrement", tempholder);
373 0 : rec.define("usesymmetricbeam", usesymmetricbeam);
374 :
375 0 : if(dopb){
376 0 : auto key = rec.asString(rec.fieldNumber("telescope"));
377 0 : if ( vplistdefaults_p.find(key) != vplistdefaults_p.end( ) ) vplistdefaults_p[key] = vplist_p.nfields( );
378 0 : else vplistdefaults_p.insert(std::pair<casacore::String, casacore::Int >(key, vplist_p.nfields( )));
379 : }
380 :
381 0 : vplist_p.defineRecord(vplist_p.nfields(), rec);
382 :
383 0 : return true;
384 : }
385 :
386 0 : Bool VPManager::setpbairy(const String& tel,
387 : const String& other,
388 : const Bool dopb, const Quantity& dishdiam,
389 : const Quantity& blockagediam,
390 : const Quantity& maxrad,
391 : const Quantity& reffreq,
392 : MDirection& squintdir,
393 : const Quantity& squintreffreq,const Bool dosquint,
394 : const Quantity& paincrement,
395 : const Bool usesymmetricbeam,
396 : Record& rec){
397 :
398 0 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
399 :
400 0 : rec=Record();
401 0 : rec.define("name", "AIRY");
402 0 : rec.define("isVP", PBMathInterface::AIRY);
403 0 : if(tel=="OTHER"){
404 0 : rec.define("telescope", other);
405 : }
406 : else{
407 0 : rec.define("telescope", tel);
408 : }
409 0 : rec.define("dopb", dopb);
410 0 : String error;
411 0 : {Record tempholder;
412 0 : QuantumHolder(dishdiam).toRecord(error, tempholder);
413 0 : rec.defineRecord("dishdiam", tempholder);
414 : }
415 0 : {Record tempholder;
416 0 : QuantumHolder(blockagediam).toRecord(error, tempholder);
417 0 : rec.defineRecord("blockagediam", tempholder);
418 : }
419 0 : { Record tempholder;
420 0 : QuantumHolder(maxrad).toRecord(error, tempholder);
421 0 : rec.defineRecord("maxrad", tempholder);
422 : }
423 0 : {Record tempholder;
424 0 : QuantumHolder(reffreq).toRecord(error, tempholder);
425 0 : rec.defineRecord("reffreq", tempholder);
426 : }
427 0 : rec.define("isthisvp", false);
428 0 : {Record tempholder;
429 0 : MeasureHolder(squintdir).toRecord(error, tempholder);
430 0 : rec.defineRecord("squintdir", tempholder);
431 : }
432 0 : {Record tempholder;
433 0 : QuantumHolder(squintreffreq).toRecord(error, tempholder);
434 0 : rec.defineRecord("squintreffreq", tempholder);
435 : }
436 0 : rec.define("dosquint", dosquint);
437 0 : {Record tempholder;
438 0 : QuantumHolder(paincrement).toRecord(error, tempholder);
439 0 : rec.defineRecord("paincrement", tempholder);
440 : }
441 0 : rec.define("usesymmetricbeam", usesymmetricbeam);
442 :
443 0 : if(dopb){
444 0 : auto key = rec.asString(rec.fieldNumber("telescope"));
445 0 : if ( vplistdefaults_p.find(key) != vplistdefaults_p.end( ) ) vplistdefaults_p[key] = vplist_p.nfields( );
446 0 : else vplistdefaults_p.insert(std::pair<casacore::String, casacore::Int >(key, vplist_p.nfields( )));
447 : }
448 :
449 0 : vplist_p.defineRecord(vplist_p.nfields(), rec);
450 :
451 0 : return true;
452 :
453 : }
454 :
455 0 : Bool VPManager::setpbcospoly(const String& telescope,
456 : const String& othertelescope,
457 : const Bool dopb, const Vector<Double>& coeff,
458 : const Vector<Double>& scale,
459 : const Quantity& maxrad,
460 : const Quantity& reffreq,
461 : const String& isthispb,
462 : MDirection& squintdir,
463 : const Quantity& squintreffreq,
464 : const Bool dosquint,
465 : const Quantity& paincrement,
466 : const Bool usesymmetricbeam,
467 : Record& rec) {
468 :
469 0 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
470 :
471 0 : rec=Record();
472 0 : rec.define("name", "COSPOLY");
473 0 : rec.define("isVP", PBMathInterface::COSPOLY);
474 0 : if(telescope=="OTHER"){
475 0 : rec.define("telescope", othertelescope);
476 : }
477 : else{
478 0 : rec.define("telescope", telescope);
479 : }
480 0 : rec.define("dopb", dopb);
481 0 : rec.define("coeff", coeff);
482 0 : rec.define("scale", scale);
483 0 : String error;
484 0 : { Record tempholder;
485 0 : QuantumHolder(maxrad).toRecord(error, tempholder);
486 0 : rec.defineRecord("maxrad", tempholder);
487 : }
488 0 : {Record tempholder;
489 0 : QuantumHolder(reffreq).toRecord(error, tempholder);
490 0 : rec.defineRecord("reffreq", tempholder);
491 : }
492 0 : if(isthispb=="PB" || isthispb=="pb"){
493 0 : rec.define("isthisvp", false);
494 : }
495 : else{
496 0 : rec.define("isthisvp", true);
497 : }
498 0 : {Record tempholder;
499 0 : MeasureHolder(squintdir).toRecord(error, tempholder);
500 0 : rec.defineRecord("squintdir", tempholder);
501 : }
502 0 : {Record tempholder;
503 0 : QuantumHolder(squintreffreq).toRecord(error, tempholder);
504 0 : rec.defineRecord("squintreffreq", tempholder);
505 : }
506 0 : rec.define("dosquint", dosquint);
507 0 : {Record tempholder;
508 0 : QuantumHolder(paincrement).toRecord(error, tempholder);
509 0 : rec.defineRecord("paincrement", tempholder);
510 : }
511 0 : rec.define("usesymmetricbeam", usesymmetricbeam);
512 :
513 0 : if(dopb){
514 0 : auto key = rec.asString(rec.fieldNumber("telescope"));
515 0 : if ( vplistdefaults_p.find(key) != vplistdefaults_p.end( ) ) vplistdefaults_p[key] = vplist_p.nfields( );
516 0 : else vplistdefaults_p.insert(std::pair<casacore::String, casacore::Int >(key, vplist_p.nfields( )));
517 : }
518 :
519 0 : vplist_p.defineRecord(vplist_p.nfields(), rec);
520 :
521 0 : return true;
522 :
523 : }
524 :
525 0 : Bool VPManager::setpbgauss(const String& tel, const String& other,
526 : const Bool dopb,
527 : const Quantity& halfwidth, const Quantity maxrad,
528 : const Quantity& reffreq,
529 : const String& isthispb,
530 : MDirection& squintdir,
531 : const Quantity& squintreffreq,
532 : const Bool dosquint,
533 : const Quantity& paincrement,
534 : const Bool usesymmetricbeam,
535 : Record& rec){
536 :
537 0 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
538 :
539 0 : rec=Record();
540 0 : rec.define("name", "GAUSS");
541 0 : rec.define("isVP", PBMathInterface::GAUSS);
542 0 : if(tel=="OTHER"){
543 0 : rec.define("telescope", other);
544 : }
545 : else{
546 0 : rec.define("telescope", tel);
547 : }
548 0 : rec.define("dopb", dopb);
549 0 : Quantity hpw=halfwidth;
550 0 : if(isthispb=="PB" || isthispb=="pb"){
551 0 : hpw.setValue(hpw.getValue()/2.0);
552 : }
553 0 : String error;
554 0 : {Record tempholder;
555 0 : QuantumHolder(hpw).toRecord(error, tempholder);
556 0 : rec.defineRecord("halfwidth", tempholder);
557 : }
558 0 : {Record tempholder;
559 0 : QuantumHolder(maxrad).toRecord(error, tempholder);
560 0 : rec.defineRecord("maxrad", tempholder);
561 : }
562 0 : {Record tempholder;
563 0 : QuantumHolder(reffreq).toRecord(error, tempholder);
564 0 : rec.defineRecord("reffreq", tempholder);
565 : }
566 0 : if(isthispb=="PB" || isthispb=="pb"){
567 0 : rec.define("isthisvp", false);
568 : }
569 : else{
570 0 : rec.define("isthisvp", true);
571 : }
572 0 : {Record tempholder;
573 0 : MeasureHolder(squintdir).toRecord(error, tempholder);
574 0 : rec.defineRecord("squintdir", tempholder);
575 : }
576 0 : {Record tempholder;
577 0 : QuantumHolder(squintreffreq).toRecord(error, tempholder);
578 0 : rec.defineRecord("squintreffreq", tempholder);
579 : }
580 0 : rec.define("dosquint", dosquint);
581 0 : {Record tempholder;
582 0 : QuantumHolder(paincrement).toRecord(error, tempholder);
583 0 : rec.defineRecord("paincrement", tempholder);
584 : }
585 0 : rec.define("usesymmetricbeam", usesymmetricbeam);
586 :
587 0 : if(dopb){
588 0 : auto key = rec.asString(rec.fieldNumber("telescope"));
589 0 : if ( vplistdefaults_p.find(key) != vplistdefaults_p.end( ) ) vplistdefaults_p[key] = vplist_p.nfields( );
590 0 : else vplistdefaults_p.insert(std::pair<casacore::String, casacore::Int >(key, vplist_p.nfields( )));
591 : }
592 :
593 0 : vplist_p.defineRecord(vplist_p.nfields(), rec);
594 :
595 0 : return true;
596 :
597 : }
598 :
599 0 : Bool VPManager::setpbinvpoly(const String& telescope,
600 : const String& othertelescope,
601 : const Bool dopb, const Vector<Double>& coeff,
602 : const Quantity& maxrad,
603 : const Quantity& reffreq,
604 : const String& isthispb,
605 : MDirection& squintdir,
606 : const Quantity& squintreffreq,
607 : const Bool dosquint,
608 : const Quantity& paincrement,
609 : const Bool usesymmetricbeam,
610 : Record& rec) {
611 :
612 0 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
613 :
614 0 : rec=Record();
615 0 : rec.define("name", "IPOLY");
616 0 : rec.define("isVP", PBMathInterface::IPOLY);
617 0 : if(telescope=="OTHER"){
618 0 : rec.define("telescope", othertelescope);
619 : }
620 : else{
621 0 : rec.define("telescope", telescope);
622 : }
623 0 : rec.define("dopb", dopb);
624 0 : rec.define("coeff", coeff);
625 0 : String error;
626 0 : {Record tempholder;
627 0 : QuantumHolder(maxrad).toRecord(error, tempholder);
628 0 : rec.defineRecord("maxrad", tempholder);
629 : }
630 0 : {Record tempholder;
631 0 : QuantumHolder(reffreq).toRecord(error, tempholder);
632 0 : rec.defineRecord("reffreq", tempholder);
633 : }
634 0 : if(isthispb=="PB" || isthispb=="pb"){
635 0 : rec.define("isthisvp", false);
636 : }
637 : else{
638 0 : rec.define("isthisvp", true);
639 : }
640 0 : {Record tempholder;
641 0 : MeasureHolder(squintdir).toRecord(error, tempholder);
642 0 : rec.defineRecord("squintdir", tempholder);
643 : }
644 0 : {Record tempholder;
645 0 : QuantumHolder(squintreffreq).toRecord(error, tempholder);
646 0 : rec.defineRecord("squintreffreq", tempholder);
647 : }
648 0 : rec.define("dosquint", dosquint);
649 0 : {Record tempholder;
650 0 : QuantumHolder(paincrement).toRecord(error, tempholder);
651 0 : rec.defineRecord("paincrement", tempholder);
652 : }
653 0 : rec.define("usesymmetricbeam", usesymmetricbeam);
654 :
655 0 : if(dopb){
656 0 : auto key = rec.asString(rec.fieldNumber("telescope"));
657 0 : if ( vplistdefaults_p.find(key) != vplistdefaults_p.end( ) ) vplistdefaults_p[key] = vplist_p.nfields( );
658 0 : else vplistdefaults_p.insert(std::pair<casacore::String, casacore::Int >(key, vplist_p.nfields( )));
659 : }
660 :
661 0 : vplist_p.defineRecord(vplist_p.nfields(), rec);
662 :
663 0 : return true;
664 :
665 :
666 : }
667 0 : Bool VPManager::setpbinvpoly(const String& telescope,
668 : const String& othertelescope,
669 : const Bool dopb, const Matrix<Double>& coeff,
670 : const Vector<Double>& freqs,
671 : const Quantity& maxrad,
672 : const Quantity& reffreq,
673 : const String& isthispb,
674 : MDirection& squintdir,
675 : const Quantity& squintreffreq,
676 : const Bool dosquint,
677 : const Quantity& paincrement,
678 : const Bool usesymmetricbeam,
679 : Record& rec) {
680 :
681 0 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
682 :
683 0 : rec=Record();
684 0 : rec.define("name", "IPOLY");
685 0 : rec.define("isVP", PBMathInterface::IPOLY);
686 0 : if(telescope=="OTHER"){
687 0 : rec.define("telescope", othertelescope);
688 : }
689 : else{
690 0 : rec.define("telescope", telescope);
691 : }
692 0 : rec.define("dopb", dopb);
693 0 : rec.define("coeff", coeff);
694 0 : rec.define("fitfreqs", freqs);
695 0 : String error;
696 0 : {Record tempholder;
697 0 : QuantumHolder(maxrad).toRecord(error, tempholder);
698 0 : rec.defineRecord("maxrad", tempholder);
699 : }
700 0 : {Record tempholder;
701 0 : QuantumHolder(reffreq).toRecord(error, tempholder);
702 0 : rec.defineRecord("reffreq", tempholder);
703 : }
704 0 : if(isthispb=="PB" || isthispb=="pb"){
705 0 : rec.define("isthisvp", false);
706 : }
707 : else{
708 0 : rec.define("isthisvp", true);
709 : }
710 0 : {Record tempholder;
711 0 : MeasureHolder(squintdir).toRecord(error, tempholder);
712 0 : rec.defineRecord("squintdir", tempholder);
713 : }
714 0 : {Record tempholder;
715 0 : QuantumHolder(squintreffreq).toRecord(error, tempholder);
716 0 : rec.defineRecord("squintreffreq", tempholder);
717 : }
718 0 : rec.define("dosquint", dosquint);
719 0 : {Record tempholder;
720 0 : QuantumHolder(paincrement).toRecord(error, tempholder);
721 0 : rec.defineRecord("paincrement", tempholder);
722 : }
723 0 : rec.define("usesymmetricbeam", usesymmetricbeam);
724 :
725 0 : if(dopb){
726 0 : auto key = rec.asString(rec.fieldNumber("telescope"));
727 0 : if ( vplistdefaults_p.find(key) != vplistdefaults_p.end( ) ) vplistdefaults_p[key] = vplist_p.nfields( );
728 0 : else vplistdefaults_p.insert(std::pair<casacore::String, casacore::Int >(key, vplist_p.nfields( )));
729 : }
730 :
731 0 : vplist_p.defineRecord(vplist_p.nfields(), rec);
732 :
733 0 : return true;
734 :
735 :
736 : }
737 :
738 0 : Bool VPManager::setpbnumeric(const String& telescope,
739 : const String& othertelescope,
740 : const Bool dopb, const Vector<Double>& vect,
741 : const Quantity& maxrad,
742 : const Quantity& reffreq,
743 : const String& isthispb,
744 : MDirection& squintdir,
745 : const Quantity& squintreffreq,
746 : const Bool dosquint,
747 : const Quantity& paincrement,
748 : const Bool usesymmetricbeam,
749 : Record &rec) {
750 :
751 0 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
752 :
753 0 : rec=Record();
754 0 : rec.define("name", "NUMERIC");
755 0 : rec.define("isVP", PBMathInterface::NUMERIC);
756 0 : if(telescope=="OTHER"){
757 0 : rec.define("telescope", othertelescope);
758 : }
759 : else{
760 0 : rec.define("telescope", telescope);
761 : }
762 0 : rec.define("dopb", dopb);
763 0 : rec.define("vect", vect);
764 0 : String error;
765 0 : {Record tempholder;
766 0 : QuantumHolder(maxrad).toRecord(error, tempholder);
767 0 : rec.defineRecord("maxrad", tempholder);
768 : }
769 0 : {Record tempholder;
770 0 : QuantumHolder(reffreq).toRecord(error, tempholder);
771 0 : rec.defineRecord("reffreq", tempholder);
772 : }
773 0 : if(isthispb=="PB" || isthispb=="pb"){
774 0 : rec.define("isthisvp", false);
775 : }
776 : else{
777 0 : rec.define("isthisvp", true);
778 : }
779 0 : {Record tempholder;
780 0 : MeasureHolder(squintdir).toRecord(error, tempholder);
781 0 : rec.defineRecord("squintdir", tempholder);
782 : }
783 0 : {Record tempholder;
784 0 : QuantumHolder(squintreffreq).toRecord(error, tempholder);
785 0 : rec.defineRecord("squintreffreq", tempholder);
786 : }
787 0 : rec.define("dosquint", dosquint);
788 0 : {Record tempholder;
789 0 : QuantumHolder(paincrement).toRecord(error, tempholder);
790 0 : rec.defineRecord("paincrement", tempholder);
791 : }
792 0 : rec.define("usesymmetricbeam", usesymmetricbeam);
793 :
794 0 : if(dopb){
795 0 : auto key = rec.asString(rec.fieldNumber("telescope"));
796 0 : if ( vplistdefaults_p.find(key) != vplistdefaults_p.end( ) ) vplistdefaults_p[key] = vplist_p.nfields( );
797 0 : else vplistdefaults_p.insert(std::pair<casacore::String, casacore::Int >(key, vplist_p.nfields( )));
798 : }
799 :
800 0 : vplist_p.defineRecord(vplist_p.nfields(), rec);
801 :
802 0 : return true;
803 :
804 : }
805 :
806 0 : Bool VPManager::setpbimage(const String& tel,
807 : const String& other,
808 : const Bool dopb,
809 : const String& realimage,
810 : const String& imagimage,
811 : const String& compleximage,
812 : const Vector<String>& antennaNames,
813 : Record& rec){
814 :
815 0 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
816 :
817 0 : rec=Record();
818 0 : rec.define("name", "IMAGE");
819 0 : rec.define("isVP", PBMathInterface::IMAGE);
820 0 : if(tel=="OTHER"){
821 0 : rec.define("telescope", other);
822 : }
823 : else{
824 0 : rec.define("telescope", tel);
825 : }
826 0 : rec.define("dopb", dopb);
827 0 : rec.define("isthisvp", true);
828 0 : if(compleximage==""){
829 0 : rec.define("realimage", realimage);
830 0 : rec.define("imagimage", imagimage);
831 : }
832 : else{
833 0 : rec.define("compleximage", compleximage);
834 : }
835 0 : rec.define("antennanames", antennaNames);
836 :
837 0 : if(dopb){
838 0 : auto key = rec.asString(rec.fieldNumber("telescope"));
839 0 : if ( vplistdefaults_p.find(key) != vplistdefaults_p.end( ) ) vplistdefaults_p[key] = vplist_p.nfields( );
840 0 : else vplistdefaults_p.insert(std::pair<casacore::String, casacore::Int >(key, vplist_p.nfields( )));
841 : }
842 :
843 0 : vplist_p.defineRecord(vplist_p.nfields(), rec);
844 :
845 0 : return true;
846 :
847 : }
848 :
849 0 : Bool VPManager::imagepbinfo(Vector<Vector<String> >& names, Vector<Record>& imagebeams){
850 :
851 0 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
852 0 : Int nRec=vplist_p.nfields();
853 0 : Bool retval=false;
854 0 : names=Vector<Vector<String> >(0);
855 0 : imagebeams=Vector<Record>(0);
856 0 : for(Int k=0; k< nRec; ++k){
857 0 : Record elrec=vplist_p.asRecord(k);
858 :
859 0 : if(elrec.isDefined("name") && elrec.asString("name")== String("IMAGE")){
860 0 : names.resize(names.nelements()+1, true);
861 0 : imagebeams.resize(names.nelements(), true);
862 0 : Vector<String> localstr;
863 0 : elrec.get("antennanames", localstr);
864 0 : names(names.nelements()-1)=localstr;
865 0 : imagebeams[names.nelements()-1]=elrec;
866 0 : retval=true;
867 : }
868 : }
869 :
870 0 : return retval;
871 :
872 :
873 : }
874 :
875 7 : Bool VPManager::setpbpoly(const String& telescope,
876 : const String& othertelescope,
877 : const Bool dopb, const Vector<Double>& coeff,
878 : const Quantity& maxrad,
879 : const Quantity& reffreq,
880 : const String& isthispb,
881 : MDirection& squintdir,
882 : const Quantity& squintreffreq, const Bool dosquint,
883 : const Quantity& paincrement,
884 : const Bool usesymmetricbeam,
885 : Record &rec) {
886 :
887 14 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
888 :
889 7 : rec=Record();
890 7 : rec.define("name", "POLY");
891 7 : rec.define("isVP", PBMathInterface::POLY);
892 7 : if(telescope=="OTHER"){
893 0 : rec.define("telescope", othertelescope);
894 : }
895 : else{
896 7 : rec.define("telescope", telescope);
897 : }
898 7 : rec.define("dopb", dopb);
899 7 : rec.define("coeff", coeff);
900 7 : String error;
901 : {
902 7 : Record tempholder;
903 7 : QuantumHolder(maxrad).toRecord(error, tempholder);
904 7 : rec.defineRecord("maxrad", tempholder);
905 : }
906 7 : {Record tempholder;
907 7 : QuantumHolder(reffreq).toRecord(error, tempholder);
908 7 : rec.defineRecord("reffreq", tempholder);
909 : }
910 7 : if(isthispb=="PB" || isthispb=="pb"){
911 7 : rec.define("isthisvp", false);
912 : }
913 : else{
914 0 : rec.define("isthisvp", true);
915 : }
916 7 : {Record tempholder;
917 7 : MeasureHolder(squintdir).toRecord(error, tempholder);
918 7 : rec.defineRecord("squintdir", tempholder);
919 : }
920 7 : {Record tempholder;
921 7 : QuantumHolder(squintreffreq).toRecord(error, tempholder);
922 7 : rec.defineRecord("squintreffreq", tempholder);
923 : }
924 7 : rec.define("dosquint", dosquint);
925 7 : {Record tempholder;
926 7 : QuantumHolder(paincrement).toRecord(error, tempholder);
927 7 : rec.defineRecord("paincrement", tempholder);
928 : }
929 7 : rec.define("usesymmetricbeam", usesymmetricbeam);
930 :
931 7 : if(dopb){
932 21 : auto key = rec.asString(rec.fieldNumber("telescope"));
933 7 : if ( vplistdefaults_p.find(key) != vplistdefaults_p.end( ) ) vplistdefaults_p[key] = vplist_p.nfields( );
934 0 : else vplistdefaults_p.insert(std::pair<casacore::String, casacore::Int >(key, vplist_p.nfields( )));
935 : }
936 :
937 7 : vplist_p.defineRecord(vplist_p.nfields(), rec);
938 :
939 14 : return true;
940 :
941 : }
942 :
943 0 : Bool VPManager::setpbantresptable(const String& telescope, const String& othertelescope,
944 : const Bool dopb, const String& tablepath){
945 :
946 0 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
947 :
948 0 : Record rec;
949 0 : rec.define("name", "REFERENCE");
950 0 : rec.define("isVP", PBMathInterface::NONE);
951 0 : if(telescope=="OTHER"){
952 0 : rec.define("telescope", othertelescope);
953 : }
954 : else{
955 0 : rec.define("telescope", telescope);
956 : }
957 0 : rec.define("dopb", dopb);
958 0 : rec.define("antresppath", tablepath);
959 :
960 0 : if(dopb){
961 0 : auto key = rec.asString(rec.fieldNumber("telescope"));
962 0 : if ( vplistdefaults_p.find(key) != vplistdefaults_p.end( ) ) vplistdefaults_p[key] = vplist_p.nfields( );
963 0 : else vplistdefaults_p.insert(std::pair<casacore::String, casacore::Int >(key, vplist_p.nfields( )));
964 : }
965 :
966 0 : vplist_p.defineRecord(vplist_p.nfields(), rec);
967 :
968 0 : return true;
969 :
970 : }
971 :
972 :
973 0 : Bool VPManager::setuserdefault(const Int vplistfield, // (-1 = reset to standard default, -2 = unset)
974 : const String& telescope,
975 : const String& antennatype){
976 :
977 0 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
978 :
979 0 : LogIO os;
980 0 : os << LogOrigin("VPManager", "setuserdefault");
981 :
982 0 : if((vplistfield < -2) || ((Int)(vplist_p.nfields()) <= vplistfield)){
983 : os << LogIO::SEVERE << "Vplist number " << vplistfield << " invalid." << endl
984 0 : << "Valid entries are -2 (none), -1 (default), up to " << vplist_p.nfields()-1
985 0 : << LogIO::POST;
986 0 : return false;
987 : }
988 :
989 0 : String antennaDesc = antennaDescription(telescope, antennatype);
990 :
991 0 : if(vplistfield>=0){
992 0 : const Record rec = vplist_p.subRecord(vplistfield);
993 : // check if this is a valid VP for this telescope
994 0 : String telName;
995 0 : const Int telFieldNumber=rec.fieldNumber("telescope");
996 0 : if (telFieldNumber!=-1){
997 0 : telName = rec.asString(telFieldNumber);
998 0 : if(telFromAntDesc(telName)!=telescope){
999 : os << LogIO::SEVERE << " entry " << vplistfield << " does not point ot a valid VP for " << telescope
1000 0 : << LogIO::POST;
1001 0 : return false;
1002 : }
1003 : }
1004 0 : Record srec = vplist_p.rwSubRecord(vplistfield);
1005 0 : srec.define("dopb", true);
1006 : }
1007 : // unset set an existing default
1008 0 : auto vpiter = vplistdefaults_p.find(antennaDesc);
1009 0 : if( vpiter != vplistdefaults_p.end( ) ){
1010 0 : vplistdefaults_p.erase(vpiter);
1011 : }
1012 0 : if(vplistfield>-2){ // (-2 means don't set a default)
1013 0 : if ( vplistdefaults_p.find(antennaDesc) != vplistdefaults_p.end( ) ) vplistdefaults_p[antennaDesc] = vplistfield;
1014 0 : else vplistdefaults_p.insert(std::pair<casacore::String, casacore::Int >(antennaDesc,vplistfield));
1015 : }
1016 :
1017 0 : return true;
1018 :
1019 : }
1020 :
1021 441 : Bool VPManager::getuserdefault(Int& vplistfield,
1022 : const String& telescope,
1023 : const String& antennatype){
1024 :
1025 882 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
1026 :
1027 882 : String antDesc = antennaDescription(telescope, antennatype);
1028 :
1029 441 : auto vpiter = vplistdefaults_p.find(antDesc);
1030 441 : if( vpiter != vplistdefaults_p.end( ) ){
1031 428 : vplistfield = vpiter->second;
1032 : } else {
1033 13 : auto tpiter = vplistdefaults_p.find(telescope);
1034 13 : if( tpiter != vplistdefaults_p.end( ) ) { // found global entry
1035 0 : vplistfield = tpiter->second;
1036 : } else {
1037 13 : vplistfield = -2;
1038 13 : return false;
1039 : }
1040 : }
1041 :
1042 428 : return true;
1043 :
1044 : }
1045 :
1046 :
1047 :
1048 : // fill vector with the names of the antenna types with available voltage patterns satisfying the given constraints
1049 0 : Bool VPManager::getanttypes(Vector<String>& anttypes,
1050 : const String& telescope,
1051 : const MEpoch& obstime,
1052 : const MFrequency& freq,
1053 : const MDirection& obsdirection // default: Zenith
1054 : ){
1055 :
1056 0 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
1057 :
1058 0 : LogIO os;
1059 0 : os << LogOrigin("VPManager", "getanttypes");
1060 :
1061 0 : Bool rval=false;
1062 :
1063 0 : anttypes.resize(0);
1064 :
1065 0 : Int ifield = -2;
1066 0 : Bool isReference = true;
1067 :
1068 : // check for global response
1069 0 : getuserdefault(ifield,telescope,"");
1070 :
1071 0 : if(ifield==-1){ // internally defined PB does not distinguish antenna types
1072 0 : anttypes.resize(1);
1073 0 : anttypes(0) = "";
1074 0 : rval = true;
1075 : }
1076 0 : else if(ifield>=0){ // externally defined PB
1077 0 : TableRecord antRec(vplist_p.asRecord(ifield));
1078 0 : String thename = antRec.asString("name");
1079 0 : if(thename=="REFERENCE"){ // points to an AntennaResponses table
1080 : // query the antenna responses
1081 0 : String antRespPath = antRec.asString("antresppath");
1082 0 : if(!aR_p.init(antRespPath)){
1083 : os << LogIO::SEVERE
1084 : << "Invalid path defined in vpmanager for \"" << telescope << "\":" << endl
1085 : << antRespPath << endl
1086 0 : << LogIO::POST;
1087 : }
1088 : else{ // init successful
1089 :
1090 : // construct a proper MFrequency
1091 : MFrequency::Types fromFrameType;
1092 0 : MFrequency::getType(fromFrameType, freq.getRefString());
1093 0 : MPosition obsPos;
1094 0 : MFrequency::Ref fromFrame;
1095 0 : MFrequency mFreq = freq;
1096 0 : if(fromFrameType!=MFrequency::TOPO){
1097 0 : if(!MeasTable::Observatory(obsPos,telescope)){
1098 : os << LogIO::SEVERE << "\"" << telescope << "\" is not listed in the Observatories table."
1099 0 : << LogIO::POST;
1100 0 : return false;
1101 : }
1102 0 : fromFrame = MFrequency::Ref(fromFrameType, MeasFrame(obsdirection, obsPos, obstime));
1103 0 : mFreq = MFrequency(freq.get(Unit("Hz")), fromFrame);
1104 : }
1105 :
1106 0 : if(aR_p.getAntennaTypes(anttypes,
1107 : telescope, // (the observatory name, e.g. "ALMA" or "ACA")
1108 : obstime,
1109 : mFreq,
1110 0 : AntennaResponses::ANY, // the requested function type
1111 : obsdirection)){ // success
1112 0 : rval = true;
1113 : }
1114 : }
1115 : }
1116 : else{ // we don't have a reference response
1117 0 : isReference = false;
1118 : }
1119 : }
1120 :
1121 0 : if(ifield==-2 or !isReference){ // no global response or reference
1122 0 : uInt count = 0;
1123 0 : for( auto iter = vplistdefaults_p.begin( ); iter != vplistdefaults_p.end( ); ++iter ) {
1124 0 : String aDesc = iter->first;
1125 0 : if(telescope == telFromAntDesc(aDesc)){
1126 0 : String aType = antTypeFromAntDesc(aDesc);
1127 0 : Bool tFound = false;
1128 0 : for(uInt j=0; j<anttypes.size(); j++){
1129 0 : if(aType==anttypes(j)){ // already in list?
1130 0 : tFound = true;
1131 0 : break;
1132 : }
1133 : }
1134 0 : if(!tFound){
1135 0 : rval = true;
1136 0 : count++;
1137 0 : anttypes.resize(count, true);
1138 0 : anttypes(count-1) = aType;
1139 : }
1140 : }
1141 : } // end for i
1142 : }
1143 :
1144 0 : return rval;
1145 :
1146 : }
1147 :
1148 :
1149 : // return number of voltage patterns satisfying the given constraints
1150 0 : Int VPManager::numvps(const String& telescope,
1151 : const MEpoch& obstime,
1152 : const MFrequency& freq,
1153 : const MDirection& obsdirection // default: Zenith
1154 : ){
1155 :
1156 0 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
1157 :
1158 0 : LogIO os;
1159 0 : os << LogOrigin("VPManager", "numvps");
1160 :
1161 0 : Vector<String> antTypes;
1162 :
1163 0 : getanttypes(antTypes, telescope, obstime, freq, obsdirection);
1164 :
1165 0 : return antTypes.size();
1166 :
1167 : }
1168 :
1169 :
1170 : // get the voltage pattern satisfying the given constraints
1171 207 : Bool VPManager::getvp(Record &rec,
1172 : const String& telescope,
1173 : const MEpoch& obstime,
1174 : const MFrequency& freq,
1175 : const String& antennatype, // default: ""
1176 : const MDirection& obsdirection){ // default is the Zenith
1177 :
1178 414 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
1179 :
1180 414 : LogIO os;
1181 207 : os << LogOrigin("VPManager", "getvp");
1182 :
1183 207 : Int ifield = -2;
1184 207 : if(!getuserdefault(ifield,telescope,antennatype)){
1185 0 : return false;
1186 : }
1187 :
1188 207 : rec = Record();
1189 207 : Int rval=false;
1190 :
1191 414 : String antDesc = antennaDescription(telescope, antennatype);
1192 :
1193 207 : if(ifield==-1){ // internally defined PB does not distinguish antenna types
1194 :
1195 207 : rec.define("name", "COMMONPB");
1196 207 : rec.define("isVP", PBMathInterface::COMMONPB);
1197 207 : rec.define("telescope", telescope);
1198 207 : rec.define("dopb", true);
1199 207 : rec.define("commonpb", telescope);
1200 207 : rec.define("dosquint", false);
1201 414 : String error;
1202 207 : Record tempholder;
1203 207 : QuantumHolder(Quantity(10.,"deg")).toRecord(error, tempholder);
1204 207 : rec.defineRecord("paincrement", tempholder);
1205 207 : rec.define("usesymmetricbeam", false);
1206 :
1207 207 : rval = true;
1208 :
1209 : }
1210 0 : else if(ifield>=0){ // externally defined PB
1211 0 : TableRecord antRec(vplist_p.asRecord(ifield));
1212 0 : String thename = antRec.asString("name");
1213 0 : if(thename=="REFERENCE"){ // points to an AntennaResponses table
1214 :
1215 : // query the antenna responses
1216 0 : String antRespPath = antRec.asString("antresppath");
1217 0 : if(!aR_p.init(antRespPath)){
1218 : os << LogIO::SEVERE
1219 : << "Invalid path defined in vpmanager for \"" << telescope << "\":" << endl
1220 : << antRespPath << endl
1221 0 : << LogIO::POST;
1222 0 : return false;
1223 : }
1224 : // init successful
1225 0 : String functionImageName;
1226 : uInt funcChannel;
1227 0 : MFrequency nomFreq;
1228 0 : MFrequency loFreq;
1229 0 : MFrequency hiFreq;
1230 : AntennaResponses::FuncTypes fType;
1231 0 : MVAngle rotAngOffset;
1232 :
1233 : // construct a proper MFrequency
1234 : MFrequency::Types fromFrameType;
1235 0 : MFrequency::getType(fromFrameType, freq.getRefString());
1236 0 : MPosition obsPos;
1237 0 : MFrequency::Ref fromFrame;
1238 0 : MFrequency mFreq = freq;
1239 0 : if(fromFrameType!=MFrequency::TOPO){
1240 0 : if(!MeasTable::Observatory(obsPos,telescope)){
1241 : os << LogIO::SEVERE << "\"" << telescope << "\" is not listed in the Observatories table."
1242 0 : << LogIO::POST;
1243 0 : return false;
1244 : }
1245 0 : fromFrame = MFrequency::Ref(fromFrameType, MeasFrame(obsdirection, obsPos, obstime));
1246 0 : mFreq = MFrequency(freq.get(Unit("Hz")), fromFrame);
1247 : }
1248 :
1249 0 : if(!aR_p.getImageName(functionImageName, // the path to the image
1250 : funcChannel, // the channel to use in the image
1251 : nomFreq, // nominal frequency of the image (in the given channel)
1252 : loFreq, // lower end of the validity range
1253 : hiFreq, // upper end of the validity range
1254 : fType, // the function type of the image
1255 : rotAngOffset, // the response rotation angle offset
1256 : /////////////////////
1257 : telescope,
1258 : obstime,
1259 : mFreq,
1260 0 : AntennaResponses::ANY, // the requested function type
1261 : antennatype,
1262 : obsdirection)
1263 : ){
1264 0 : rec = Record();
1265 0 : return false;
1266 : }
1267 :
1268 : // getImageName was successful
1269 :
1270 : // construct record
1271 0 : rec = Record();
1272 0 : Unit uHz("Hz");
1273 0 : switch(fType){
1274 0 : case AntennaResponses::AIF: // complex aperture illumination function
1275 : os << LogIO::WARN << "Responses type AIF provided for " << telescope << " in " << endl
1276 : << antRespPath << endl
1277 : << " not yet supported."
1278 0 : << LogIO::POST;
1279 0 : rval = false;
1280 0 : break;
1281 0 : case AntennaResponses::EFP: // complex electric field pattern
1282 0 : rec.define("name", "IMAGE");
1283 0 : rec.define("isVP", PBMathInterface::IMAGE);
1284 0 : rec.define("telescope", telescope);
1285 0 : rec.define("dopb", true);
1286 0 : rec.define("isthisvp", true);
1287 0 : rec.define("compleximage", functionImageName);
1288 0 : rec.define("channel", funcChannel);
1289 0 : rec.define("reffreq", nomFreq.get(uHz).getValue());
1290 0 : rec.define("minvalidfreq", loFreq.get(uHz).getValue());
1291 0 : rec.define("maxvalidfreq", hiFreq.get(uHz).getValue());
1292 0 : rval = true;
1293 0 : break;
1294 0 : case AntennaResponses::VP: // real voltage pattern
1295 0 : rec.define("name", "IMAGE");
1296 0 : rec.define("isVP", PBMathInterface::IMAGE);
1297 0 : rec.define("telescope", telescope);
1298 0 : rec.define("dopb", true);
1299 0 : rec.define("isthisvp", true);
1300 0 : rec.define("realimage", functionImageName);
1301 0 : rec.define("channel", funcChannel);
1302 0 : rec.define("reffreq", nomFreq.get(uHz).getValue());
1303 0 : rec.define("minvalidfreq", loFreq.get(uHz).getValue());
1304 0 : rec.define("maxvalidfreq", hiFreq.get(uHz).getValue());
1305 0 : rval = true;
1306 0 : break;
1307 0 : case AntennaResponses::VPMAN: // the function is available in casa via the vp manager, i.e. use COMMONPB
1308 : // same as if ifield == -1
1309 0 : rec.define("name", "COMMONPB");
1310 0 : rec.define("isVP", PBMathInterface::COMMONPB);
1311 0 : rec.define("telescope", telescope);
1312 0 : rec.define("dopb", true);
1313 0 : rec.define("isthisvp", false);
1314 0 : rec.define("commonpb", telescope);
1315 0 : rec.define("dosquint", false);
1316 : {
1317 0 : String error;
1318 0 : Record tempholder;
1319 0 : QuantumHolder(Quantity(10.,"deg")).toRecord(error, tempholder);
1320 0 : rec.defineRecord("paincrement", tempholder);
1321 : }
1322 0 : rec.define("usesymmetricbeam", false);
1323 0 : rval = true;
1324 0 : break;
1325 0 : case AntennaResponses::INTERNAL: // the function is generated using the BeamCalc class
1326 : {
1327 0 : String antRayPath = functionImageName;
1328 :
1329 0 : Double refFreqHz = 0.; // the TOPO ref freq in Hz
1330 :
1331 : // determine TOPO reference frequency
1332 0 : if(fromFrameType!=MFrequency::TOPO){
1333 0 : MFrequency::Ref toFrame = MFrequency::Ref(MFrequency::TOPO, MeasFrame(obsdirection, obsPos, obstime));
1334 0 : MFrequency::Convert freqTrans(uHz, fromFrame, toFrame);
1335 0 : refFreqHz = freqTrans(mFreq.get(uHz).getValue()).get(uHz).getValue();
1336 0 : cout << "old freq " << mFreq.get(uHz).getValue() << ", new freq " << refFreqHz << endl;
1337 : }
1338 : else{
1339 0 : refFreqHz = mFreq.get(uHz).getValue();
1340 : }
1341 :
1342 0 : String beamCalcedImagePath = "./BeamCalcTmpImage_"+telescope+"_"+antennatype+"_"
1343 0 : +String::toString(refFreqHz/1E6)+"MHz";
1344 :
1345 : // calculate the beam
1346 :
1347 0 : if(!(telescope=="ALMA" || telescope=="ACA" || telescope =="OSF")){
1348 : os << LogIO::WARN << "Responses type INTERNAL provided for \"" << telescope << " in " << endl
1349 : << antRespPath << endl
1350 : << " not yet supported."
1351 0 : << LogIO::POST;
1352 0 : rval = false;
1353 : }
1354 : else{ // telescope=="ALMA" || telescope=="ACA" || telescope =="OSF"
1355 : try{
1356 : // handle preexisting beam image
1357 0 : Directory f(beamCalcedImagePath);
1358 0 : if(f.exists()){
1359 0 : os << LogIO::NORMAL << "Will re-use VP image \"" << beamCalcedImagePath << "\"" << LogIO::POST;
1360 : }
1361 : else{
1362 0 : CoordinateSystem coordsys;
1363 :
1364 : // DirectionCoordinate
1365 0 : Matrix<Double> xform(2,2);
1366 0 : xform = 0.0; xform.diagonal() = 1.0;
1367 : DirectionCoordinate dirCoords(MDirection::AZELGEO,
1368 0 : Projection(Projection::SIN),
1369 : 0.0, 0.0,
1370 0 : -0.5*C::pi/180.0/3600.0 * 5E11/refFreqHz,
1371 0 : 0.5*C::pi/180.0/3600.0 * 5E11/refFreqHz,
1372 : xform,
1373 0 : 128., 128.); // 256/2.
1374 0 : Vector<String> units(2);
1375 : //units = "deg";
1376 : //dirCoords.setWorldAxisUnits(units);
1377 :
1378 : // StokesCoordinate
1379 0 : Vector<Int> stoks(4);
1380 0 : stoks(0) = Stokes::XX;
1381 0 : stoks(1) = Stokes::XY;
1382 0 : stoks(2) = Stokes::YX;
1383 0 : stoks(3) = Stokes::YY;
1384 0 : StokesCoordinate stokesCoords(stoks);
1385 :
1386 : // SpectralCoordinate
1387 : SpectralCoordinate spectralCoords(MFrequency::TOPO,
1388 : refFreqHz,
1389 : 1.0E+3, // dummy increment
1390 : 0,
1391 0 : refFreqHz);
1392 0 : units.resize(1);
1393 0 : units = "Hz";
1394 0 : spectralCoords.setWorldAxisUnits(units);
1395 :
1396 0 : coordsys.addCoordinate(dirCoords);
1397 0 : coordsys.addCoordinate(stokesCoords);
1398 0 : coordsys.addCoordinate(spectralCoords);
1399 :
1400 0 : TiledShape ts(IPosition(4,256,256,4,1));
1401 0 : PagedImage<Complex> im(ts, coordsys, beamCalcedImagePath);
1402 0 : im.set(Complex(1.0,0.0));
1403 : // set XY and YX to zero
1404 0 : IPosition windowShape(4,im.shape()(0), im.shape()(1), 1, im.shape()(3));
1405 0 : LatticeStepper stepper(im.shape(), windowShape);
1406 0 : LatticeIterator<Complex> it(im, stepper);
1407 0 : Int planeNumber = 0;
1408 0 : for (it.reset(); !it.atEnd(); it++) {
1409 0 : if(planeNumber==1 || planeNumber==2){
1410 0 : it.woCursor() = Complex(0.,0.);
1411 : }
1412 0 : planeNumber++;
1413 : }
1414 :
1415 : // perform the ray tracing
1416 0 : ALMACalcIlluminationConvFunc almaPB;
1417 0 : Long cachesize=(HostInfo::memoryTotal(true)/8)*1024;
1418 0 : almaPB.setMaximumCacheSize(cachesize);
1419 0 : almaPB.setAntRayPath(antRayPath);
1420 0 : almaPB.applyVP(im, telescope, obstime, antennatype, antennatype,
1421 0 : MVFrequency(refFreqHz),
1422 : rotAngOffset.radian(), // the parallactic angle offset
1423 : true); // doSquint
1424 : } // endif exists
1425 0 : } catch (AipsError x) {
1426 : os << LogIO::SEVERE
1427 : << "BeamCalc failed with message " << endl
1428 : << " " << x.getMesg()
1429 0 : << LogIO::POST;
1430 0 : return false;
1431 : }
1432 :
1433 : // construct record
1434 0 : rec.define("name", "IMAGE");
1435 0 : rec.define("isVP", PBMathInterface::IMAGE);
1436 0 : rec.define("isthisvp", true);
1437 0 : rec.define("telescope", telescope);
1438 0 : rec.define("dopb", true);
1439 0 : rec.define("compleximage", beamCalcedImagePath);
1440 0 : rec.define("channel", 0);
1441 0 : rec.define("antennatype", antennatype);
1442 0 : rec.define("reffreq", refFreqHz);
1443 0 : rec.define("minvalidfreq", refFreqHz-refFreqHz/10.*5.);
1444 0 : rec.define("maxvalidfreq", refFreqHz+refFreqHz/10.*5.);
1445 0 : rval = true;
1446 : }
1447 : }
1448 0 : break;
1449 0 : case AntennaResponses::NA: // not available
1450 : default:
1451 0 : rval = false;
1452 0 : break;
1453 : } // end switch(ftype)
1454 : }
1455 : else{ // we have a PBMath response
1456 :
1457 0 : rec = vplist_p.subRecord(ifield);
1458 0 : rval = true;
1459 :
1460 : } // end if internally defined
1461 : }
1462 :
1463 207 : return rval;
1464 :
1465 : }
1466 :
1467 : // get the voltage pattern without giving observation parameters
1468 234 : Bool VPManager::getvp(Record &rec,
1469 : const String& telescope,
1470 : const String& antennatype // default: ""
1471 : ){
1472 :
1473 468 : std::lock_guard<std::recursive_mutex> locker(mutex_p);
1474 :
1475 468 : LogIO os;
1476 234 : os << LogOrigin("VPManager", "getvp2");
1477 :
1478 468 : MEpoch obstime;
1479 468 : MFrequency freq;
1480 468 : MDirection obsdirection;
1481 :
1482 234 : Int ifield = -2;
1483 234 : if(!getuserdefault(ifield,telescope,antennatype)){
1484 13 : return false;
1485 : }
1486 :
1487 221 : rec = Record();
1488 221 : Int rval=false;
1489 :
1490 221 : if(ifield==-1){ // internally defined PB, obs parameters ignored
1491 207 : rval = getvp(rec, telescope, obstime, freq, antennatype, obsdirection);
1492 : }
1493 14 : else if(ifield>=0){ // externally defined PB
1494 14 : TableRecord antRec(vplist_p.asRecord(ifield));
1495 14 : String thename = antRec.asString("name");
1496 14 : if(thename=="REFERENCE"){ // points to an AntennaResponses table
1497 : os << LogIO::SEVERE
1498 : << "Need to provide observation parameters time, frequency, and direction to access AntennaResponses table."
1499 0 : << LogIO::POST;
1500 0 : return false;
1501 : }
1502 : else{ // we have a PBMath response
1503 14 : rec = vplist_p.subRecord(ifield);
1504 14 : rval = true;
1505 : }
1506 : }// end if internally defined
1507 :
1508 221 : return rval;
1509 :
1510 : }
1511 :
1512 0 : Bool VPManager::getvps(Vector<Record> & unique_out_rec_list, // the list of unique beam records
1513 : Vector<Vector<uInt> >& beam_index, // indices to the above vector in sync with AntennaNames
1514 : const String& telescope,
1515 : const Vector<MEpoch>& inpTimeRange, // only elements 0 and 1 are used; if 1 is not present it is assumed to be inf
1516 : const Vector<MFrequency>& inpFreqRange, // must contain at least one element; beams will be provided for each element
1517 : const Vector<String>& AntennaNames, // characters 0 and 1 are used for ALMA to determine the antenna type
1518 : const MDirection& obsdirection){ // default is the Zenith
1519 :
1520 0 : LogIO os;
1521 0 : os << LogOrigin("VPManager", "getvp3");
1522 :
1523 0 : Bool rval=true;
1524 :
1525 0 : if (inpTimeRange.size()==0){
1526 0 : os << LogIO::WARN << "No observation time provided." << LogIO::POST;
1527 0 : return false;
1528 : }
1529 0 : if (inpFreqRange.size()==0){
1530 0 : os << LogIO::WARN << "No observation frequency provided." << LogIO::POST;
1531 0 : return false;
1532 : }
1533 0 : if (AntennaNames.size()==0){
1534 0 : os << LogIO::WARN << "No antenna names provided." << LogIO::POST;
1535 0 : return false;
1536 : }
1537 :
1538 0 : Bool checkTimeDep = false;
1539 0 : Unit uS("s");
1540 0 : if(inpTimeRange.size()>1 && (inpTimeRange[0].get(uS) < inpTimeRange[1].get(uS))){
1541 0 : checkTimeDep = true;
1542 : }
1543 :
1544 0 : unique_out_rec_list.resize(0);
1545 0 : beam_index.resize(0);
1546 :
1547 0 : vector<Record> recList;
1548 0 : vector<vector<uInt> > beamIndex;
1549 :
1550 0 : if(!(telescope=="ALMA" || telescope=="ACA" || telescope=="OSF")){
1551 0 : os << LogIO::WARN << "Don't know how to determine antenna type from antenna name for telescope " << telescope << LogIO::POST;
1552 0 : os << LogIO::WARN << "Will try to use whole name ..." << LogIO::POST;
1553 : }
1554 :
1555 0 : for(uInt k=0; k<AntennaNames.size(); k++){
1556 :
1557 0 : vector<uInt> index;
1558 :
1559 0 : for(uInt i=0; i<inpFreqRange.size(); i++){
1560 :
1561 0 : Record rec0;
1562 :
1563 0 : MEpoch obstime = inpTimeRange[0];
1564 0 : MFrequency freq=inpFreqRange[i];
1565 :
1566 0 : String antennatype=AntennaNames[k];
1567 0 : if(telescope=="ALMA" || telescope=="ACA" || telescope=="OSF"){
1568 0 : antennatype=antennatype(0,2);
1569 : }
1570 : //else {} // warning is given above, outside of the loop
1571 :
1572 0 : if (getvp(rec0, telescope, obstime, freq, antennatype, obsdirection)){
1573 :
1574 0 : if(checkTimeDep){ // check if there is a step in time
1575 0 : Record recOther;
1576 0 : if (getvp(recOther, telescope, inpTimeRange[1], freq, antennatype, obsdirection)){
1577 0 : if(!vpRecIsIdentical(rec0,recOther)){
1578 0 : os << LogIO::WARN << "Primary beam changes during observation time range." << LogIO::POST;
1579 0 : return false;
1580 : }
1581 : }
1582 : }
1583 :
1584 0 : Bool notPresent=true;
1585 0 : for(uInt j=0; j<recList.size(); j++){ // check if we already have this PB in our collection
1586 0 : if(vpRecIsIdentical(rec0,recList[j])){
1587 0 : index.push_back(j);
1588 0 : notPresent=false;
1589 0 : break;
1590 : }
1591 : }
1592 0 : if(notPresent){ // put new beam into the list
1593 0 : index.push_back(recList.size());
1594 0 : recList.push_back(rec0);
1595 : }
1596 :
1597 : }
1598 : else{
1599 0 : rval=false;
1600 0 : break;
1601 : }
1602 :
1603 0 : if(!rval) break;
1604 :
1605 : } // end for i
1606 :
1607 0 : beamIndex.push_back(index);
1608 :
1609 : } // end for k
1610 :
1611 0 : if(rval){
1612 0 : unique_out_rec_list=Vector<Record>(recList);
1613 0 : beam_index.resize(beamIndex.size());
1614 0 : for(uInt i=0; i<beamIndex.size(); i++){
1615 0 : beam_index[i] = Vector<uInt>(beamIndex[i]);
1616 : }
1617 : }
1618 :
1619 0 : return rval;
1620 :
1621 : }
1622 :
1623 0 : Bool VPManager::vpRecIsIdentical(const Record& rec0, const Record& rec1){
1624 :
1625 0 : Bool rval = false;
1626 0 : if( (rec0.nfields() == rec1.nfields()) &&
1627 0 : rec0.asString("name") == rec1.asString("name")
1628 : ){
1629 :
1630 0 : if(rec0.isDefined("compleximage")){
1631 0 : if(rec1.isDefined("compleximage") &&
1632 0 : rec0.asString("compleximage") == rec1.asString("compleximage") &&
1633 0 : rec0.asDouble("reffreq") == rec1.asDouble("reffreq")){
1634 0 : rval = true;
1635 : }
1636 : }
1637 0 : else if(rec0.isDefined("realimage")){
1638 0 : if(rec1.isDefined("realimage") &&
1639 0 : rec0.asString("realimage") == rec1.asString("realimage") &&
1640 0 : rec0.asDouble("reffreq") == rec1.asDouble("reffreq")){
1641 0 : rval = true;
1642 : }
1643 : }
1644 : else{
1645 0 : rval=true;
1646 : }
1647 :
1648 0 : Vector<String> fnames(5);
1649 0 : fnames[0]="dosquint";
1650 0 : fnames[1]="usesymmetricbeam";
1651 0 : fnames[2]="isthisvp";
1652 0 : fnames[3]="isVP";
1653 0 : fnames[4]="dopb";
1654 0 : for(uInt i=0; i<fnames.size(); i++){
1655 0 : if( rval && rec0.isDefined(fnames[i])){
1656 0 : rval=false;
1657 0 : if(rec1.isDefined(fnames[i]) &&
1658 0 : (rec0.asBool(fnames[i]) == rec1.asBool(fnames[i]))
1659 : ){
1660 0 : rval = true;
1661 : }
1662 : }
1663 : }
1664 : } // end if
1665 :
1666 0 : return rval;
1667 : }
1668 :
1669 :
1670 : } //# NAMESPACE CASA - END
|
