Line data Source code
1 : //# AntennaResponses.h: AntennaResponses provides access to antenna response data
2 : //# Copyright (C) 1995-1999,2000-2004
3 : //# Associated Universities, Inc. Washington DC, USA
4 : //# Copyright by ESO (in the framework of the ALMA collaboration)
5 : //#
6 : //# This library is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU Library General Public License as published by
8 : //# the Free Software Foundation; either version 2 of the License, or (at your
9 : //# option) any later version.
10 : //#
11 : //# This library 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 Library General Public
14 : //# License for more details.
15 : //#
16 : //# You should have received a copy of the GNU Library General Public License
17 : //# along with this library; if not, write to the Free Software Foundation,
18 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning CASA should be addressed as follows:
21 : //# Internet email: aips2-request@nrao.edu.
22 : //# Postal address: CASA Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : //#
27 : //#
28 : //# $Id: $
29 :
30 : //# Includes
31 : #include <casacore/casa/aips.h>
32 : #include <synthesis/MeasurementEquations/AntennaResponses.h>
33 : #include <casacore/casa/Arrays/Vector.h>
34 : #include <casacore/casa/Exceptions/Error.h>
35 : #include <casacore/casa/Quanta/MVTime.h>
36 : #include <casacore/casa/Quanta/QLogical.h>
37 : #include <casacore/casa/Quanta/Quantum.h>
38 : #include <casacore/casa/OS/Time.h>
39 : #include <casacore/casa/Logging/LogIO.h>
40 : #include <casacore/casa/Utilities/GenSort.h>
41 : #include <casacore/casa/System/AipsrcValue.h>
42 : #include <casacore/casa/BasicSL/String.h>
43 : #include <iostream>
44 : #include <casacore/measures/TableMeasures/ScalarMeasColumn.h>
45 : #include <casacore/measures/TableMeasures/ArrayMeasColumn.h>
46 : #include <casacore/measures/TableMeasures/TableMeasValueDesc.h>
47 : #include <casacore/measures/TableMeasures/TableMeasOffsetDesc.h>
48 : #include <casacore/measures/TableMeasures/TableMeasRefDesc.h>
49 : #include <casacore/measures/TableMeasures/TableMeasDesc.h>
50 : #include <casacore/measures/TableMeasures/ArrayQuantColumn.h>
51 : #include <casacore/measures/Measures/MEpoch.h>
52 : #include <casacore/measures/Measures/MDirection.h>
53 : #include <casacore/measures/Measures/MFrequency.h>
54 : #include <casacore/measures/Measures/MCFrequency.h>
55 : #include <casacore/measures/Measures/MCEpoch.h>
56 : #include <casacore/measures/Measures/MCDirection.h>
57 : #include <casacore/measures/Measures/MCPosition.h>
58 : #include <casacore/measures/Measures/MeasTable.h>
59 : #include <casacore/measures/Measures/MeasData.h>
60 : #include <casacore/measures/Measures/MeasRef.h>
61 : #include <casacore/measures/Measures/MeasFrame.h>
62 : #include <casacore/measures/Measures/MeasConvert.h>
63 : #include <casacore/tables/Tables/Table.h>
64 : #include <casacore/tables/Tables/TableDesc.h>
65 : #include <casacore/tables/Tables/ArrayColumn.h>
66 : #include <casacore/tables/Tables/ColumnDesc.h>
67 : #include <casacore/tables/Tables/ArrColDesc.h>
68 : #include <casacore/tables/Tables/SetupNewTab.h>
69 : #include <casacore/tables/Tables/ScaColDesc.h>
70 : #include <casacore/tables/Tables/ScalarColumn.h>
71 :
72 : using namespace casacore;
73 : namespace casa { //# NAMESPACE CASA - BEGIN
74 :
75 0 : AntennaResponses::AntennaResponses(const String& path){
76 0 : init(path);
77 0 : }
78 :
79 0 : Bool AntennaResponses::init(const String& path){
80 : // reset members to empty
81 :
82 0 : paths_p.resize();
83 :
84 0 : numRows_p = 0;
85 0 : ObsName_p.resize();
86 0 : StartTime_p.resize();
87 0 : AntennaType_p.resize();
88 0 : ReceiverType_p.resize();
89 0 : BeamId_p.resize();
90 0 : BeamNumber_p.resize();
91 0 : ValidCenter_p.resize();
92 0 : ValidCenterMin_p.resize();
93 0 : ValidCenterMax_p.resize();
94 0 : NumSubbands_p.resize();
95 0 : BandName_p.resize();
96 0 : SubbandMinFreq_p.resize();
97 0 : SubbandMaxFreq_p.resize();
98 0 : FuncType_p.resize();
99 0 : FuncName_p.resize();
100 0 : FuncChannel_p.resize();
101 0 : NomFreq_p.resize();
102 0 : RotAngOffset_p.resize();
103 :
104 0 : pathIndex_p.resize();
105 :
106 0 : if(path==""){
107 0 : paths_p.resize(1, "");
108 0 : return true;
109 : }
110 : else{
111 : // fill members from table on disk
112 0 : return append(path);
113 : }
114 : }
115 :
116 :
117 0 : Bool AntennaResponses::append(const String& path){
118 :
119 0 : if(isInit(path)){
120 : // Returns false if the path was already read before.
121 : //cout << "Path has been read before." << endl;
122 0 : return false;
123 : }
124 :
125 : // open table
126 0 : Table tab(path);
127 :
128 : // get number of rows
129 0 : uInt numRows = tab.nrow();
130 :
131 0 : uInt nPaths = paths_p.nelements();
132 :
133 0 : if(numRows>0){
134 : // read columns and append to vectors;
135 :
136 0 : ScalarColumn<Int> beamIdCol(tab, "BEAM_ID");
137 0 : ScalarColumn<String> nameCol(tab, "NAME");
138 0 : ScalarColumn<Int> beamNumberCol(tab, "BEAM_NUMBER");
139 0 : ROScalarMeasColumn<MEpoch> startTimeCol(tab, "START_TIME");
140 0 : ScalarColumn<String> antennaTypeCol(tab, "ANTENNA_TYPE");
141 0 : ScalarColumn<String> receiverTypeCol(tab, "RECEIVER_TYPE");
142 0 : ScalarColumn<Int> numSubbandsCol(tab, "NUM_SUBBANDS");
143 0 : ArrayColumn<String> bandNameCol(tab, "BAND_NAME");
144 0 : ROArrayQuantColumn<Double> subbandMinFreqCol(tab, "SUBBAND_MIN_FREQ");
145 0 : ROArrayQuantColumn<Double> subbandMaxFreqCol(tab, "SUBBAND_MAX_FREQ");
146 0 : ROScalarMeasColumn<MDirection> centerCol(tab, "CENTER");
147 0 : ROScalarMeasColumn<MDirection> validCenterMinCol(tab, "VALID_CENTER_MIN");
148 0 : ROScalarMeasColumn<MDirection> validCenterMaxCol(tab, "VALID_CENTER_MAX");
149 0 : ArrayColumn<Int> functionTypeCol(tab, "FUNCTION_TYPE");
150 0 : ArrayColumn<String> functionNameCol(tab, "FUNCTION_NAME");
151 0 : ArrayColumn<uInt> functionChannelCol(tab, "FUNCTION_CHANNEL");
152 0 : ROArrayQuantColumn<Double> nomFreqCol(tab, "NOMINAL_FREQ");
153 0 : ROArrayQuantColumn<Double> rotAngOffsetCol(tab, "RESPONSE_ROTATION_OFFSET");
154 :
155 0 : numRows_p += numRows;
156 :
157 0 : ObsName_p.resize(numRows_p,true);
158 0 : StartTime_p.resize(numRows_p,true);
159 0 : AntennaType_p.resize(numRows_p,true);
160 0 : ReceiverType_p.resize(numRows_p,true);
161 0 : BeamId_p.resize(numRows_p,true);
162 0 : BeamNumber_p.resize(numRows_p,true);
163 0 : ValidCenter_p.resize(numRows_p,true);
164 0 : ValidCenterMin_p.resize(numRows_p,true);
165 0 : ValidCenterMax_p.resize(numRows_p,true);
166 0 : NumSubbands_p.resize(numRows_p,true);
167 0 : BandName_p.resize(numRows_p,true);
168 0 : SubbandMinFreq_p.resize(numRows_p,true);
169 0 : SubbandMaxFreq_p.resize(numRows_p,true);
170 0 : FuncType_p.resize(numRows_p,true);
171 0 : FuncName_p.resize(numRows_p,true);
172 0 : FuncChannel_p.resize(numRows_p,true);
173 0 : NomFreq_p.resize(numRows_p,true);
174 0 : RotAngOffset_p.resize(numRows_p,true);
175 :
176 0 : pathIndex_p.resize(numRows_p,true);
177 :
178 0 : for(uInt i=0; i<numRows; i++){
179 0 : uInt j = i + numRows_p - numRows;
180 0 : BeamId_p(j) = beamIdCol(i);
181 0 : ObsName_p(j) = nameCol(i);
182 0 : BeamNumber_p(j) = beamNumberCol(i);
183 0 : StartTime_p(j) = startTimeCol(i);
184 0 : AntennaType_p(j) = antennaTypeCol(i);
185 0 : ReceiverType_p(j) = receiverTypeCol(i);
186 0 : NumSubbands_p(j) = numSubbandsCol(i);
187 0 : BandName_p(j).assign(bandNameCol(i));
188 :
189 0 : Vector<Quantity> tQ;
190 0 : tQ = subbandMinFreqCol(i);
191 0 : uInt nSubB = tQ.nelements();
192 0 : SubbandMinFreq_p(j).resize(nSubB);
193 0 : for(uInt k=0; k<nSubB; k++){
194 0 : SubbandMinFreq_p(j)(k) = MVFrequency(tQ(k));
195 : }
196 0 : tQ = subbandMaxFreqCol(i);
197 0 : SubbandMaxFreq_p(j).resize(nSubB);
198 0 : for(uInt k=0; k<nSubB; k++){
199 0 : SubbandMaxFreq_p(j)(k) = MVFrequency(tQ(k));
200 : }
201 0 : tQ = nomFreqCol(i);
202 0 : NomFreq_p(j).resize(nSubB);
203 0 : for(uInt k=0; k<nSubB; k++){
204 0 : NomFreq_p(j)(k) = MVFrequency(tQ(k));
205 : }
206 0 : tQ = rotAngOffsetCol(i);
207 0 : RotAngOffset_p(j).resize(nSubB);
208 0 : for(uInt k=0; k<nSubB; k++){
209 0 : RotAngOffset_p(j)(k) = MVAngle(tQ(k));
210 : }
211 :
212 0 : ValidCenter_p(j) = centerCol(i);
213 0 : ValidCenterMin_p(j) = validCenterMinCol(i);
214 0 : ValidCenterMax_p(j) = validCenterMaxCol(i);
215 :
216 0 : Vector<Int> tFType;
217 0 : tFType = functionTypeCol(i);
218 0 : FuncType_p(j).resize(nSubB);
219 0 : for(uInt k=0; k<nSubB; k++){
220 0 : FuncType_p(j)(k) = FuncType(tFType(k));
221 : }
222 :
223 0 : FuncName_p(j).assign(functionNameCol(i));
224 0 : FuncChannel_p(j).assign(functionChannelCol(i));
225 :
226 0 : pathIndex_p(j) = nPaths;
227 :
228 : }
229 :
230 : } // end if
231 :
232 0 : paths_p.resize(nPaths+1, true);
233 0 : String tempS = path;
234 0 : while(tempS.lastchar()=='/' && tempS.size()>1){ // don't want trailing "/"
235 0 : tempS.erase(tempS.size()-1,1);
236 : }
237 0 : paths_p(nPaths) = path;
238 :
239 0 : return true;
240 :
241 : }
242 :
243 0 : Bool AntennaResponses::isInit(){
244 :
245 0 : return (paths_p.nelements()!=0);
246 : }
247 :
248 0 : Bool AntennaResponses::isInit(const String& path){
249 0 : Bool found = false;
250 0 : for (uInt i=0; i<paths_p.nelements(); i++){
251 0 : if(paths_p(i) == path){
252 0 : found = true;
253 : }
254 : }
255 0 : return found;
256 : }
257 :
258 0 : Bool AntennaResponses::getRowAndIndex(uInt& row, uInt& subBand,
259 : const String& obsName,
260 : const MEpoch& obsTime,
261 : const MFrequency& freq,
262 : const FuncTypes& requFType,
263 : const String& antennaType,
264 : const MDirection& center,
265 : const String& receiverType,
266 : const Int& beamNumber){
267 0 : Bool rval = false;
268 0 : Unit uS("s");
269 0 : Unit uHz("Hz");
270 :
271 : // for the combination INTERNAL + 0 freq, return the row and index for the first freq found
272 0 : Bool matchFreq = !((freq.get(uHz).getValue() == 0.) && requFType==INTERNAL);
273 :
274 : // calculate azimuth, elevation, and topo frequency
275 : // first, get the reference frame
276 0 : MPosition mp;
277 0 : if (!MeasTable::Observatory(mp,obsName)) {
278 : // unknown observatory
279 0 : LogIO os(LogOrigin("AntennaResponses",
280 0 : String("getRowAndIndex"),
281 0 : WHERE));
282 0 : os << LogIO::NORMAL << String("Unknown observatory ") << obsName
283 0 : << LogIO::POST;
284 0 : return false;
285 : }
286 0 : mp=MPosition::Convert(mp, MPosition::ITRF)();
287 0 : MeasFrame frame(mp, obsTime);
288 0 : Vector<Double> azel = MDirection::Convert(center,
289 0 : MDirection::Ref(MDirection::AZEL,
290 : frame)
291 0 : )().getAngle("deg").getValue();
292 :
293 0 : Quantity f(0., uHz);
294 0 : if(matchFreq){
295 0 : f = MFrequency::Convert(freq, MFrequency::TOPO)().get(uHz);
296 : }
297 :
298 :
299 : // loop over rows
300 : uInt i,j;
301 0 : vector<uInt> rowV, subBandV;
302 0 : vector<Quantity> timeV;
303 0 : for(i=0; i<numRows_p; i++){
304 0 : if(ObsName_p(i) == obsName
305 0 : && StartTime_p(i).get(uS) <= obsTime.get(uS)
306 0 : && AntennaType_p(i) == antennaType
307 0 : && ReceiverType_p(i) == receiverType
308 0 : && BeamNumber_p(i) == beamNumber
309 : ){
310 : // remains to test center and freq and functype
311 : // first freq and functype
312 0 : Bool found = false;
313 0 : for(j=0; j<NumSubbands_p(i); j++){
314 : //cout << "f " << f << " min " << SubbandMinFreq_p(i)(j).get()
315 : // << " max " << SubbandMaxFreq_p(i)(j).get() << endl;
316 0 : if( (FuncType_p(i)(j) == requFType
317 0 : || requFType == AntennaResponses::ANY)
318 0 : && (
319 0 : !matchFreq // if matchFreq is false, any frequency is accepted
320 0 : || (SubbandMinFreq_p(i)(j).get() <= f && f <= SubbandMaxFreq_p(i)(j).get())
321 : )
322 : ){
323 0 : found = true;
324 0 : break;
325 : }
326 : }
327 0 : if(found){ // now test center
328 0 : Vector<Double> azelMin = MDirection::Convert(ValidCenterMin_p(i),
329 0 : MDirection::Ref(MDirection::AZEL,
330 : frame)
331 0 : )().getAngle("deg").getValue();
332 0 : Vector<Double> azelMax = MDirection::Convert(ValidCenterMax_p(i),
333 0 : MDirection::Ref(MDirection::AZEL,
334 : frame)
335 0 : )().getAngle("deg").getValue();
336 :
337 : // need to accomodate the ambiguity of the AZ
338 0 : Double modAz = fmod(azel(0) - azelMin(0),360.);
339 0 : Double modAzMax = fmod(azelMax(0) - azelMin(0), 360.);
340 :
341 : // cout << " i, j, " << i << ", " << j << " az min actual max" << azelMin(0) << " "
342 : // << azel(0) << " " << azelMax(0) << endl;
343 : // cout << " i, j, " << i << ", " << j << " cannonised az actual max" << modAz << " " << modAzMax << endl;
344 : // cout << " i, j, " << i << ", " << j << " el min actual max" << azelMin(1) << " " << azel(1) << " " << azelMax(1) << endl;
345 :
346 0 : if((fabs(azelMin(0)-azelMax(0))<1e-5 // all AZ are valid (accomodate numerical problems at 360 deg)
347 0 : || ( 0. <= modAz && modAz <= modAzMax))
348 0 : && azelMin(1) <= azel(1) && azel(1) <= azelMax(1)){
349 : // memorize the applicable row, sub band, and time
350 0 : rval = true;
351 0 : rowV.push_back(i);
352 0 : subBandV.push_back(j);
353 0 : timeV.push_back(StartTime_p(i).get(uS));
354 : }
355 : }
356 : }
357 : } // end for
358 :
359 0 : if(rval){
360 0 : Vector<uInt> sortIndex;
361 0 : GenSortIndirect<Quantity>::sort(sortIndex, Vector<Quantity>(timeV));
362 : // take the latest row
363 0 : row = rowV[sortIndex(sortIndex.nelements()-1)];
364 0 : subBand = subBandV[sortIndex(sortIndex.nelements()-1)];
365 : }
366 0 : return rval;
367 :
368 : }
369 :
370 0 : Bool AntennaResponses::getRowAndIndex(uInt& row, uInt& subBand,
371 : const String& obsName,
372 : const Int& beamId,
373 : const MFrequency& freq){
374 :
375 0 : Bool rval = false;
376 0 : Unit uHz("Hz");
377 :
378 : // calculate topo frequency
379 0 : Quantity f = MFrequency::Convert(freq, MFrequency::TOPO)().get(uHz);
380 :
381 : // loop over rows
382 : uInt i,j;
383 0 : for(i=0; i<numRows_p; i++){
384 0 : if(ObsName_p(i) == obsName
385 0 : && BeamNumber_p(i) == beamId
386 : ){
387 0 : for(j=0; j<NumSubbands_p(i); j++){
388 0 : if(SubbandMinFreq_p(i)(j).get() <= f
389 0 : && f <= SubbandMaxFreq_p(i)(j).get()
390 : ){
391 0 : rval = true;
392 0 : break;
393 : }
394 : }
395 0 : if(rval){
396 0 : break;
397 : }
398 : }
399 : } // end for
400 :
401 0 : if(rval){
402 0 : row = i;
403 0 : subBand = j;
404 : }
405 0 : return rval;
406 :
407 : }
408 :
409 :
410 0 : Bool AntennaResponses::getImageName(String& functionImageName, // the path to the image
411 : uInt& functionChannel, // the channel to use
412 : MFrequency& nomFreq, // nominal frequency of the image (in the given channel)
413 : FuncTypes& fType, // the function type of the image
414 : MVAngle& rotAngOffset, // response rotation angle offset
415 : const String& obsName, // (the observatory name, e.g. "ALMA" or "ACA")
416 : const MEpoch& obsTime,
417 : const MFrequency& freq,
418 : const FuncTypes& requFType, // the requested function type
419 : const String& antennaType,
420 : const MDirection& center,
421 : const String& receiverType,
422 : const Int& beamNumber){
423 :
424 : uInt row, subBand;
425 :
426 0 : if(!getRowAndIndex(row, subBand,
427 : obsName, obsTime, freq,
428 : requFType, antennaType,
429 : center, receiverType,
430 : beamNumber)){
431 0 : return false;
432 : }
433 : else{
434 0 : functionImageName = FuncName_p(row)(subBand);
435 0 : if(functionImageName.firstchar()!='/'){ // need to prepend the path
436 0 : String tempS = paths_p(pathIndex_p(row));
437 0 : string::size_type p = tempS.find_last_of('/',tempS.size());
438 0 : functionImageName = tempS.substr(0,p) + "/" + functionImageName;
439 : }
440 0 : functionChannel = FuncChannel_p(row)(subBand);
441 0 : nomFreq = MFrequency(NomFreq_p(row)(subBand),MFrequency::TOPO);
442 0 : fType = FuncType_p(row)(subBand);
443 0 : rotAngOffset = RotAngOffset_p(row)(subBand);
444 0 : return true;
445 : }
446 :
447 : }
448 :
449 : // overloaded method with additional output hi and lo ends of validity range
450 0 : Bool AntennaResponses::getImageName(String& functionImageName, // the path to the image
451 : uInt& functionChannel, // the channel to use
452 : MFrequency& nomFreq, // nominal frequency of the image (in the given channel)
453 : MFrequency& loFreq, // lower end of the frequency validity range
454 : MFrequency& hiFreq, // upper end of the frequency validity range
455 : FuncTypes& fType, // the function type of the image
456 : MVAngle& rotAngOffset, // response rotation angle offset
457 : const String& obsName, // (the observatory name, e.g. "ALMA" or "ACA")
458 : const MEpoch& obsTime,
459 : const MFrequency& freq,
460 : const FuncTypes& requFType, // the requested function type
461 : const String& antennaType,
462 : const MDirection& center,
463 : const String& receiverType,
464 : const Int& beamNumber){
465 :
466 : uInt row, subBand;
467 :
468 0 : if(!getRowAndIndex(row, subBand,
469 : obsName, obsTime, freq,
470 : requFType, antennaType,
471 : center, receiverType,
472 : beamNumber)){
473 0 : return false;
474 : }
475 : else{
476 0 : functionImageName = FuncName_p(row)(subBand);
477 0 : if(functionImageName.firstchar()!='/'){ // need to prepend the path
478 0 : String tempS = paths_p(pathIndex_p(row));
479 0 : string::size_type p = tempS.find_last_of('/',tempS.size());
480 0 : functionImageName = tempS.substr(0,p) + "/" + functionImageName;
481 : }
482 0 : functionChannel = FuncChannel_p(row)(subBand);
483 0 : nomFreq = MFrequency(NomFreq_p(row)(subBand),MFrequency::TOPO);
484 0 : loFreq = MFrequency(SubbandMinFreq_p(row)(subBand),MFrequency::TOPO);
485 0 : hiFreq = MFrequency(SubbandMaxFreq_p(row)(subBand),MFrequency::TOPO);
486 0 : fType = FuncType_p(row)(subBand);
487 0 : rotAngOffset = RotAngOffset_p(row)(subBand);
488 0 : return true;
489 : }
490 :
491 :
492 : }
493 :
494 :
495 : // overloaded method: similar to previous but using beamId
496 : // (instead of obs. time, ant. type, rec. type, and center)
497 0 : Bool AntennaResponses::getImageName(String& functionImageName,
498 : uInt& functionChannel,
499 : MFrequency& nomFreq, // nominal frequency of the image
500 : FuncTypes& fType, // the function type of the image
501 : MVAngle& rotAngOffset, // response rotation angle offset
502 : const String& obsName, // (the observatory name, e.g. "ALMA" or "ACA")
503 : const Int& beamId,
504 : const MFrequency& freq){
505 : uInt row, subBand;
506 :
507 0 : if(!getRowAndIndex(row, subBand,
508 : obsName, beamId, freq)){
509 0 : return false;
510 : }
511 : else{
512 : // cout << "row " << row << " subband " << subBand << endl;
513 : // cout << "numrows " << numRows_p << endl;
514 : // for(uInt i=0; i<FuncName_p.nelements(); i++){
515 : // for(uInt j=0; j< FuncName_p(i).nelements(); j++){
516 : // cout << "i, j " << i << ", " << j << " " << FuncName_p(i)(j) << endl;
517 : // }
518 : // }
519 0 : functionImageName = FuncName_p(row)(subBand);
520 0 : if(functionImageName.firstchar()!='/'){ // need to prepend the path
521 0 : String tempS = paths_p(pathIndex_p(row));
522 0 : string::size_type p = tempS.find_last_of('/',tempS.size());
523 0 : functionImageName = tempS.substr(0,p) + "/" + functionImageName;
524 : }
525 0 : functionChannel = FuncChannel_p(row)(subBand);
526 0 : nomFreq = MFrequency(NomFreq_p(row)(subBand),MFrequency::TOPO);
527 0 : fType = FuncType_p(row)(subBand);
528 0 : rotAngOffset = RotAngOffset_p(row)(subBand);
529 0 : return true;
530 : }
531 : }
532 :
533 0 : Bool AntennaResponses::getAntennaTypes(Vector<String>& antTypes,
534 : const String& obsName, // (the observatory name, e.g. "ALMA" or "ACA")
535 : const MEpoch& obsTime,
536 : const MFrequency& freq,
537 : const FuncTypes& requFType,
538 : const MDirection& center,
539 : const String& receiverType,
540 : const Int& beamNumber){
541 0 : Bool rval = false;
542 0 : Unit uS("s");
543 0 : Unit uHz("Hz");
544 :
545 0 : antTypes.resize(0);
546 :
547 : // for the combination INTERNAL + 0 freq, return the row and index for the first freq found
548 0 : Bool matchFreq = !((freq.get(uHz).getValue() == 0.) && requFType==INTERNAL);
549 :
550 : // calculate azimuth, elevation, and topo frequency
551 : // first, get the reference frame
552 0 : MPosition mp;
553 0 : if (!MeasTable::Observatory(mp,obsName)) {
554 : // unknown observatory
555 0 : LogIO os(LogOrigin("AntennaResponses",
556 0 : String("getAntennaTypes"),
557 0 : WHERE));
558 0 : os << LogIO::NORMAL << String("Unknown observatory ") << obsName
559 0 : << LogIO::POST;
560 0 : return false;
561 : }
562 0 : mp=MPosition::Convert(mp, MPosition::ITRF)();
563 0 : MeasFrame frame(mp, obsTime);
564 0 : Vector<Double> azel = MDirection::Convert(center,
565 0 : MDirection::Ref(MDirection::AZEL,
566 : frame)
567 0 : )().getAngle("deg").getValue();
568 :
569 0 : Quantity f(0., uHz);
570 0 : if(matchFreq){
571 0 : f = MFrequency::Convert(freq, MFrequency::TOPO)().get(uHz);
572 : }
573 :
574 :
575 : // loop over rows
576 : uInt i,j;
577 0 : std::map<String, Int > antTypeMap;
578 0 : for(i=0; i<numRows_p; i++){
579 0 : if(ObsName_p(i) == obsName
580 0 : && StartTime_p(i).get(uS) <= obsTime.get(uS)
581 0 : && ReceiverType_p(i) == receiverType
582 0 : && BeamNumber_p(i) == beamNumber
583 : ){
584 : // remains to test center and freq and functype
585 : // first freq and functype
586 0 : Bool found = false;
587 0 : for(j=0; j<NumSubbands_p(i); j++){
588 : //cout << "f " << f << " min " << SubbandMinFreq_p(i)(j).get()
589 : // << " max " << SubbandMaxFreq_p(i)(j).get() << endl;
590 0 : if( (FuncType_p(i)(j) == requFType
591 0 : || requFType == AntennaResponses::ANY)
592 0 : && (
593 0 : !matchFreq // if matchFreq is false, any frequency is accepted
594 0 : || (SubbandMinFreq_p(i)(j).get() <= f && f <= SubbandMaxFreq_p(i)(j).get())
595 : )
596 : ){
597 0 : found = true;
598 0 : break;
599 : }
600 : }
601 0 : if(found){ // now test center
602 0 : Vector<Double> azelMin = MDirection::Convert(ValidCenterMin_p(i),
603 0 : MDirection::Ref(MDirection::AZEL,
604 : frame)
605 0 : )().getAngle("deg").getValue();
606 0 : Vector<Double> azelMax = MDirection::Convert(ValidCenterMax_p(i),
607 0 : MDirection::Ref(MDirection::AZEL,
608 : frame)
609 0 : )().getAngle("deg").getValue();
610 :
611 : // need to accomodate the ambiguity of the AZ
612 0 : Double modAz = fmod(azel(0) - azelMin(0),360.);
613 0 : Double modAzMax = fmod(azelMax(0) - azelMin(0), 360.);
614 :
615 0 : if((fabs(azelMin(0)-azelMax(0))<1e-5 // all AZ are valid (accomodate numerical problems at 360 deg)
616 0 : || ( 0. <= modAz && modAz <= modAzMax))
617 0 : && azelMin(1) <= azel(1) && azel(1) <= azelMax(1)){
618 : // memorize the antenna type
619 0 : if( antTypeMap.find(AntennaType_p(i)) == antTypeMap.end( ) ){
620 0 : rval = true;
621 0 : antTypeMap.insert(std::pair<String, Int >(AntennaType_p(i),1));
622 : }
623 : }
624 : }
625 : }
626 : } // end for
627 :
628 0 : if(rval){
629 0 : antTypes.resize(antTypeMap.size( ));
630 0 : uInt count = 0;
631 0 : for( auto iter = antTypeMap.begin( ); iter != antTypeMap.end( ); ++iter, ++count ){
632 0 : antTypes(count) = iter->first;
633 : }
634 : }
635 :
636 0 : return rval;
637 :
638 : }
639 :
640 :
641 0 : Bool AntennaResponses::putRow(uInt& row,
642 : const String& obsName,
643 : const Int& beamId,
644 : const Vector<String>& bandName,
645 : const Vector<MVFrequency>& subbandMinFreq,
646 : const Vector<MVFrequency>& subbandMaxFreq,
647 : const Vector<FuncTypes>& funcType,
648 : const Vector<String>& funcName,
649 : const Vector<uInt>& funcChannel,
650 : const Vector<MVFrequency>& nomFreq,
651 : const Vector<MVAngle>& rotAngOffset,
652 : const String& antennaType,
653 : const MEpoch& startTime,
654 : const MDirection& center,
655 : const MDirection& validCenterMin,
656 : const MDirection& validCenterMax,
657 : const String& receiverType,
658 : const Int& beamNumber){
659 : // Put the given row into the present antenna reponses table (in memory).
660 :
661 : // Returns false, if the table was not initialised or the given data was
662 : // not consistent.
663 0 : if(paths_p.nelements()==0){
664 0 : LogIO os(LogOrigin("AntennaResponses",
665 0 : String("putRow"),
666 0 : WHERE));
667 0 : os << LogIO::NORMAL << String("Table not initialized.") << obsName
668 0 : << LogIO::POST;
669 0 : return false;
670 : }
671 : // Consistency checks:
672 : // - all vectors have same dimension which is then used to set numSubbands
673 0 : uInt tNumSubbands = bandName.nelements();
674 0 : if(!(
675 0 : tNumSubbands == subbandMinFreq.nelements() &&
676 0 : tNumSubbands == subbandMaxFreq.nelements() &&
677 0 : tNumSubbands == funcType.nelements() &&
678 0 : tNumSubbands == funcName.nelements() &&
679 0 : tNumSubbands == funcChannel.nelements() &&
680 0 : tNumSubbands == nomFreq.nelements() &&
681 0 : tNumSubbands == rotAngOffset.nelements())
682 : ){
683 0 : LogIO os(LogOrigin("AntennaResponses", String("putRow"), WHERE));
684 0 : os << LogIO::NORMAL << String("Inconsistent vector dimensions.") << obsName
685 0 : << LogIO::POST;
686 0 : return false;
687 : }
688 : // - beamId is unique for the given observatory
689 0 : Bool isUnique = true;
690 0 : for(uInt i=0; i<numRows_p; i++){
691 0 : if(ObsName_p(i)==obsName && BeamId_p(i)==beamId && row!=i){
692 0 : isUnique = false;
693 0 : break;
694 : }
695 : }
696 0 : if(!isUnique){
697 0 : LogIO os(LogOrigin("AntennaResponses", String("putRow"), WHERE));
698 : os << LogIO::WARN << "Beam id " << beamId << " not unique."
699 0 : << LogIO::POST;
700 0 : return false;
701 : }
702 : // - center, validCenterMin, and validCenterMax have the same MDirection type
703 0 : String dirRef = center.getRefString();
704 0 : if(!(dirRef == validCenterMin.getRefString() &&
705 0 : dirRef == validCenterMax.getRefString())
706 : ){
707 0 : LogIO os(LogOrigin("AntennaResponses",
708 0 : String("putRow"),
709 0 : WHERE));
710 : os << LogIO::WARN << "Inconsistent direction type."
711 0 : << LogIO::POST;
712 0 : return false;
713 : }
714 :
715 0 : uInt theRow = 0;
716 :
717 : // If the row exists at the position given by uInt row, it is overwritten.
718 0 : if(row<numRows_p){
719 0 : theRow = row;
720 : }
721 : else{ // If it doesn't exist, the table is resized by one in memory and the new
722 : // row is added at the last position. The variable "row" then contains the
723 : // actual row that was filled.
724 0 : theRow = row = numRows_p;
725 0 : numRows_p++;
726 0 : ObsName_p.resize(numRows_p,true);
727 0 : StartTime_p.resize(numRows_p,true);
728 0 : AntennaType_p.resize(numRows_p,true);
729 0 : ReceiverType_p.resize(numRows_p,true);
730 0 : BeamId_p.resize(numRows_p,true);
731 0 : BeamNumber_p.resize(numRows_p,true);
732 0 : ValidCenter_p.resize(numRows_p,true);
733 0 : ValidCenterMin_p.resize(numRows_p,true);
734 0 : ValidCenterMax_p.resize(numRows_p,true);
735 0 : NumSubbands_p.resize(numRows_p,true);
736 0 : BandName_p.resize(numRows_p,true);
737 0 : SubbandMinFreq_p.resize(numRows_p,true);
738 0 : SubbandMaxFreq_p.resize(numRows_p,true);
739 0 : FuncType_p.resize(numRows_p,true);
740 0 : FuncName_p.resize(numRows_p,true);
741 0 : FuncChannel_p.resize(numRows_p,true);
742 0 : NomFreq_p.resize(numRows_p,true);
743 0 : RotAngOffset_p.resize(numRows_p,true);
744 : }
745 :
746 0 : ObsName_p(theRow) = obsName;
747 0 : StartTime_p(theRow) = startTime;
748 0 : AntennaType_p(theRow) = antennaType;
749 0 : ReceiverType_p(theRow) = receiverType;
750 0 : BeamId_p(theRow) = beamId;
751 0 : BeamNumber_p(theRow) = beamNumber;
752 0 : ValidCenter_p(theRow) = center;
753 0 : ValidCenterMin_p(theRow) = validCenterMin;
754 0 : ValidCenterMax_p(theRow) = validCenterMax;
755 0 : NumSubbands_p(theRow) = tNumSubbands;
756 0 : BandName_p(theRow).assign(bandName);
757 0 : SubbandMinFreq_p(theRow).assign(subbandMinFreq);
758 0 : SubbandMaxFreq_p(theRow).assign(subbandMaxFreq);
759 0 : FuncType_p(theRow).assign(funcType);
760 0 : FuncName_p(theRow).assign(funcName);
761 0 : FuncChannel_p(theRow).assign(funcChannel);
762 0 : NomFreq_p(theRow).assign(nomFreq);
763 0 : RotAngOffset_p(theRow).assign(rotAngOffset);
764 :
765 0 : return true;
766 :
767 : }
768 :
769 :
770 0 : void AntennaResponses::create(const String& path){
771 :
772 : // set up table description
773 :
774 0 : TableDesc tD("AntennaResponsesDesc", TableDesc::New);
775 0 : tD.comment() = "antenna responses table";
776 :
777 0 : tD.addColumn(ScalarColumnDesc<Int> ("BEAM_ID", "unique for the given observatory name"));
778 0 : tD.addColumn(ScalarColumnDesc<String> ("NAME", "name of the observatory as in the Observatories table"));
779 0 : tD.addColumn(ScalarColumnDesc<Int> ("BEAM_NUMBER", "for observataories which support several simultaneous beams, zero-based"));
780 0 : tD.addColumn(ScalarColumnDesc<Double> ("START_TIME", "the time from which onwards this table row is valid, measure fixed to UTC"));
781 0 : tD.addColumn(ScalarColumnDesc<String> ("ANTENNA_TYPE", "for heterogeneous arrays: indication of the antenna type"));
782 0 : tD.addColumn(ScalarColumnDesc<String> ("RECEIVER_TYPE", "permits multiple receivers per band"));
783 0 : tD.addColumn(ScalarColumnDesc<Int> ("NUM_SUBBANDS", "number of elements in the array columns in this table"));
784 0 : tD.addColumn(ArrayColumnDesc<String> ("BAND_NAME", "name of the frequency band"));
785 0 : tD.addColumn(ArrayColumnDesc<Double> ("SUBBAND_MIN_FREQ", "minimum frequency of the subband in the observatory frame"));
786 0 : tD.addColumn(ArrayColumnDesc<Double> ("SUBBAND_MAX_FREQ", "maximum frequency of the subband in the observatory frame"));
787 0 : tD.addColumn(ArrayColumnDesc<Double> ("CENTER", "the nominal center sky position where this row is valid"));
788 0 : tD.addColumn(ScalarColumnDesc<Int> ("CENTER_REF", ColumnDesc::Direct));
789 0 : tD.addColumn(ArrayColumnDesc<Double> ("VALID_CENTER_MIN", "sky position validity range min values"));
790 0 : tD.addColumn(ArrayColumnDesc<Double> ("VALID_CENTER_MAX", "sky position validity range max values"));
791 0 : tD.addColumn(ArrayColumnDesc<Int> ("FUNCTION_TYPE"));
792 0 : tD.addColumn(ArrayColumnDesc<String> ("FUNCTION_NAME"));
793 0 : tD.addColumn(ArrayColumnDesc<uInt> ("FUNCTION_CHANNEL"));
794 0 : tD.addColumn(ArrayColumnDesc<Double> ("NOMINAL_FREQ"));
795 0 : tD.addColumn(ArrayColumnDesc<Double> ("RESPONSE_ROTATION_OFFSET"));
796 :
797 : // Add TableMeasures information for designated Measures/Quanta columns
798 :
799 0 : TableMeasValueDesc timeMeasVal(tD, "START_TIME");
800 0 : TableMeasRefDesc timeMeasRef(MEpoch::DEFAULT);
801 0 : TableMeasDesc<MEpoch> timeMeasCol(timeMeasVal, timeMeasRef);
802 0 : timeMeasCol.write(tD);
803 :
804 0 : TableQuantumDesc timeQuantDesc(tD, "START_TIME", Unit ("s"));
805 0 : timeQuantDesc.write(tD);
806 :
807 0 : TableQuantumDesc freqQuantDesc(tD, "SUBBAND_MIN_FREQ", Unit ("Hz"));
808 0 : freqQuantDesc.write (tD);
809 0 : TableQuantumDesc freqQuantDesc2(tD, "SUBBAND_MAX_FREQ", Unit ("Hz"));
810 0 : freqQuantDesc2.write (tD);
811 0 : TableQuantumDesc freqQuantDesc3(tD, "NOMINAL_FREQ", Unit ("Hz"));
812 0 : freqQuantDesc3.write (tD);
813 0 : TableQuantumDesc angQuantDesc(tD, "RESPONSE_ROTATION_OFFSET", Unit ("deg"));
814 0 : angQuantDesc.write (tD);
815 :
816 0 : TableMeasValueDesc refDirMeasVal (tD, "CENTER");
817 0 : TableMeasRefDesc refDirMeasRef (tD, "CENTER_REF");
818 0 : TableMeasDesc<MDirection> refDirMeasCol (refDirMeasVal, refDirMeasRef);
819 0 : refDirMeasCol.write(tD);
820 0 : TableMeasValueDesc refDirMeasValMin (tD, "VALID_CENTER_MIN");
821 0 : TableMeasDesc<MDirection> refDirMeasColMin (refDirMeasValMin, refDirMeasRef);
822 0 : refDirMeasColMin.write(tD);
823 0 : TableMeasValueDesc refDirMeasValMax (tD, "VALID_CENTER_MAX");
824 0 : TableMeasDesc<MDirection> refDirMeasColMax (refDirMeasValMax, refDirMeasRef);
825 0 : refDirMeasColMax.write(tD);
826 :
827 : // create the table
828 0 : SetupNewTable newtab (path, tD, Table::New);
829 0 : Table tab(newtab, numRows_p);
830 :
831 0 : if(numRows_p>0){
832 : // fill the table
833 :
834 0 : ScalarColumn<Int> beamIdCol(tab, "BEAM_ID");
835 0 : ScalarColumn<String> nameCol(tab, "NAME");
836 0 : ScalarColumn<Int> beamNumberCol(tab, "BEAM_NUMBER");
837 0 : ScalarMeasColumn<MEpoch> startTimeCol(tab, "START_TIME");
838 0 : ScalarColumn<String> antennaTypeCol(tab, "ANTENNA_TYPE");
839 0 : ScalarColumn<String> receiverTypeCol(tab, "RECEIVER_TYPE");
840 0 : ScalarColumn<Int> numSubbandsCol(tab, "NUM_SUBBANDS");
841 0 : ArrayColumn<String> bandNameCol(tab, "BAND_NAME");
842 0 : ArrayQuantColumn<Double> subbandMinFreqCol(tab, "SUBBAND_MIN_FREQ");
843 0 : ArrayQuantColumn<Double> subbandMaxFreqCol(tab, "SUBBAND_MAX_FREQ");
844 0 : ScalarMeasColumn<MDirection> centerCol(tab, "CENTER");
845 0 : ScalarMeasColumn<MDirection> validCenterMinCol(tab, "VALID_CENTER_MIN");
846 0 : ScalarMeasColumn<MDirection> validCenterMaxCol(tab, "VALID_CENTER_MAX");
847 0 : ArrayColumn<Int> functionTypeCol(tab, "FUNCTION_TYPE");
848 0 : ArrayColumn<String> functionNameCol(tab, "FUNCTION_NAME");
849 0 : ArrayColumn<uInt> functionChannelCol(tab, "FUNCTION_CHANNEL");
850 0 : ArrayQuantColumn<Double> nomFreqCol(tab, "NOMINAL_FREQ");
851 0 : ArrayQuantColumn<Double> rotAngOffsetCol(tab, "RESPONSE_ROTATION_OFFSET");
852 :
853 0 : for(uInt i=0; i<numRows_p; i++){
854 0 : beamIdCol.put(i, BeamId_p(i));
855 0 : nameCol.put(i, ObsName_p(i));
856 0 : beamNumberCol.put(i, BeamNumber_p(i));
857 0 : startTimeCol.put(i, StartTime_p(i));
858 0 : antennaTypeCol.put(i, AntennaType_p(i));
859 0 : receiverTypeCol.put(i, ReceiverType_p(i));
860 0 : numSubbandsCol.put(i, NumSubbands_p(i));
861 0 : bandNameCol.put(i, BandName_p(i));
862 :
863 0 : Vector<Quantity> bMF(SubbandMinFreq_p(i).nelements());
864 0 : for(uInt k=0; k<bMF.nelements(); k++){
865 0 : bMF(k) = (SubbandMinFreq_p(i)(k)).get(); // convert MVFrequency to Quantity in Hz
866 : }
867 0 : subbandMinFreqCol.put(i, bMF);
868 :
869 0 : for(uInt k=0; k<bMF.nelements(); k++){
870 0 : bMF(k) = (SubbandMaxFreq_p(i)(k)).get(); // convert MVFrequency to Quantity in Hz
871 : }
872 0 : subbandMaxFreqCol.put(i, bMF);
873 :
874 0 : for(uInt k=0; k<bMF.nelements(); k++){
875 0 : bMF(k) = (NomFreq_p(i)(k)).get(); // convert MVFrequency to Quantity in Hz
876 : }
877 0 : nomFreqCol.put(i, bMF);
878 :
879 0 : for(uInt k=0; k<bMF.nelements(); k++){
880 0 : bMF(k) = (RotAngOffset_p(i)(k)).get(Unit("deg")); // convert MVAngle to Quantity in deg
881 : }
882 0 : rotAngOffsetCol.put(i, bMF);
883 :
884 0 : centerCol.put(i, ValidCenter_p(i));
885 0 : validCenterMinCol.put(i, ValidCenterMin_p(i));
886 0 : validCenterMaxCol.put(i, ValidCenterMax_p(i));
887 :
888 0 : Vector<Int> iFT(FuncType_p(i).nelements());
889 0 : for(uInt k=0; k<iFT.nelements(); k++){
890 0 : iFT(k) = static_cast<Int>(FuncType_p(i)(k));
891 : }
892 0 : functionTypeCol.put(i, iFT);
893 :
894 0 : functionNameCol.put(i, FuncName_p(i));
895 0 : functionChannelCol.put(i, FuncChannel_p(i));
896 : }
897 :
898 : }
899 :
900 0 : return;
901 :
902 : }
903 :
904 0 : AntennaResponses::FuncTypes AntennaResponses::FuncType(Int i){
905 0 : if(-1 < i && i < static_cast<Int>(AntennaResponses::N_FuncTypes) ){
906 0 : return static_cast<AntennaResponses::FuncTypes>(i);
907 : }
908 : else{
909 0 : return AntennaResponses::INVALID;
910 : }
911 : }
912 :
913 0 : AntennaResponses::FuncTypes AntennaResponses::FuncType(const String& sftyp){
914 :
915 0 : if(sftyp=="NA") return AntennaResponses::NA;
916 0 : if(sftyp=="AIF") return AntennaResponses::AIF;
917 0 : if(sftyp=="EFP") return AntennaResponses::EFP;
918 0 : if(sftyp=="VP") return AntennaResponses::VP;
919 0 : if(sftyp=="VPMAN") return AntennaResponses::VPMAN;
920 0 : if(sftyp=="INTERNAL") return AntennaResponses::INTERNAL;
921 0 : return AntennaResponses::INVALID;
922 :
923 : }
924 :
925 0 : Bool AntennaResponses::getBandName(String& bandName,
926 : const String& obsName,
927 : const MVFrequency& freq){
928 : // brute force search
929 0 : Quantity f = freq.get();
930 0 : bandName = "";
931 0 : Bool rval = false;
932 :
933 : uInt i, j;
934 :
935 0 : for(i=0; i<numRows_p; i++){
936 0 : if(obsName == ObsName_p(i)){
937 0 : for(j=0; j<NumSubbands_p(i); j++){
938 0 : if(SubbandMinFreq_p(i)(j).get() <= f
939 0 : && f <= SubbandMaxFreq_p(i)(j).get()){
940 0 : rval = true;
941 0 : break;
942 : }
943 : }
944 0 : if(rval){
945 0 : break;
946 : }
947 : }
948 : }
949 0 : if(rval){
950 0 : bandName = BandName_p(i)(j);
951 : }
952 0 : return rval;
953 :
954 : }
955 :
956 : } //# NAMESPACE CASA - END
957 :
|
