Line data Source code
1 : //# PBMath.cc: Implementation for PBMath
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2002,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be adressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //#
27 : //# $Id$
28 :
29 : #include <casacore/casa/aips.h>
30 : #include <casacore/casa/BasicSL/Complex.h>
31 : #include <synthesis/TransformMachines/PBMath.h>
32 : #include <synthesis/TransformMachines/PBMath1DAiry.h>
33 : #include <synthesis/TransformMachines/PBMath1DCosPoly.h>
34 : #include <synthesis/TransformMachines/PBMath1DGauss.h>
35 : #include <synthesis/TransformMachines/PBMath1DNumeric.h>
36 : #include <synthesis/TransformMachines/PBMath1DPoly.h>
37 : #include <synthesis/TransformMachines/PBMath1DEVLA.h>
38 :
39 : #include <casacore/measures/Measures/Stokes.h>
40 : #include <casacore/casa/BasicSL/Constants.h>
41 :
42 : #include <components/ComponentModels/Flux.h>
43 : #include <components/ComponentModels/ComponentShape.h>
44 :
45 : #include <casacore/images/Images/PagedImage.h>
46 : #include <casacore/images/Images/TempImage.h>
47 : #include <casacore/images/Images/ImageInterface.h>
48 : #include <casacore/images/Regions/ImageRegion.h>
49 :
50 : #include <casacore/lattices/Lattices/LatticeStepper.h>
51 : #include <casacore/lattices/Lattices/LatticeIterator.h>
52 : #include <casacore/lattices/LEL/LatticeExpr.h>
53 :
54 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
55 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
56 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
57 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
58 : #include <casacore/coordinates/Coordinates/Projection.h>
59 : #include <casacore/casa/Logging/LogIO.h>
60 :
61 : #include <casacore/casa/Utilities/Assert.h>
62 : #include <casacore/casa/Containers/RecordInterface.h>
63 : #include <casacore/casa/Quanta/QuantumHolder.h>
64 : #include <casacore/measures/Measures/MeasureHolder.h>
65 :
66 : #include <casacore/scimath/Mathematics/MathFunc.h>
67 :
68 : using namespace casacore;
69 : namespace casa { //# NAMESPACE CASA - BEGIN
70 :
71 1496 : PBMath::PBMath()
72 : {
73 1496 : pb_pointer_p = 0;
74 1496 : };
75 :
76 23 : PBMath::PBMath(const RecordInterface& rec){
77 23 : pb_pointer_p=pbMathInterfaceFromRecord(rec);
78 :
79 :
80 23 : }
81 26 : PBMathInterface* PBMath::pbMathInterfaceFromRecord(const RecordInterface& rec){
82 26 : String name;
83 26 : PBMathInterface *pbPoint=nullptr;
84 26 : const Int nameFieldNumber=rec.fieldNumber("name");
85 26 : if (nameFieldNumber!=-1)
86 26 : name = rec.asString(nameFieldNumber);
87 :
88 : //cerr << "Name " << name << endl;
89 :
90 26 : if (name == "COMMONPB") {
91 13 : String commonpb;
92 : PBMath::CommonPB commonpbEnum;
93 13 : commonpb = rec.asString (rec.fieldNumber("commonpb"));
94 13 : PBMath::enumerateCommonPB(commonpb, commonpbEnum);
95 13 : pbPoint=pbMathInterfaceForCommonPB(commonpbEnum, False);
96 26 : } else if (name == "AIRY") {
97 :
98 :
99 0 : Quantity dishdiam;
100 0 : getQuantity(rec, "dishdiam", dishdiam);
101 0 : Quantity blockagediam;
102 0 : getQuantity(rec, "blockagediam", blockagediam);
103 0 : Quantity reffreq;
104 0 : getQuantity(rec, "reffreq", reffreq);
105 0 : Quantity maxrad;
106 0 : getQuantity(rec, "maxrad", maxrad);
107 0 : MDirection squintdir;
108 0 : getMDirection(rec, "squintdir", squintdir);
109 0 : Quantity squintreffreq;
110 0 : getQuantity(rec, "squintreffreq", squintreffreq);
111 : Bool usesymmetricbeam;
112 0 : usesymmetricbeam=rec.asBool( rec.fieldNumber("usesymmetricbeam"));
113 :
114 0 : pbPoint = new PBMath1DAiry( dishdiam, blockagediam,
115 : maxrad, reffreq,
116 0 : BeamSquint(squintdir, squintreffreq),
117 0 : usesymmetricbeam);
118 :
119 13 : } else if (name == "GAUSS") {
120 :
121 :
122 0 : Quantity halfwidth;
123 0 : getQuantity(rec, "halfwidth", halfwidth);
124 0 : Quantity reffreq;
125 0 : getQuantity(rec, "reffreq", reffreq);
126 0 : Quantity maxrad;
127 0 : getQuantity(rec, "maxrad", maxrad);
128 0 : MDirection squintdir;
129 0 : getMDirection(rec, "squintdir", squintdir);
130 0 : Quantity squintreffreq;
131 0 : getQuantity(rec, "squintreffreq", squintreffreq);
132 : Bool usesymmetricbeam;
133 0 : usesymmetricbeam=rec.asBool( rec.fieldNumber("usesymmetricbeam"));
134 : Bool isthisvp;
135 0 : isthisvp=rec.asBool( rec.fieldNumber("isthisvp"));
136 :
137 0 : pbPoint = new PBMath1DGauss( halfwidth, maxrad, reffreq,
138 : isthisvp,
139 0 : BeamSquint(squintdir, squintreffreq),
140 0 : usesymmetricbeam);
141 :
142 13 : } else if (name == "POLY") {
143 :
144 :
145 13 : Vector<Double> coeff;
146 13 : coeff=rec.asArrayDouble( rec.fieldNumber("coeff"));
147 13 : Quantity reffreq;
148 13 : getQuantity(rec, "reffreq", reffreq);
149 13 : Quantity maxrad;
150 13 : getQuantity(rec, "maxrad", maxrad);
151 13 : MDirection squintdir;
152 13 : getMDirection(rec, "squintdir", squintdir);
153 13 : Quantity squintreffreq;
154 13 : getQuantity(rec, "squintreffreq", squintreffreq);
155 : Bool usesymmetricbeam;
156 13 : usesymmetricbeam=rec.asBool( rec.fieldNumber("usesymmetricbeam"));
157 : Bool isthisvp;
158 13 : isthisvp=rec.asBool( rec.fieldNumber("isthisvp"));
159 :
160 26 : pbPoint = new PBMath1DPoly( coeff, maxrad, reffreq,
161 : isthisvp,
162 26 : BeamSquint(squintdir, squintreffreq),
163 13 : usesymmetricbeam);
164 :
165 :
166 13 : } else if (name == "IPOLY") {
167 :
168 0 : Quantity reffreq;
169 0 : getQuantity(rec, "reffreq", reffreq);
170 0 : Quantity maxrad;
171 0 : getQuantity(rec, "maxrad", maxrad);
172 0 : MDirection squintdir;
173 0 : getMDirection(rec, "squintdir", squintdir);
174 0 : Quantity squintreffreq;
175 0 : getQuantity(rec, "squintreffreq", squintreffreq);
176 : Bool usesymmetricbeam;
177 0 : usesymmetricbeam=rec.asBool( rec.fieldNumber("usesymmetricbeam"));
178 : Bool isthisvp;
179 0 : isthisvp=rec.asBool( rec.fieldNumber("isthisvp"));
180 0 : Bool wide = rec.isDefined("fitfreqs");
181 0 : if (wide) {
182 0 : Matrix <Double> coeff;
183 0 : coeff=rec.asArrayDouble( rec.fieldNumber("coeff"));
184 0 : Vector<Double> freqs;
185 0 : freqs=rec.asArrayDouble( rec.fieldNumber("fitfreqs"));
186 0 : pbPoint = new PBMath1DIPoly( coeff, freqs, maxrad, reffreq,
187 : isthisvp,
188 0 : BeamSquint(squintdir, squintreffreq),
189 0 : usesymmetricbeam);
190 0 : } else {
191 0 : Vector<Double> coeff;
192 0 : coeff=rec.asArrayDouble( rec.fieldNumber("coeff"));
193 0 : pbPoint = new PBMath1DIPoly( coeff, maxrad, reffreq,
194 : isthisvp,
195 0 : BeamSquint(squintdir, squintreffreq),
196 0 : usesymmetricbeam);
197 :
198 0 : }
199 :
200 0 : } else if (name == "COSPOLY") {
201 :
202 :
203 0 : Vector<Double> coeff;
204 0 : coeff=rec.asArrayDouble( rec.fieldNumber("coeff"));
205 0 : Vector<Double> scale;
206 0 : scale=rec.asArrayDouble( rec.fieldNumber("scale"));
207 0 : Quantity reffreq;
208 0 : getQuantity(rec, "reffreq", reffreq);
209 0 : Quantity maxrad;
210 0 : getQuantity(rec, "maxrad", maxrad);
211 0 : MDirection squintdir;
212 0 : getMDirection(rec, "squintdir", squintdir);
213 0 : Quantity squintreffreq;
214 0 : getQuantity(rec, "squintreffreq", squintreffreq);
215 : Bool usesymmetricbeam;
216 0 : usesymmetricbeam=rec.asBool( rec.fieldNumber("usesymmetricbeam"));
217 : Bool isthisvp;
218 0 : isthisvp=rec.asBool( rec.fieldNumber("isthisvp"));
219 :
220 0 : pbPoint = new PBMath1DCosPoly( coeff, scale, maxrad, reffreq,
221 : isthisvp,
222 0 : BeamSquint(squintdir, squintreffreq),
223 0 : usesymmetricbeam);
224 :
225 0 : } else if (name == "NUMERIC") {
226 0 : Vector<Double> vect;
227 0 : vect=rec.asArrayDouble( rec.fieldNumber("vect"));
228 0 : Quantity reffreq;
229 0 : getQuantity(rec, "reffreq", reffreq);
230 0 : Quantity maxrad;
231 0 : getQuantity(rec, "maxrad", maxrad);
232 0 : MDirection squintdir;
233 0 : getMDirection(rec, "squintdir", squintdir);
234 0 : Quantity squintreffreq;
235 0 : getQuantity(rec, "squintreffreq", squintreffreq);
236 : Bool usesymmetricbeam;
237 0 : usesymmetricbeam=rec.asBool( rec.fieldNumber("usesymmetricbeam"));
238 : Bool isthisvp;
239 0 : isthisvp=rec.asBool( rec.fieldNumber("isthisvp"));
240 :
241 0 : Vector<Float> aTempVector;
242 0 : aTempVector.resize(vect.nelements());
243 0 : for (uInt i = 0; i< vect.nelements(); i++) {
244 0 : aTempVector(i) = (Float)vect(i);
245 : }
246 : //cerr << "reffreq "<< reffreq << " maxrad " << maxrad << " isvp " << isthisvp << endl;
247 :
248 0 : pbPoint = new PBMath1DNumeric( aTempVector, maxrad, reffreq,
249 : isthisvp,
250 0 : BeamSquint(squintdir, squintreffreq),
251 0 : usesymmetricbeam);
252 0 : } else if (name == "IMAGE") {
253 0 : if(rec.isDefined("compleximage")){
254 0 : String complexIm;
255 0 : rec.get("compleximage", complexIm);
256 0 : if(Table::isReadable(complexIm)) {
257 0 : PagedImage<Complex> cim(complexIm);
258 0 : pbPoint=new PBMath2DImage(cim);
259 0 : }
260 : else{
261 0 : throw(AipsError("Complex Image "+ complexIm + " is not readable"));
262 : }
263 0 : }
264 : else{
265 0 : String realImageName, imagImageName;
266 0 : realImageName=rec.asString(rec.fieldNumber("realimage"));
267 0 : imagImageName=rec.asString(rec.fieldNumber("imagimage"));
268 :
269 0 : PagedImage<Float> realImage(realImageName);
270 0 : if(Table::isReadable(imagImageName)) {
271 0 : PagedImage<Float> imagImage(imagImageName);
272 0 : pbPoint = new PBMath2DImage(realImage, imagImage);
273 0 : }
274 : else {
275 0 : if(!Table::isReadable(imagImageName)) {
276 0 : pbPoint = new PBMath2DImage(realImage);
277 : }
278 : else {
279 0 : throw(AipsError("Image "+ realImageName + " is not readable or does not exist"));
280 : }
281 : }
282 0 : }
283 : }
284 : else{
285 0 : throw AipsError("Unrecognized PB name: " + name);
286 : }
287 26 : return pbPoint;
288 26 : };
289 :
290 0 : PBMath::PBMath(PBMath::CommonPB myPBType, Bool useSymmetricBeam)
291 : {
292 :
293 :
294 0 : initByTelescope(myPBType, useSymmetricBeam);
295 :
296 0 : };
297 :
298 843 : PBMath::PBMath(String& telescopeName, Bool useSymmetricBeam, Quantity freq){
299 :
300 843 : String band;
301 : PBMath::CommonPB whichPB;
302 843 : String pbName;
303 843 : whichCommonPBtoUse(telescopeName, freq, band, whichPB, pbName);
304 843 : initByTelescope(whichPB, useSymmetricBeam, freq.getValue("Hz"));
305 :
306 :
307 843 : }
308 :
309 0 : PBMath::PBMath(Double dishDiam, Bool useSymmetricBeam, Quantity freq){
310 :
311 0 : initByDiameter(dishDiam, useSymmetricBeam, freq.getValue("Hz"));
312 :
313 0 : }
314 :
315 : // Explicitly call each letter class's constructor
316 : // PBClass is for cases where we cannot distinquish the
317 : // PB's to be made from the constructor arguments
318 0 : PBMath::PBMath(PBMathInterface::PBClass theClass, Quantity halfWidth,
319 : Quantity maxRad, Quantity refFreq,
320 : Bool isThisVP,
321 0 : BeamSquint squint, Bool useSymmetricBeam )
322 : {
323 0 : LogIO os(LogOrigin("PBMath", "PBMath"));
324 0 : if (theClass == PBMathInterface::GAUSS) {
325 : pb_pointer_p = new PBMath1DGauss(halfWidth, maxRad, refFreq,
326 0 : isThisVP, squint, useSymmetricBeam);
327 : } else {
328 0 : os << "Unrecognized voltage pattern class type for this constructor:1" << LogIO::EXCEPTION;
329 : }
330 0 : };
331 :
332 :
333 0 : PBMath::PBMath(PBMathInterface::PBClass theClass, const Vector<Double>& avector, Quantity maxRad,
334 0 : Quantity refFreq, Bool isThisVP, BeamSquint squint, Bool useSymmetricBeam )
335 : {
336 0 : LogIO os(LogOrigin("PBMath", "PBMath"));
337 0 : if (theClass == PBMathInterface::POLY) {
338 : pb_pointer_p = new PBMath1DPoly(avector, maxRad, refFreq,
339 0 : isThisVP, squint, useSymmetricBeam);
340 0 : } else if (theClass == PBMathInterface::IPOLY) {
341 : pb_pointer_p = new PBMath1DIPoly(avector, maxRad, refFreq,
342 0 : isThisVP, squint, useSymmetricBeam);
343 0 : } else if (theClass == PBMathInterface::NUMERIC) {
344 0 : Vector<Float> aTempVector;
345 0 : aTempVector.resize(avector.nelements());
346 0 : for (uInt i = 0; i< avector.nelements(); i++) {
347 0 : aTempVector(i) = (Float)avector(i);
348 : }
349 : pb_pointer_p = new PBMath1DNumeric(aTempVector, maxRad, refFreq,
350 0 : isThisVP, squint, useSymmetricBeam);
351 0 : } else {
352 0 : os << "Unrecognized voltage pattern class type for this constructor:2" << LogIO::EXCEPTION;
353 : }
354 0 : };
355 :
356 0 : PBMath::PBMath(PBMathInterface::PBClass theClass, const Vector<Float>& avector, Quantity maxRad,
357 0 : Quantity refFreq, Bool isThisVP, BeamSquint squint, Bool useSymmetricBeam )
358 : {
359 0 : LogIO os(LogOrigin("PBMath", "PBMath"));
360 0 : if (theClass == PBMathInterface::NUMERIC) {
361 : pb_pointer_p = new PBMath1DNumeric(avector, maxRad, refFreq,
362 0 : isThisVP, squint, useSymmetricBeam);
363 : } else {
364 0 : os << "Unrecognized voltage pattern class type for this constructor:2" << LogIO::EXCEPTION;
365 : }
366 0 : };
367 :
368 0 : PBMath::PBMath(PBMathInterface::PBClass theClass,
369 0 : ImageInterface<Float>& reJones)
370 : {
371 0 : LogIO os(LogOrigin("PBMath", "PBMath"));
372 0 : if (theClass == PBMathInterface::IMAGE) {
373 0 : pb_pointer_p = new PBMath2DImage(reJones);
374 : } else {
375 0 : os << "Unrecognized voltage pattern class type for this constructor:2" << LogIO::EXCEPTION;
376 : }
377 0 : };
378 :
379 0 : PBMath::PBMath(PBMathInterface::PBClass theClass,
380 : ImageInterface<Float>& reJones,
381 0 : ImageInterface<Float>& imJones)
382 : {
383 0 : LogIO os(LogOrigin("PBMath", "PBMath"));
384 0 : if (theClass == PBMathInterface::IMAGE) {
385 0 : pb_pointer_p = new PBMath2DImage(reJones, imJones);
386 : } else {
387 0 : os << "Unrecognized voltage pattern class type for this constructor:2" << LogIO::EXCEPTION;
388 : }
389 0 : };
390 :
391 :
392 0 : PBMath::PBMath(PBMathInterface::PBClass theClass, const Vector<Double>& coeff,
393 : const Vector<Double>& scales, Quantity maxRad, Quantity refFreq,
394 : Bool isThisVP,
395 : BeamSquint squint,
396 0 : Bool useSymmetricBeam )
397 : {
398 0 : LogIO os(LogOrigin("PBMath", "PBMath"));
399 0 : if (theClass == PBMathInterface::COSPOLY) {
400 : pb_pointer_p = new PBMath1DCosPoly(coeff, scales, maxRad, refFreq,
401 0 : isThisVP, squint, useSymmetricBeam);
402 : } else {
403 0 : os << "Unrecognized voltage pattern class type for this constructor:3" << LogIO::EXCEPTION;
404 : }
405 0 : };
406 :
407 0 : PBMath::PBMath(PBMathInterface::PBClass theClass, Quantity dishDiam,
408 : Quantity blockDiam, Quantity maxRad,
409 : Quantity refFreq, BeamSquint squint,
410 0 : Bool useSymmetricBeam )
411 : {
412 0 : LogIO os(LogOrigin("PBMath", "PBMath"));
413 0 : if (theClass == PBMathInterface::AIRY) {
414 : pb_pointer_p = new PBMath1DAiry(dishDiam, blockDiam, maxRad,
415 0 : refFreq, squint, useSymmetricBeam);
416 : } else {
417 0 : os << "Unrecognized voltage pattern class type for this constructor:4" << LogIO::EXCEPTION;
418 : }
419 0 : };
420 :
421 0 : PBMath::PBMath(const PBMath &other)
422 : {
423 0 : pb_pointer_p = other.pb_pointer_p;
424 0 : };
425 :
426 :
427 2340 : PBMath::~PBMath()
428 : {
429 2340 : };
430 :
431 1496 : PBMath& PBMath::operator=(const PBMath& other)
432 : {
433 1496 : if (this != &other) {
434 1496 : pb_pointer_p = other.pb_pointer_p;
435 : }
436 1496 : return *this;
437 : };
438 :
439 0 : Bool PBMath::operator==(const PBMath& other) const
440 : {
441 0 : return (pb_pointer_p == other.pb_pointer_p ? true : false);
442 : };
443 :
444 :
445 0 : Bool PBMath::operator!=(const PBMath& other) const
446 : {
447 0 : return (pb_pointer_p != other.pb_pointer_p ? true : false);
448 : };
449 :
450 :
451 : ImageRegion*
452 0 : PBMath::extent(const ImageInterface<Complex> & im, const MDirection& pointing,
453 : const Int row, const Float fPad, const Int iChan,
454 : const SkyJones::SizeType sizeType)
455 : {
456 0 : return (pb_pointer_p->extent(im, pointing, row, fPad, iChan, sizeType));
457 : };
458 :
459 : ImageRegion*
460 0 : PBMath::extent(const ImageInterface<Float> & im, const MDirection& pointing,
461 : const Int row, const Float fPad, const Int iChan,
462 : const SkyJones::SizeType sizeType)
463 : {
464 0 : return (pb_pointer_p->extent(im, pointing, row, fPad, iChan, sizeType));
465 : };
466 :
467 : Int
468 68 : PBMath::support(const CoordinateSystem& cs)
469 : {
470 68 : return (pb_pointer_p->support(cs));
471 : };
472 :
473 :
474 : //Bool PBMath::flushToTable(Table& beamSubTable, Int iRow)
475 : //{
476 : // pb_pointer_p->flushToTable(Table& beamSubTable, Int iRow);
477 : //};
478 :
479 : //PBMath PBMath::copy() const
480 : //{
481 : // PBMath tmp (*this);
482 : // tmp.pb_pointer_p = pb_pointer_p->clone();
483 : // return tmp;
484 : //};
485 :
486 :
487 : void
488 17 : PBMath::summary(Int nValues)
489 : {
490 17 : pb_pointer_p->summary(nValues);
491 17 : };
492 :
493 :
494 :
495 : ImageInterface<Complex>&
496 0 : PBMath::applyVP(const ImageInterface<Complex>& in,
497 : ImageInterface<Complex>& out,
498 : const MDirection& sp,
499 : const Quantity parAngle,
500 : const BeamSquint::SquintType doSquint,
501 : Bool inverse,
502 : Bool conjugate,
503 : Float cutoff,
504 : Bool forward)
505 : {
506 0 : return pb_pointer_p->applyVP(in, out, sp, parAngle, doSquint, inverse,
507 0 : conjugate, cutoff, forward);
508 : };
509 :
510 : ImageInterface<Complex>&
511 814 : PBMath::applyPB(const ImageInterface<Complex>& in,
512 : ImageInterface<Complex>& out,
513 : const MDirection& sp,
514 : const Quantity parAngle,
515 : const BeamSquint::SquintType doSquint,
516 : Bool inverse,
517 : Float cutoff,
518 : Bool forward)
519 : {
520 814 : return pb_pointer_p->applyPB(in, out, sp, parAngle, doSquint, inverse, cutoff, forward);
521 : };
522 :
523 :
524 : ImageInterface<Float>&
525 0 : PBMath::applyPB(const ImageInterface<Float>& in,
526 : ImageInterface<Float>& out,
527 : const MDirection& sp,
528 : const Quantity parAngle,
529 : const BeamSquint::SquintType doSquint,
530 : Float cutoff)
531 : {
532 0 : return pb_pointer_p->applyPB(in, out, sp, parAngle, doSquint, cutoff);
533 : };
534 :
535 : ImageInterface<Float>&
536 40 : PBMath::applyPB2(const ImageInterface<Float>& in,
537 : ImageInterface<Float>& out,
538 : const MDirection& sp,
539 : const Quantity parAngle,
540 : const BeamSquint::SquintType doSquint,
541 : Float cutoff)
542 : {
543 40 : return pb_pointer_p->applyPB2(in, out, sp, parAngle, doSquint, cutoff);
544 : };
545 : SkyComponent&
546 0 : PBMath::applyVP(SkyComponent& in,
547 : SkyComponent& out,
548 : const MDirection& sp,
549 : const Quantity frequency,
550 : const Quantity parAngle,
551 : const BeamSquint::SquintType doSquint,
552 : Bool inverse,
553 : Bool conjugate,
554 : Float cutoff,
555 : Bool forward)
556 : {
557 0 : return pb_pointer_p->applyVP(in, out, sp, frequency, parAngle, doSquint, inverse, conjugate,
558 0 : cutoff, forward);
559 : };
560 :
561 : SkyComponent&
562 2491 : PBMath::applyPB(SkyComponent& in,
563 : SkyComponent& out,
564 : const MDirection& sp,
565 : const Quantity frequency,
566 : const Quantity parAngle,
567 : const BeamSquint::SquintType doSquint,
568 : Bool inverse,
569 : Float cutoff,
570 : Bool forward)
571 : {
572 4982 : return pb_pointer_p->applyPB(in, out, sp, frequency, parAngle, doSquint, inverse, cutoff,
573 4982 : forward);
574 : };
575 :
576 : SkyComponent&
577 0 : PBMath::applyPB2(SkyComponent& in,
578 : SkyComponent& out,
579 : const MDirection& sp,
580 : const Quantity frequency,
581 : const Quantity parAngle,
582 : const BeamSquint::SquintType doSquint)
583 : {
584 0 : return pb_pointer_p->applyPB2(in, out, sp, frequency, parAngle, doSquint);
585 : };
586 :
587 :
588 :
589 : ImageInterface<Complex>&
590 0 : PBMath::apply(const ImageInterface<Complex>& in,
591 : ImageInterface<Complex>& out,
592 : const MDirection& sp,
593 : const Quantity parAngle,
594 : const BeamSquint::SquintType doSquint,
595 : Bool inverse,
596 : Bool conjugate,
597 : Int ipower,
598 : Float cutoff,
599 : Bool forward)
600 : {
601 0 : return pb_pointer_p->apply(in, out, sp, parAngle, doSquint, inverse,
602 0 : conjugate, ipower, cutoff, forward);
603 : };
604 :
605 :
606 : SkyComponent&
607 0 : PBMath::apply(SkyComponent& in,
608 : SkyComponent& out,
609 : const MDirection& sp,
610 : const Quantity frequency,
611 : const Quantity parAngle,
612 : const BeamSquint::SquintType doSquint,
613 : Bool inverse,
614 : Bool conjugate,
615 : Int ipower, // ie, 1=VP, 2=PB, 4=PB^2
616 : Float cutoff,
617 : Bool forward)
618 : {
619 0 : return pb_pointer_p->apply(in, out, sp, frequency, parAngle, doSquint, inverse,
620 0 : conjugate, ipower, cutoff, forward);
621 : };
622 :
623 :
624 :
625 636 : Bool PBMath::ok() const
626 : {
627 636 : if (pb_pointer_p.null()) {
628 0 : return false;
629 : } else {
630 636 : return pb_pointer_p->ok();
631 : }
632 : };
633 :
634 :
635 : void
636 1198 : PBMath::whichCommonPBtoUse(String &telescope, Quantity &freq,
637 : String &band, PBMath::CommonPB &whichPB, String &PBName)
638 : {
639 2396 : LogIO os(LogOrigin("PBMath", "PBMath"));
640 :
641 : // note: these bands are fairly fast and loose,
642 : // and owe a lot to the fact that the band coverage is sparse!
643 1198 : Double freqGHz = freq.getValue("GHz");
644 1198 : if (telescope(0,3)=="VLA") {
645 552 : if (freqGHz > 35.0 && freqGHz < 55.0) {
646 8 : whichPB = PBMath::VLA_Q;
647 8 : band = "Q";
648 544 : } else if (freqGHz > 35.0 && freqGHz < 55.0) {
649 0 : whichPB = PBMath::VLA_Q;
650 0 : band = "Q";
651 544 : } else if (freqGHz > 19.0 && freqGHz < 35.0) {
652 0 : whichPB = PBMath::VLA_K;
653 0 : band = "K";
654 544 : } else if (freqGHz > 11.0 && freqGHz < 19.0) {
655 0 : whichPB = PBMath::VLA_U;
656 0 : band = "U";
657 544 : } else if (freqGHz > 7.0 && freqGHz < 11.0) {
658 0 : whichPB = PBMath::VLA_X;
659 0 : band = "X";
660 544 : } else if (freqGHz > 4.0 && freqGHz < 7.0) {
661 0 : whichPB = PBMath::VLA_C;
662 0 : band = "C";
663 544 : } else if (freqGHz > 1.0 && freqGHz < 2.0) {
664 294 : whichPB = PBMath::VLA_L;
665 294 : band = "L";
666 250 : } else if (freqGHz > 0.2 && freqGHz < 0.4) {
667 0 : whichPB = PBMath::VLA_P;
668 0 : band = "P";
669 250 : } else if (freqGHz < 0.1) {
670 0 : whichPB = PBMath::VLA_4;
671 0 : band = "4";
672 : } else {
673 250 : whichPB = PBMath::VLA_NVSS;
674 250 : band = "UNKNOWN";
675 : }
676 646 : } else if(telescope(0,4)=="EVLA") {
677 412 : whichPB = PBMath::EVLA;
678 412 : band = "UNKNOWN";
679 234 : } else if(telescope(0,5)=="NGVLA") {
680 0 : whichPB = PBMath::NGVLA;
681 0 : band = "UNKNOWN";
682 234 : } else if (telescope(0,4)=="WSRT") {
683 0 : if (freqGHz > 3.0 && freqGHz < 6.0) {
684 0 : whichPB = PBMath::WSRT;
685 0 : band = "C";
686 0 : } else if (freqGHz > 0.80 && freqGHz < 2.0) {
687 0 : whichPB = PBMath::WSRT;
688 0 : band = "L";
689 0 : } else if (freqGHz > 0.5 && freqGHz < 0.8) {
690 0 : whichPB = PBMath::WSRT_LOW;
691 0 : band = "49"; // Hey, I just don't know its name: 49cm
692 0 : } else if (freqGHz < 0.5 ) {
693 0 : whichPB = PBMath::WSRT_LOW;
694 0 : band = "P";
695 : } else {
696 0 : whichPB = PBMath::WSRT;
697 0 : band = "UNKNOWN";
698 : }
699 234 : } else if (telescope(0,4)=="ATCA") {
700 0 : if (freqGHz >= 7.0 && freqGHz < 12.5) {
701 0 : whichPB = PBMath::ATCA_X;
702 0 : band = "X";
703 0 : } else if (freqGHz > 4.0 && freqGHz < 7.0) {
704 0 : whichPB = PBMath::ATCA_C;
705 0 : band = "C";
706 0 : } else if (freqGHz > 1.0 && freqGHz < 3.5) {
707 0 : whichPB = PBMath::ATCA_16;
708 0 : band = "16"; // 16cm band - combines L+S
709 0 : } else if (freqGHz > 15.0 && freqGHz < 26.0) {
710 0 : whichPB = PBMath::ATCA_K;
711 0 : band = "K";
712 0 : } else if (freqGHz > 30.0 && freqGHz < 50.0) {
713 0 : whichPB = PBMath::ATCA_Q;
714 0 : band = "Q";
715 0 : } else if (freqGHz > 80 && freqGHz < 120.0) {
716 0 : band = "W";
717 0 : whichPB = PBMath::ATCA_W;
718 : } else {
719 0 : whichPB = PBMath::ATCA_L1;
720 0 : band = "UNKNOWN";
721 : }
722 0 : if (telescope(0,9)=="ATCA_C_RI") {
723 : // Remy Indebetouw measured PB at 5.5GHz 2011/10/20
724 0 : whichPB = PBMath::ATCA_C_RI;
725 : }
726 0 : os << "Selecting primary beam for ATCA band: "+band<<LogIO::NORMAL;
727 234 : } else if (telescope(0,8)=="HATCREEK") {
728 0 : whichPB = PBMath::HATCREEK;
729 0 : band = "UNKNOWN";
730 234 : } else if (telescope(0,4)=="BIMA") { // BIMA is a synonym for HATCREEK
731 0 : whichPB = PBMath::HATCREEK;
732 0 : band = "UNKNOWN";
733 234 : } else if (telescope(0,3)=="GBT") {
734 80 : whichPB = PBMath::GBT;
735 80 : band = "UNKNOWN";
736 154 : } else if (telescope(0,4)=="GMRT") {
737 0 : whichPB = PBMath::GMRT;
738 0 : band = "UNKNOWN";
739 154 : } else if (telescope(0,4)=="OVRO") {
740 0 : whichPB = PBMath::OVRO;
741 0 : band = "UNKNOWN";
742 154 : } else if (telescope(0,7)=="NRAO12M") {
743 0 : whichPB = PBMath::NRAO12M;
744 0 : band = "UNKNOWN";
745 154 : } else if (telescope(0,7)=="IRAMPDB" || (telescope.contains("IRAM") && telescope.contains("PDB"))) {
746 0 : whichPB = PBMath::IRAMPDB;
747 0 : band = "UNKNOWN";
748 154 : } else if (telescope(0,7)=="IRAM30M") {
749 0 : whichPB = PBMath::IRAM30M;
750 0 : band = "UNKNOWN";
751 154 : } else if (telescope(0,9)=="NRAO140FT") {
752 0 : whichPB = PBMath::NRAO140FT;
753 0 : band = "UNKNOWN";
754 154 : } else if (telescope(0,4)=="ALMA") {
755 152 : whichPB = PBMath::ALMA;
756 152 : band = "UNKNOWN";
757 2 : } else if (telescope(0,6)=="ALMASD") {
758 0 : whichPB = PBMath::ALMASD;
759 0 : band = "UNKNOWN";
760 2 : } else if (telescope(0,3)=="ACA") {
761 2 : whichPB = PBMath::ACA;
762 2 : band = "UNKNOWN";
763 0 : } else if (telescope(0,3)=="SMA"){
764 0 : whichPB = PBMath::SMA;
765 0 : band= "UNKNOWN";
766 0 : } else if (telescope(0,3)=="ATA") {
767 0 : whichPB = PBMath::ATA;
768 0 : band = "UNKNOWN";
769 0 : } else if (telescope(0,3)=="ATF") {
770 0 : whichPB = PBMath::ALMA;
771 0 : band = "UNKNOWN";
772 0 : } else if (telescope(0,4)=="NONE") {
773 0 : whichPB = PBMath::NONE;
774 0 : band = "UNKNOWN";
775 : } else {
776 0 : whichPB = PBMath::UNKNOWN;
777 0 : band = "UNKNOWN";
778 : }
779 :
780 1198 : nameCommonPB(whichPB, PBName);
781 1198 : };
782 :
783 :
784 : // converts the enumerated type into a string
785 10437 : void PBMath::nameCommonPB(const PBMath::CommonPB iPB, String & str)
786 : {
787 :
788 10437 : switch (iPB) {
789 0 : case PBMath::UNKNOWN:
790 0 : str = "UNKNOWN";
791 0 : break;
792 0 : case PBMath::DEFAULT:
793 0 : str = "DEFAULT";
794 0 : break;
795 224 : case PBMath::ATCA_L1:
796 224 : str = "ATCA_L1";
797 224 : break;
798 224 : case PBMath::ATCA_L2:
799 224 : str = "ATCA_L2";
800 224 : break;
801 224 : case PBMath::ATCA_L3:
802 224 : str = "ATCA_L3";
803 224 : break;
804 224 : case PBMath::ATCA_16:
805 224 : str = "ATCA_16";
806 224 : break;
807 224 : case PBMath::ATCA_S:
808 224 : str = "ATCA_S";
809 224 : break;
810 224 : case PBMath::ATCA_C:
811 224 : str = "ATCA_C";
812 224 : break;
813 224 : case PBMath::ATCA_C_RI:
814 224 : str = "ATCA_C_RI";
815 224 : break;
816 224 : case PBMath::ATCA:
817 224 : str = "ATCA";
818 224 : break;
819 224 : case PBMath::ATCA_X:
820 224 : str = "ATCA_X";
821 224 : break;
822 224 : case PBMath::ATCA_K:
823 224 : str = "ATCA_K";
824 224 : break;
825 224 : case PBMath::ATCA_Q:
826 224 : str = "ATCA_Q";
827 224 : break;
828 224 : case PBMath::ATCA_W:
829 224 : str = "ATCA_W";
830 224 : break;
831 322 : case PBMath::GBT:
832 322 : str = "GBT";
833 322 : break;
834 224 : case PBMath::GMRT:
835 224 : str = "GMRT";
836 224 : break;
837 224 : case PBMath::HATCREEK:
838 224 : str = "HATCREEK";
839 224 : break;
840 224 : case PBMath::NRAO12M:
841 224 : str = "NRAO12M";
842 224 : break;
843 224 : case PBMath::IRAMPDB:
844 224 : str = "IRAMPDB";
845 224 : break;
846 224 : case PBMath::IRAM30M:
847 224 : str = "IRAM30M";
848 224 : break;
849 224 : case PBMath::OVRO:
850 224 : str = "OVRO";
851 224 : break;
852 636 : case PBMath::EVLA:
853 636 : str = "EVLA";
854 636 : break;
855 224 : case PBMath::NGVLA:
856 224 : str = "NGVLA";
857 224 : break;
858 224 : case PBMath::VLA:
859 224 : str = "VLA";
860 224 : break;
861 224 : case PBMath::VLA_INVERSE:
862 224 : str = "VLA_INVERSE";
863 224 : break;
864 474 : case PBMath::VLA_NVSS:
865 474 : str = "VLA_NVSS";
866 474 : break;
867 224 : case PBMath::VLA_2NULL:
868 224 : str = "VLA_2NULL";
869 224 : break;
870 224 : case PBMath::VLA_4:
871 224 : str = "VLA_4";
872 224 : break;
873 224 : case PBMath::VLA_P:
874 224 : str = "VLA_P";
875 224 : break;
876 518 : case PBMath::VLA_L:
877 518 : str = "VLA_L";
878 518 : break;
879 224 : case PBMath::VLA_C:
880 224 : str = "VLA_C";
881 224 : break;
882 224 : case PBMath::VLA_X:
883 224 : str = "VLA_X";
884 224 : break;
885 224 : case PBMath::VLA_U:
886 224 : str = "VLA_U";
887 224 : break;
888 224 : case PBMath::VLA_K:
889 224 : str = "VLA_K";
890 224 : break;
891 235 : case PBMath::VLA_Q:
892 235 : str = "VLA_Q";
893 235 : break;
894 224 : case PBMath::WSRT:
895 224 : str = "WSRT";
896 224 : break;
897 224 : case PBMath::WSRT_LOW:
898 224 : str = "WSRT_LOW";
899 224 : break;
900 399 : case PBMath::ALMA:
901 399 : str = "ALMA";
902 399 : break;
903 234 : case PBMath::ALMASD:
904 234 : str = "ALMASD";
905 234 : break;
906 227 : case PBMath::ACA:
907 227 : str = "ACA";
908 227 : break;
909 224 : case PBMath::SMA:
910 224 : str = "SMA";
911 224 : break;
912 224 : case PBMath::ATA:
913 224 : str = "ATA";
914 224 : break;
915 0 : case PBMath::NONE:
916 0 : str = "NONE";
917 0 : break;
918 224 : default:
919 224 : str = "UNKNOWN";
920 224 : break;
921 : }
922 10437 : };
923 :
924 :
925 839 : void PBMath::enumerateCommonPB(const String & str, PBMath::CommonPB& ipb)
926 : {
927 839 : if (str == "UNKNOWN") {
928 0 : ipb = PBMath::UNKNOWN;
929 839 : } else if (str == "DEFAULT") {
930 0 : ipb = PBMath::DEFAULT;
931 839 : } else if (str == "ATCA_L1") {
932 0 : ipb = PBMath::ATCA_L1;
933 839 : } else if (str == "ATCA_L2") {
934 0 : ipb = PBMath::ATCA_L2;
935 839 : } else if (str == "ATCA_L3") {
936 0 : ipb = PBMath::ATCA_L3;
937 839 : } else if (str == "ATCA_16") {
938 0 : ipb = PBMath::ATCA_16;
939 839 : } else if (str == "ATCA_S") {
940 0 : ipb = PBMath::ATCA_S;
941 839 : } else if (str == "ATCA_C") {
942 0 : ipb = PBMath::ATCA_C;
943 839 : } else if (str == "ATCA_C_RI") {
944 0 : ipb = PBMath::ATCA_C_RI;
945 839 : } else if (str == "ATCA") {
946 0 : ipb = PBMath::ATCA;
947 839 : } else if (str == "ATCA_X") {
948 0 : ipb = PBMath::ATCA_X;
949 839 : } else if (str == "ATCA_K") {
950 0 : ipb = PBMath::ATCA_K;
951 839 : } else if (str == "ATCA_Q") {
952 0 : ipb = PBMath::ATCA_Q;
953 839 : } else if (str == "ATCA_W") {
954 0 : ipb = PBMath::ATCA_W;
955 839 : } else if (str == "HATCREEK") {
956 0 : ipb = PBMath::HATCREEK;
957 839 : } else if (str == "BIMA") { // BIMA is a synonym for HATCREEK
958 0 : ipb = PBMath::HATCREEK;
959 839 : } else if (str == "EVLA"){
960 202 : ipb = PBMath::EVLA;
961 637 : } else if (str == "NGVLA"){
962 0 : ipb = PBMath::NGVLA;
963 637 : }else if (str == "VLA" ) {
964 523 : ipb = PBMath::VLA;
965 114 : } else if (str == "VLA_INVERSE") {
966 0 : ipb = PBMath::VLA_INVERSE;
967 114 : } else if (str == "VLA_NVSS") {
968 0 : ipb = PBMath::VLA_NVSS;
969 114 : } else if (str == "VLA_2NULL") {
970 0 : ipb = PBMath::VLA_2NULL;
971 114 : } else if (str == "VLA_4") {
972 0 : ipb = PBMath::VLA_4;
973 114 : } else if (str == "VLA_P") {
974 0 : ipb = PBMath::VLA_P;
975 114 : } else if (str == "VLA_L") {
976 0 : ipb = PBMath::VLA_L;
977 114 : } else if (str == "VLA_C") {
978 0 : ipb = PBMath::VLA_C;
979 114 : } else if (str == "VLA_X") {
980 0 : ipb = PBMath::VLA_X;
981 114 : } else if (str == "VLA_U") {
982 0 : ipb = PBMath::VLA_U;
983 114 : } else if (str == "VLA_K") {
984 0 : ipb = PBMath::VLA_K;
985 114 : } else if (str == "VLA_Q") {
986 0 : ipb = PBMath::VLA_Q;
987 114 : } else if (str == "WSRT") {
988 0 : ipb = PBMath::WSRT;
989 114 : } else if (str == "WSRT_LOW") {
990 0 : ipb = PBMath::WSRT_LOW;
991 114 : } else if (str == "OVRO") {
992 0 : ipb = PBMath::OVRO;
993 114 : } else if (str == "GBT") {
994 0 : ipb = PBMath::GBT;
995 114 : } else if (str == "GMRT") {
996 0 : ipb = PBMath::GMRT;
997 114 : } else if (str == "NRAO12M") {
998 0 : ipb = PBMath::NRAO12M;
999 114 : } else if (str == "IRAMPDB") {
1000 0 : ipb = PBMath::IRAMPDB;
1001 114 : } else if (str == "IRAM30M") {
1002 0 : ipb = PBMath::IRAM30M;
1003 114 : } else if (str == "NRAO140FT") {
1004 0 : ipb = PBMath::NRAO140FT;
1005 114 : } else if (str == "ALMA") {
1006 77 : ipb = PBMath::ALMA;
1007 37 : } else if (str == "ALMASD") {
1008 10 : ipb = PBMath::ALMASD;
1009 27 : } else if (str == "ACA") {
1010 1 : ipb = PBMath::ACA;
1011 26 : } else if (str == "SMA"){
1012 0 : ipb = PBMath::SMA;
1013 26 : } else if (str == "ATA"){
1014 0 : ipb = PBMath::ATA;
1015 26 : } else if (str == "ATF") {
1016 0 : ipb = PBMath::ALMA;
1017 26 : } else if (str == "NONE") {
1018 13 : ipb = PBMath::NONE;
1019 : } else {
1020 13 : ipb = PBMath::UNKNOWN;
1021 : }
1022 839 : };
1023 :
1024 : Bool
1025 39 : PBMath::getQuantity(const RecordInterface& rec, const String& item,
1026 : Quantity& returnedQuantity)
1027 : {
1028 39 : String error;
1029 39 : QuantumHolder h;
1030 39 : const RecordInterface& quantRec = rec.asRecord(item);
1031 :
1032 39 : if (!h.fromRecord(error, quantRec)) {
1033 0 : throw (AipsError ("PBMath::getQuantity - could not recover "+item+" because "+error));
1034 : }
1035 39 : if (!h.isQuantity()) {
1036 0 : throw (AipsError ("PBMath::getQuantity - could not recover "+item+
1037 0 : " because it isnt a Quantity"));
1038 : }
1039 39 : returnedQuantity = h.asQuantumDouble();
1040 39 : return true;
1041 39 : };
1042 :
1043 : Bool
1044 13 : PBMath::getMDirection(const RecordInterface& rec, const String& item,
1045 : MDirection& returnedMDirection)
1046 : {
1047 13 : String error;
1048 13 : MeasureHolder h;
1049 13 : const RecordInterface& measureRec = rec.asRecord(item);
1050 :
1051 13 : if (!h.fromRecord(error, measureRec)) {
1052 0 : throw (AipsError ("PBMath::getMDirection - could not recover "+item+" because "+error));
1053 : }
1054 13 : if (!h.isMDirection()) {
1055 0 : throw (AipsError ("PBMath::getMDirection - could not recover "+item+
1056 0 : " because it isnt a MDirection"));
1057 : }
1058 13 : returnedMDirection = h.asMDirection();
1059 13 : return true;
1060 13 : };
1061 :
1062 :
1063 0 : void PBMath::initByDiameter(Double diameter, Bool /*useSymmetricBeam*/,
1064 : Double /*frequency*/){
1065 :
1066 0 : LogIO os(LogOrigin("PBMath", "initByDiameter", WHERE));
1067 0 : os << "PBMath init to Airy scaled to diameter = "<<diameter<<LogIO::POST;
1068 :
1069 : // This attempts to reproduce the AIRY pattern VLA PB
1070 0 : Vector<Float> vlanum(19);
1071 0 : vlanum(0) = 1.000000;
1072 0 : vlanum(1) = 0.997634;
1073 0 : vlanum(2) = 0.972516;
1074 0 : vlanum(3) = 0.913722;
1075 0 : vlanum(4) = 0.837871;
1076 0 : vlanum(5) = 0.750356;
1077 0 : vlanum(6) = 0.651549;
1078 0 : vlanum(7) = 0.549903;
1079 0 : vlanum(8) = 0.449083;
1080 0 : vlanum(9) = 0.352819;
1081 0 : vlanum(10) = 0.266025;
1082 0 : vlanum(11) = 0.190533;
1083 0 : vlanum(12) = 0.128047;
1084 0 : vlanum(13) = 0.0794855;
1085 0 : vlanum(14) = 0.0438381;
1086 0 : vlanum(15) = 0.0201386;
1087 0 : vlanum(16) = 0.0065117;
1088 0 : vlanum(17) = 0.000690041;
1089 0 : vlanum(18) = 8.87288e-05;
1090 :
1091 0 : if(diameter>0.){
1092 0 : Double scalesize = 1.1998662 * 25.0/diameter;
1093 0 : pb_pointer_p = new PBMath1DNumeric(vlanum, Quantity(scalesize,"'"),
1094 0 : Quantity(43.0,"GHz"), false);
1095 : }
1096 0 : else if(diameter==0.){
1097 0 : throw(AipsError("Attempt to initialize PBMath with zero dish diameter."));
1098 : }
1099 : else{
1100 0 : throw(AipsError("Attempt to initialize PBMath with negative dish diameter."));
1101 : }
1102 :
1103 0 : }
1104 :
1105 950 : PBMathInterface* PBMath::pbMathInterfaceForCommonPB(const PBMath::CommonPB ipb, Bool useSymmetricBeam){
1106 :
1107 950 : PBMathInterface* thepbMath=nullptr;
1108 1900 : LogIO os(LogOrigin("PBMath", "pbMathInterfaceForCommonPB"));
1109 950 : Vector<Double> vlacoef(4);
1110 950 : vlacoef(0)= 1.0;
1111 950 : vlacoef(1)= -1.300633e-03;
1112 950 : vlacoef(2)= 6.480550e-07;
1113 950 : vlacoef(3)= -1.267928e-10;
1114 :
1115 : // This attempts to reproduce the AIRY pattern VLA PB
1116 950 : Vector<Float> vlanum(19);
1117 950 : vlanum(0) = 1.000000;
1118 950 : vlanum(1) = 0.997634;
1119 950 : vlanum(2) = 0.972516;
1120 950 : vlanum(3) = 0.913722;
1121 950 : vlanum(4) = 0.837871;
1122 950 : vlanum(5) = 0.750356;
1123 950 : vlanum(6) = 0.651549;
1124 950 : vlanum(7) = 0.549903;
1125 950 : vlanum(8) = 0.449083;
1126 950 : vlanum(9) = 0.352819;
1127 950 : vlanum(10) = 0.266025;
1128 950 : vlanum(11) = 0.190533;
1129 950 : vlanum(12) = 0.128047;
1130 950 : vlanum(13) = 0.0794855;
1131 950 : vlanum(14) = 0.0438381;
1132 950 : vlanum(15) = 0.0201386;
1133 950 : vlanum(16) = 0.0065117;
1134 950 : vlanum(17) = 0.000690041;
1135 950 : vlanum(18) = 8.87288e-05;
1136 :
1137 :
1138 950 : switch (ipb) {
1139 286 : case EVLA:
1140 :
1141 : {
1142 :
1143 286 : thepbMath = new PBMath1DEVLA( Quantity(58.0,"'"), useSymmetricBeam, 1.0e9);
1144 : }
1145 286 : break;
1146 0 : case VLA:
1147 : // No squint assumed here.
1148 : // This beam has a peak of 1.00000 (unlike Napier and Rots; a refit to their data)
1149 : {
1150 0 : Vector<Double> coef(4);
1151 0 : coef(0)= 1.0;
1152 0 : coef(1)= -1.300633e-03;
1153 0 : coef(2)= 6.480550e-07;
1154 0 : coef(3)= -1.267928e-10;
1155 :
1156 0 : thepbMath = new PBMath1DPoly( coef, Quantity(58.0,"'"), Quantity(1.0,"GHz"));
1157 0 : }
1158 0 : break;
1159 0 : case VLA_INVERSE:
1160 : // No squint assumed here
1161 : // This beam does not have a peak of unity, kept for historical purposes
1162 : // This is the beam model used by CLASSIC AIPS
1163 : {
1164 0 : Vector<Double> coef(5);
1165 0 : coef(0)= 0.9920378;
1166 0 : coef(1)= 0.9956885e-3;
1167 0 : coef(2)= 0.3814573e-5;
1168 0 : coef(3)= -0.5311695e-8;
1169 0 : coef(4)= 0.3980963e-11;
1170 :
1171 0 : thepbMath = new PBMath1DIPoly( coef, Quantity(45.0,"'"), Quantity(1.0,"GHz"));
1172 0 : }
1173 0 : break;
1174 224 : case VLA_NVSS:
1175 : // No squint assumed here
1176 896 : thepbMath = new PBMath1DAiry( Quantity(24.5,"m"), Quantity(0.0,"m"),
1177 672 : Quantity(0.8564,"deg"), Quantity(1.0,"GHz"));
1178 224 : break;
1179 0 : case VLA_2NULL:
1180 : // No squint assumed here
1181 0 : thepbMath = new PBMath1DAiry( Quantity(24.5,"m"), Quantity(0.0,"m"),
1182 0 : Quantity(1.566,"deg"), Quantity(1.0,"GHz"));
1183 0 : break;
1184 :
1185 :
1186 294 : case VLA_L:
1187 : // Includes Squint.
1188 : // "mag" is the equivalent magnitude of the "half-squint", ie, difference
1189 : // between the nominal pointing center and the Stokes RR beam, but scaled
1190 : // to a reference frequency of 1 GHz.
1191 : // "ang" is calculated this way: look at the orientation of the feed
1192 : // looking DOWN on antenna (these are recorded in memos); flip to
1193 : // BEHIND the antenna, lookup up at sky; add 90 degrees to the angle for
1194 : // the squint orientation, which is perpendicular to the feed offset.
1195 : // The trick: the squint, measured CCW from +AZ, should be the same
1196 : // angle as the feed orientation measured CW from UP.
1197 : //
1198 : {
1199 294 : Float mag = 1.21; // half-squint magnitude in arcmin at 1 GHz)
1200 294 : Float ang = -135.0*C::pi/180.0; // squint orientation, rads, North of +AZ axis
1201 588 : thepbMath = new PBMath1DAiry( Quantity(25.0,"m"), Quantity(2.36,"m"),
1202 588 : Quantity(0.8564,"deg"), Quantity(1.0,"GHz"),
1203 588 : BeamSquint(MDirection(Quantity((mag*cos(ang)), "'"),
1204 588 : Quantity((mag*sin(ang)), "'"),
1205 588 : MDirection::Ref(MDirection::AZEL)),
1206 588 : Quantity(1.0, "GHz")),
1207 294 : useSymmetricBeam);
1208 : }
1209 294 : break;
1210 0 : case VLA_C:
1211 : {
1212 0 : Float mag = 1.21; // half-squint magnitude in arcmin at 1 GHz)
1213 0 : Float ang = 25.0*C::pi/180.0; // squint orientation, rads, North of +AZ axis
1214 0 : thepbMath = new PBMath1DAiry( Quantity(25.0,"m"), Quantity(2.36,"m"),
1215 0 : Quantity(0.8564,"deg"), Quantity(1.0,"GHz"),
1216 0 : BeamSquint(MDirection(Quantity((mag*cos(ang)), "'"),
1217 0 : Quantity((mag*sin(ang)), "'"),
1218 0 : MDirection::Ref(MDirection::AZEL)),
1219 0 : Quantity(1.0, "GHz")),
1220 0 : useSymmetricBeam);
1221 : }
1222 0 : break;
1223 0 : case VLA_X:
1224 : {
1225 : // This is based on the OLD X feed position; upgrade in progress
1226 0 : Float mag = 1.21; // half-squint magnitude in arcmin at 1 GHz)
1227 0 : Float ang = 82.0*C::pi/180.0; // squint orientation, rads, North of +AZ axis
1228 0 : thepbMath = new PBMath1DAiry( Quantity(25.0,"m"), Quantity(2.36,"m"),
1229 0 : Quantity(0.8564,"deg"), Quantity(1.0,"GHz"),
1230 0 : BeamSquint(MDirection(Quantity((mag*cos(ang)), "'"),
1231 0 : Quantity((mag*sin(ang)), "'"),
1232 0 : MDirection::Ref(MDirection::AZEL)),
1233 0 : Quantity(1.0, "GHz")),
1234 0 : useSymmetricBeam);
1235 : /*
1236 : MathFunc<Float> dd(SPHEROIDAL);
1237 : Vector<Float> valsph(31);
1238 : for(Int k=0; k <31; ++k){
1239 : valsph(k)=dd.value((Float)(k)/10.0);
1240 : }
1241 : //
1242 : Quantity fulrad((52.3/2.0*3.0/1.1117),"arcsec");
1243 : Quantity lefreq(144.0, "GHz");
1244 :
1245 : pb_pointer_p = new PBMath1DNumeric(valsph,fulrad,lefreq);
1246 : */
1247 : }
1248 0 : break;
1249 0 : case VLA_U:
1250 : {
1251 0 : Float mag = 1.21; // half-squint magnitude in arcmin at 1 GHz)
1252 0 : Float ang = -25.0*C::pi/180.0; // squint orientation, rads, North of +AZ axis
1253 0 : thepbMath = new PBMath1DAiry( Quantity(25.0,"m"), Quantity(2.36,"m"),
1254 0 : Quantity(0.8564,"deg"), Quantity(1.0,"GHz"),
1255 0 : BeamSquint(MDirection(Quantity((mag*cos(ang)), "'"),
1256 0 : Quantity((mag*sin(ang)), "'"),
1257 0 : MDirection::Ref(MDirection::AZEL)),
1258 0 : Quantity(1.0, "GHz")),
1259 0 : useSymmetricBeam);
1260 : }
1261 0 : break;
1262 0 : case VLA_K:
1263 : {
1264 : // This is based on the OLD K feed position; upgrade in progress
1265 0 : Float mag = 1.21; // half-squint magnitude in arcmin at 1 GHz)
1266 0 : Float ang = -6.0*C::pi/180.0; // squint orientation, rads, North of +AZ axis
1267 0 : thepbMath = new PBMath1DAiry( Quantity(25.0,"m"), Quantity(2.36,"m"),
1268 0 : Quantity(0.8564,"deg"), Quantity(1.0,"GHz"),
1269 0 : BeamSquint(MDirection(Quantity((mag*cos(ang)), "'"),
1270 0 : Quantity((mag*sin(ang)), "'"),
1271 0 : MDirection::Ref(MDirection::AZEL)),
1272 0 : Quantity(1.0, "GHz")),
1273 0 : useSymmetricBeam);
1274 : }
1275 0 : break;
1276 5 : case VLA_Q:
1277 : {
1278 : // This is based on an approximate feed position; awaiting new numbers
1279 5 : Float mag = 1.21; // half-squint magnitude in arcmin at 1 GHz)
1280 5 : Float ang = 0.0*C::pi/180.0; // squint orientation, rads, North of +AZ axis
1281 10 : thepbMath = new PBMath1DAiry( Quantity(25.0,"m"), Quantity(2.36,"m"),
1282 10 : Quantity(0.8564,"deg"), Quantity(1.0,"GHz"),
1283 10 : BeamSquint(MDirection(Quantity((mag*cos(ang)), "'"),
1284 10 : Quantity((mag*sin(ang)), "'"),
1285 10 : MDirection::Ref(MDirection::AZEL)),
1286 10 : Quantity(1.0, "GHz")),
1287 5 : useSymmetricBeam);
1288 : }
1289 5 : break;
1290 0 : case VLA_P:
1291 : {
1292 : // This is not based on any P band measurements; assume no squint
1293 0 : thepbMath = new PBMath1DAiry( Quantity(25.0,"m"), Quantity(2.36,"m"),
1294 0 : Quantity(0.8564,"deg"), Quantity(1.0,"GHz"));
1295 : }
1296 0 : break;
1297 0 : case VLA_4:
1298 : {
1299 : // This is not based on any 4 band measurements; assume no squint
1300 0 : thepbMath = new PBMath1DAiry( Quantity(25.0,"m"), Quantity(2.36,"m"),
1301 0 : Quantity(0.8564,"deg"), Quantity(1.0,"GHz"));
1302 : }
1303 0 : break;
1304 0 : case NGVLA:
1305 : // No squint assumed here or blockage
1306 0 : thepbMath = new PBMath1DAiry( Quantity(18.0,"m"), Quantity(0.0,"m"),
1307 0 : Quantity(1.5,"deg"), Quantity(1.0,"GHz"));
1308 : // cout << "Using NGVLA PBMath = 18m dish w/o blockage" << endl;
1309 0 : break;
1310 : // These 4 cases, OVRO, GMRT, NRAO12M, and NRAO140FT, do not represent
1311 : // well determined models (though some of them have well determined models).
1312 : // We are inserting them for completeness, though we should update some of them.
1313 : // For now we just assume that they are AIRY disks (via the Numeric fix, as with
1314 : // the VLA) of the approximate size (we fudge this by setting the scale-size
1315 : // to 1.1998662' * 25M/DIAM
1316 :
1317 0 : case OVRO:
1318 : {
1319 0 : Double scalesize = 1.1998662 * 25.0/10.4;
1320 0 : thepbMath = new PBMath1DNumeric(vlanum, Quantity(scalesize,"'"),
1321 0 : Quantity(43.0,"GHz"), false);
1322 : }
1323 0 : break;
1324 40 : case GBT:
1325 : {
1326 : //Measured beam was at 340MHz and spanning 1 deg.
1327 : //Data from Ron Maddalena
1328 : //If more data is available at other frequencies it should be cased
1329 : //out here otherwise this beam is going to be scaled for all frequencies
1330 40 : Double scalesize = 60.0*3.4e8/1.0e9;
1331 40 : Double factor=1.0e9/3.4e8;
1332 40 : Vector<Double> coeff(7);
1333 :
1334 40 : coeff[0]=1.00632057;
1335 40 : coeff[1]=-0.00181068914;
1336 40 : coeff[2]=1.32948262e-06;
1337 40 : coeff[3]=-4.72454392e-10;
1338 40 : coeff[4]=7.11800547e-14;
1339 40 : coeff[5]=-3.50646851e-20;
1340 40 : coeff[6]=-7.40533706e-22;
1341 :
1342 : //These coefficient are at 340 MHz
1343 : //For some peculiar reason PBMath1DPoly seem to scale things well only if
1344 : //the ref refreq at 1 GHz..
1345 : //hence converting coeeficients to the one at 1 GHz
1346 :
1347 40 : Double coeffac=1.0;
1348 280 : for (uInt k=1; k < 7 ; ++k){
1349 240 : coeffac=coeffac*factor*factor;
1350 240 : coeff[k]=coeff[k]*coeffac;
1351 : }
1352 :
1353 :
1354 160 : thepbMath = new PBMath1DPoly( coeff, Quantity(scalesize,"'"),
1355 120 : Quantity(1.0e9,"Hz"),false);
1356 40 : }
1357 40 : break;
1358 0 : case GMRT:
1359 : {
1360 0 : Double scalesize = 1.1998662 * 25.0/45.0;
1361 0 : thepbMath = new PBMath1DNumeric(vlanum, Quantity(scalesize,"'"),
1362 0 : Quantity(43.0,"GHz"), false);
1363 : }
1364 0 : break;
1365 0 : case NRAO12M:
1366 : {
1367 0 : Double scalesize = 1.1998662 * 25.0/12.0;
1368 0 : thepbMath = new PBMath1DNumeric(vlanum, Quantity(scalesize,"'"),
1369 0 : Quantity(43.0,"GHz"), false);
1370 : }
1371 0 : break;
1372 0 : case NRAO140FT:
1373 : {
1374 0 : Double scalesize = 1.1998662 * 25.0/43.0;
1375 0 : thepbMath = new PBMath1DNumeric(vlanum, Quantity(scalesize,"'"),
1376 0 : Quantity(43.0,"GHz"), false);
1377 : }
1378 0 : break;
1379 :
1380 :
1381 :
1382 0 : case WSRT:
1383 : {
1384 : // WSRT for freq > 800 MHz
1385 0 : Bool thisIsVP = true;
1386 0 : Vector<Double> coef(4);
1387 0 : Vector<Double> cosScale(4);
1388 0 : coef(0)= 0.0;
1389 0 : coef(1)= 0.0;
1390 0 : coef(2)= 0.0;
1391 0 : coef(3)= 1.0;
1392 0 : cosScale(0)= 0.0;
1393 0 : cosScale(1)= 0.0;
1394 0 : cosScale(2)= 0.0;
1395 0 : cosScale(3)= 0.01891;
1396 : // 0.01891 = 0.0672 * 1000(MHz/GHz) /(60(arcm/deg)) * 2pi/180
1397 :
1398 0 : thepbMath = new PBMath1DCosPoly(coef, cosScale,
1399 0 : Quantity(1.0, "deg"),
1400 0 : Quantity(1.0, "GHz"), thisIsVP);
1401 0 : }
1402 0 : break;
1403 0 : case WSRT_LOW:
1404 : {
1405 : // WSRT for freq < 800 MHz
1406 0 : Bool thisIsVP = true;
1407 0 : Vector<Double> coef(4);
1408 0 : Vector<Double> cosScale(4);
1409 0 : coef(0)= 0.0;
1410 0 : coef(1)= 0.0;
1411 0 : coef(2)= 0.0;
1412 0 : coef(3)= 1.0;
1413 0 : cosScale(0)= 0.0;
1414 0 : cosScale(1)= 0.0;
1415 0 : cosScale(2)= 0.0;
1416 0 : cosScale(3)= 0.01830;
1417 : // 0.01830 = 0.065 * 1000(MHz/GHz) /(60(arcm/deg)) * 2pi/180
1418 :
1419 0 : thepbMath = new PBMath1DCosPoly(coef, cosScale, Quantity(1.0, "deg"),
1420 0 : Quantity(1.0, "GHz"), thisIsVP);
1421 0 : }
1422 0 : break;
1423 0 : case ATCA_16:
1424 : {
1425 0 : Matrix<Double> coef(5,7);
1426 0 : coef(0,0)= 1.0;
1427 0 : coef(1,0)= 1.06274e-03;
1428 0 : coef(2,0)= 1.32342e-06;
1429 0 : coef(3,0)= -8.72013e-10;
1430 0 : coef(4,0)= 1.08020e-12;
1431 0 : coef(0,1)= 1.0;
1432 0 : coef(1,1)= 9.80817e-04;
1433 0 : coef(2,1)= 1.17898e-06;
1434 0 : coef(3,1)= -7.83160e-10;
1435 0 : coef(4,1)= 8.66199e-13;
1436 0 : coef(0,2)=1.0;
1437 0 : coef(1,2)= 9.53553e-04;
1438 0 : coef(2,2)= 9.33233e-07;
1439 0 : coef(3,2)= -4.26759e-10;
1440 0 : coef(4,2)= 5.63667e-13;
1441 0 : coef(0,3)= 1.0;
1442 0 : coef(1,3)= 9.78268e-04;
1443 0 : coef(2,3)= 6.63231e-07;
1444 0 : coef(3,3)= 4.18235e-11;
1445 0 : coef(4,3)= 2.62297e-13;
1446 0 : coef(0,4)= 1.0;
1447 0 : coef(1,4)= 1.02424e-03;
1448 0 : coef(2,4)= 6.12726e-07;
1449 0 : coef(3,4)= 2.25733e-10;
1450 0 : coef(4,4)= 2.04834e-13;
1451 0 : coef(0,5)= 1.0;
1452 0 : coef(1,5)= 1.05818e-03;
1453 0 : coef(2,5)= 5.37473e-07;
1454 0 : coef(3,5)= 4.22386e-10;
1455 0 : coef(4,5)= 1.17530e-13;
1456 0 : coef(0,6)= 1.0;
1457 0 : coef(1,6)= 1.10650e-03;
1458 0 : coef(2,6)= 5.11574e-07;
1459 0 : coef(3,6)= 5.89732e-10;
1460 0 : coef(4,6)= 8.13628e-14;
1461 0 : Vector<Double> freqs(7);
1462 0 : for (Int i=0; i<7; i++) {
1463 0 : freqs(i) = (1332+i*256)*1.e6;
1464 : }
1465 0 : thepbMath = new PBMath1DIPoly(coef, freqs, Quantity(53.0,"'"),
1466 0 : Quantity(1.0,"GHz"));
1467 0 : }
1468 0 : break;
1469 0 : case ATCA:
1470 : case ATCA_L1:
1471 : {
1472 0 : Vector<Double> coef(5);
1473 0 : coef(0)= 1.0;
1474 0 : coef(1)= 8.99e-4;
1475 0 : coef(2)= 2.15e-6;
1476 0 : coef(3)= -2.23e-9;
1477 0 : coef(4)= 1.56e-12;
1478 :
1479 0 : thepbMath = new PBMath1DIPoly( coef, Quantity(53.0,"'"), Quantity(1.0,"GHz"));
1480 0 : }
1481 0 : break;
1482 0 : case ATCA_L2:
1483 : {
1484 0 : Vector<Double> coef(7);
1485 0 : coef(0)= 1.0;
1486 0 : coef(1)= -1.0781341990755E-03;
1487 0 : coef(2)= 4.6179146405726E-07;
1488 0 : coef(3)= -1.0108079576125E-10;
1489 0 : coef(4)= 1.2073518438662E-14;
1490 0 : coef(5)= -7.5132629268134E-19;
1491 0 : coef(6)= 1.9083641820123E-23;
1492 :
1493 0 : thepbMath = new PBMath1DPoly( coef, Quantity(53.0,"'"), Quantity(1.0,"GHz"));
1494 0 : }
1495 0 : break;
1496 0 : case ATCA_L3:
1497 : {
1498 0 : os << "ATCA_L3 not yet implemented defaulting to L2 version" << LogIO::WARN << LogIO::POST;
1499 0 : Vector<Double> coef(7);
1500 0 : coef(0)= 1.0;
1501 0 : coef(1)= -1.0781341990755E-03;
1502 0 : coef(2)= 4.6179146405726E-07;
1503 0 : coef(3)= -1.0108079576125E-10;
1504 0 : coef(4)= 1.2073518438662E-14;
1505 0 : coef(5)= -7.5132629268134E-19;
1506 0 : coef(6)= 1.9083641820123E-23;
1507 0 : thepbMath = new PBMath1DPoly( coef, Quantity(53.0,"'"), Quantity(1.0,"GHz"));
1508 0 : }
1509 0 : break;
1510 0 : case ATCA_S:
1511 : {
1512 0 : Vector<Double> coef(5);
1513 0 : coef(0)= 1.0;
1514 0 : coef(1)= 1.02e-3;
1515 0 : coef(2)= 9.48e-7;
1516 0 : coef(3)= -3.68e-10;
1517 0 : coef(4)= 4.88e-13;
1518 :
1519 0 : thepbMath = new PBMath1DIPoly( coef, Quantity(53.0,"'"), Quantity(1.0,"GHz"),
1520 : false,
1521 0 : BeamSquint(MDirection(Quantity(0.0, "'"),
1522 0 : Quantity(0.0, "'"),
1523 0 : MDirection::Ref(MDirection::AZEL)),
1524 0 : Quantity(1.0, "GHz")),
1525 0 : false);
1526 0 : }
1527 0 : break;
1528 0 : case ATCA_C:
1529 : case ATCA_X:
1530 : {
1531 0 : Matrix<Double> coef(5,2);
1532 0 : coef(0,0)= 1.0;
1533 0 : coef(1,0)= 1.08e-3;
1534 0 : coef(2,0)= 1.31e-6;
1535 0 : coef(3,0)= -1.17e-9;
1536 0 : coef(4,0)= 1.07e-12;
1537 0 : coef(0,1)= 1.0;
1538 0 : coef(1,1)= 1.04e-3;
1539 0 : coef(2,1)= 8.36e-7;
1540 0 : coef(3,1)= -4.68e-10;
1541 0 : coef(4,1)= 5.50e-13;
1542 0 : Vector<Double> freqs(2);
1543 0 : freqs(0)=4.8e9;
1544 0 : freqs(1)=8.64e9;
1545 0 : thepbMath = new PBMath1DIPoly( coef, freqs, Quantity(53.0,"'"), Quantity(1.0,"GHz"),
1546 : false,
1547 0 : BeamSquint(MDirection(Quantity(0.0, "'"),
1548 0 : Quantity(0.0, "'"),
1549 0 : MDirection::Ref(MDirection::AZEL)),
1550 0 : Quantity(1.0, "GHz")),
1551 0 : false);
1552 0 : }
1553 0 : break;
1554 : {
1555 : Vector<Double> coef(5);
1556 :
1557 : thepbMath = new PBMath1DIPoly( coef, Quantity(53.0,"'"), Quantity(1.0,"GHz"),
1558 : false,
1559 : BeamSquint(MDirection(Quantity(0.0, "'"),
1560 : Quantity(0.0, "'"),
1561 : MDirection::Ref(MDirection::AZEL)),
1562 : Quantity(1.0, "GHz")),
1563 : false);
1564 : }
1565 : break;
1566 0 : case ATCA_C_RI:
1567 : // Remy Indebetouw measured the PB through the second sidelobe 20111020
1568 : {
1569 0 : os << "PBMath using 2011/10/22 5.5GHz PB" << LogIO::POST;
1570 0 : Vector<Float> coef(22);
1571 0 : coef( 0) = 1.00000;
1572 0 : coef( 1) = 0.98132;
1573 0 : coef( 2) = 0.96365;
1574 0 : coef( 3) = 0.87195;
1575 0 : coef( 4) = 0.75109;
1576 0 : coef( 5) = 0.62176;
1577 0 : coef( 6) = 0.48793;
1578 0 : coef( 7) = 0.34985;
1579 0 : coef( 8) = 0.21586;
1580 0 : coef( 9) = 0.10546;
1581 0 : coef(10) = 0.03669;
1582 0 : coef(11) = -0.03556;
1583 0 : coef(12) = -0.08266;
1584 0 : coef(13) = -0.12810;
1585 0 : coef(14) = -0.15440;
1586 0 : coef(15) = -0.16090;
1587 0 : coef(16) = -0.15360;
1588 0 : coef(17) = -0.13566;
1589 0 : coef(18) = -0.10666;
1590 0 : coef(19) = -0.06847;
1591 0 : coef(20) = -0.03136;
1592 0 : coef(21) = -0.00854;
1593 :
1594 0 : thepbMath = new PBMath1DNumeric( coef, Quantity(21.07,"'"),
1595 0 : Quantity(5.5,"GHz"),
1596 : true,
1597 0 : BeamSquint(MDirection(Quantity(0.0, "'"),
1598 0 : Quantity(0.0, "'"),
1599 0 : MDirection::Ref(MDirection::AZEL)),
1600 0 : Quantity(5.5, "GHz")),
1601 0 : true);
1602 0 : }
1603 0 : break;
1604 0 : case ATCA_K:
1605 : case ATCA_Q:
1606 : {
1607 0 : Vector<Double> coef(5);
1608 0 : coef(0)= 1.0;
1609 0 : coef(1)= 9.832e-4;
1610 0 : coef(2)= 1.081e-6;
1611 0 : coef(3)= -4.676e-10;
1612 0 : coef(4)= 6.650e-13;
1613 :
1614 0 : thepbMath = new PBMath1DIPoly( coef, Quantity(53.0,"'"), Quantity(1.0,"GHz"),
1615 : false,
1616 0 : BeamSquint(MDirection(Quantity(0.0, "'"),
1617 0 : Quantity(0.0, "'"),
1618 0 : MDirection::Ref(MDirection::AZEL)),
1619 0 : Quantity(1.0, "GHz")),
1620 0 : false);
1621 0 : }
1622 0 : break;
1623 0 : case ATCA_W:
1624 : {
1625 0 : Vector<Double> coef(5);
1626 0 : coef(0)= 1.0;
1627 0 : coef(1)= 1.271e-3;
1628 0 : coef(2)=-3.040e-7;
1629 0 : coef(3)= 1.410e-9;
1630 0 : coef(4)= 0;
1631 :
1632 0 : thepbMath = new PBMath1DIPoly( coef, Quantity(53.0,"'"), Quantity(1.0,"GHz"),
1633 : false,
1634 0 : BeamSquint(MDirection(Quantity(0.0, "'"),
1635 0 : Quantity(0.0, "'"),
1636 0 : MDirection::Ref(MDirection::AZEL)),
1637 0 : Quantity(1.0, "GHz")),
1638 0 : false);
1639 0 : }
1640 0 : break;
1641 0 : case HATCREEK:
1642 0 : thepbMath = new PBMath1DGauss( Quantity((191.67/2.0),"'"), // half width==> /2
1643 0 : Quantity(215.0, "'"),
1644 0 : Quantity(1.0, "GHz"),
1645 : false,
1646 0 : BeamSquint(MDirection(Quantity(0.0, "'"),
1647 0 : Quantity(0.0, "'"),
1648 0 : MDirection::Ref(MDirection::AZEL)),
1649 0 : Quantity(1.0, "GHz")),
1650 0 : false );
1651 0 : break;
1652 0 : case IRAMPDB:
1653 0 : thepbMath = new PBMath1DAiry( Quantity(15.0,"m"), Quantity(1.0,"m"),
1654 0 : Quantity(1.784,"deg"), Quantity(1.0,"GHz") );
1655 0 : break;
1656 0 : case IRAM30M:
1657 0 : thepbMath = new PBMath1DAiry( Quantity(15.0,"m"), Quantity(1.0,"m"),
1658 0 : Quantity(1.784,"deg"), Quantity(1.0,"GHz") );
1659 0 : break;
1660 87 : case ALMA:
1661 : /////Value suggested as more appropriate for all 3 ALMA dishes
1662 : ////previously was using the physical dimension of 12 m, 12m, and 6m
1663 : /////by Todd Hunter & Crystal Brogan 2011-12-06
1664 348 : thepbMath = new PBMath1DAiry( Quantity(10.7,"m"), Quantity(0.75,"m"),
1665 261 : Quantity(1.784,"deg"), Quantity(1.0,"GHz") );
1666 87 : break;
1667 0 : case ALMASD:
1668 0 : thepbMath = new PBMath1DAiry( Quantity(10.7,"m"), Quantity(0.75,"m"),
1669 0 : Quantity(1.784,"deg"), Quantity(1.0,"GHz") );
1670 0 : break;
1671 1 : case ACA:
1672 4 : thepbMath = new PBMath1DAiry( Quantity(6.25,"m"), Quantity(0.75,"m"),
1673 3 : Quantity(3.568,"deg"), Quantity(1.0,"GHz") );
1674 1 : break;
1675 0 : case SMA:
1676 : {
1677 : //Value provided by Crystal Brogan as of Dec-12-2007
1678 : //needs updating when proper values are given
1679 : /*
1680 : pb_pointer_p = new PBMath1DGauss( Quantity((52.3/2.0),"arcsec"), // half width==> /2
1681 : Quantity(35.0, "arcsec"),
1682 : Quantity(224.0, "GHz"),
1683 : false,
1684 : BeamSquint(MDirection(Quantity(0.0, "'"),
1685 : Quantity(0.0, "'"),
1686 : MDirection::Ref(MDirection::AZEL)),
1687 : Quantity(224.0, "GHz")),
1688 : false);
1689 :
1690 : //Crap beam no finite support in FFT
1691 :
1692 : */
1693 : //Using a spheroidal beam
1694 :
1695 0 : MathFunc<Float> dd(SPHEROIDAL);
1696 0 : Vector<Float> valsph(31);
1697 0 : for(Int k=0; k <31; ++k){
1698 0 : valsph(k)=dd.value((Float)(k)/10.0);
1699 : }
1700 : //
1701 0 : Quantity fulrad((52.3/2.0*3.0/1.1117),"arcsec");
1702 0 : Quantity lefreq(224.0, "GHz");
1703 :
1704 0 : thepbMath = new PBMath1DNumeric(valsph,fulrad,lefreq);
1705 0 : }
1706 0 : break;
1707 :
1708 0 : case ATA:
1709 0 : thepbMath = new PBMath1DAiry( Quantity(6.0,"m"), Quantity(0.5,"m"),
1710 0 : Quantity(3.568,"deg"), Quantity(1.0,"GHz") );
1711 0 : break;
1712 :
1713 13 : case NONE:
1714 : {
1715 13 : Vector<Double> coef(1);
1716 13 : coef(0)= 1.0;
1717 :
1718 26 : thepbMath = new PBMath1DIPoly( coef, Quantity(180,"deg"), Quantity(1.0,"GHz"),
1719 : false,
1720 26 : BeamSquint(MDirection(Quantity(0.0, "'"),
1721 26 : Quantity(0.0, "'"),
1722 26 : MDirection::Ref(MDirection::AZEL)),
1723 26 : Quantity(1.0, "GHz")),
1724 13 : false);
1725 13 : }
1726 13 : break;
1727 0 : default:
1728 0 : os << "Unrecognized CommonPB Type" << LogIO::EXCEPTION;
1729 0 : break;
1730 : }
1731 950 : return thepbMath;
1732 :
1733 950 : }
1734 :
1735 :
1736 843 : void PBMath::initByTelescope(PBMath::CommonPB myPBType,
1737 : Bool useSymmetricBeam,
1738 : Double /*frequency*/)
1739 : {
1740 :
1741 :
1742 843 : pb_pointer_p = pbMathInterfaceForCommonPB(myPBType, useSymmetricBeam);
1743 : // Remember, these are fit parameters for the PB, not the PB
1744 :
1745 :
1746 843 : };
1747 :
1748 : } //# NAMESPACE CASA - END
|
