Line data Source code
1 : //# PBMath1D.cc: Implementation for PBMath1D
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be 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 <casacore/casa/Arrays/Matrix.h>
32 : #include <casacore/casa/Arrays/Vector.h>
33 : #include <synthesis/TransformMachines/PBMath1D.h>
34 :
35 : #include <casacore/images/Regions/ImageRegion.h>
36 : #include <casacore/images/Images/ImageInterface.h>
37 :
38 : #include <components/ComponentModels/SkyComponent.h>
39 : #include <components/ComponentModels/Flux.h>
40 : #include <components/ComponentModels/ComponentShape.h>
41 :
42 : #include <casacore/lattices/Lattices/LatticeIterator.h>
43 : #include <casacore/lattices/Lattices/LatticeStepper.h>
44 : #include <casacore/lattices/LRegions/LCSlicer.h>
45 : #include <casacore/casa/Arrays/IPosition.h>
46 :
47 : #include <casacore/measures/Measures.h>
48 : #include <casacore/measures/Measures/MeasConvert.h>
49 :
50 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
51 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
52 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
53 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
54 : #include <casacore/coordinates/Coordinates/Projection.h>
55 : #include <casacore/coordinates/Coordinates/CoordinateUtil.h>
56 :
57 : #include <casacore/casa/BasicSL/String.h>
58 : #include <casacore/casa/Utilities/Assert.h>
59 : #include <casacore/casa/Exceptions/Error.h>
60 :
61 :
62 :
63 : using namespace casacore;
64 : namespace casa { //# NAMESPACE CASA - BEGIN
65 :
66 0 : PBMath1D::PBMath1D()
67 0 : : composite_p(2048)
68 : {
69 0 : };
70 :
71 :
72 22 : PBMath1D::PBMath1D(Quantity maximumRadius,
73 : Quantity refFreq,
74 : Bool isThisVP,
75 : BeamSquint squint,
76 22 : Bool useSymmetricBeam) :
77 : PBMathInterface(isThisVP, squint, useSymmetricBeam),
78 22 : wideFit_p(false),maximumRadius_p(maximumRadius),
79 22 : refFreq_p(refFreq),
80 44 : composite_p(2048)
81 : {
82 22 : fScale_p = refFreq_p.getValue("GHz"); // scale is ratio of refFreq_p to 1GHz
83 22 : refFreq_p = Quantity( 1.0, "GHz"); // internal Ref Freq is now 1GHz
84 : // convert instantiation parameters to GHz*arcmin reference
85 22 : maximumRadius_p = maximumRadius_p * fScale_p;
86 22 : scale_p = 1.0/(C::arcmin * C::giga);
87 22 : };
88 :
89 22 : PBMath1D::~PBMath1D()
90 : {
91 22 : };
92 :
93 :
94 :
95 :
96 : ImageRegion*
97 0 : PBMath1D::extent (const ImageInterface<Complex>& in, const MDirection& pointDir,
98 : const Int row, const Float fPad, const Int iChan,
99 : const SkyJones::SizeType sizeType)
100 : {
101 : if (row) {} // Not used yet
102 :
103 0 : CoordinateSystem coords=in.coordinates();
104 :
105 0 : Vector<Float> blc(4);
106 0 : Vector<Float> trc(4);
107 0 : blc.set(0.0);
108 0 : trc.set(0.0);
109 : {
110 : Int stokesIndex, k1, k2;
111 0 : CoordinateUtil::findStokesAxis(stokesIndex, k1, k2, coords);
112 0 : blc(stokesIndex) = 0.0;
113 0 : trc(stokesIndex) = in.shape()(stokesIndex)-1;
114 0 : Int spectralIndex=CoordinateUtil::findSpectralAxis(coords);
115 0 : blc(spectralIndex) = 0.0;
116 0 : trc(spectralIndex) = in.shape()(spectralIndex)-1;
117 : }
118 0 : extentguts(coords, pointDir, fPad, iChan, blc, trc);
119 0 : refineSize(blc, trc, in.shape(), sizeType);
120 0 : LCSlicer lcs( blc, trc );
121 0 : return ( new ImageRegion(lcs) );
122 0 : };
123 : ImageRegion*
124 0 : PBMath1D::extent (const ImageInterface<Float>& in, const MDirection& pointDir,
125 : const Int row, const Float fPad, const Int iChan,
126 : const SkyJones::SizeType sizeType)
127 : {
128 : if (row) {} // unused
129 0 : CoordinateSystem coords=in.coordinates();
130 0 : Vector<Float> blc(4);
131 0 : Vector<Float> trc(4);
132 0 : blc.set(0.0);
133 0 : trc.set(0.0);
134 : {
135 : Int stokesIndex, k1, k2;
136 0 : CoordinateUtil::findStokesAxis(stokesIndex, k1, k2, coords);
137 0 : blc(stokesIndex) = 0.0;
138 0 : trc(stokesIndex) = in.shape()(stokesIndex)-1;
139 0 : Int spectralIndex=CoordinateUtil::findSpectralAxis(coords);
140 0 : blc(spectralIndex) = 0.0;
141 0 : trc(spectralIndex) = in.shape()(spectralIndex)-1;
142 : }
143 0 : extentguts(coords, pointDir, fPad, iChan, blc, trc);
144 0 : refineSize(blc, trc, in.shape(), sizeType);
145 0 : LCSlicer lcs( blc, trc );
146 0 : return ( new ImageRegion(lcs) );
147 0 : };
148 :
149 :
150 :
151 8 : Int PBMath1D::support(const CoordinateSystem& cs){
152 8 : Int directionIndex=cs.findCoordinate(Coordinate::DIRECTION);
153 8 : AlwaysAssert(directionIndex>=0, AipsError);
154 : DirectionCoordinate
155 8 : directionCoord=cs.directionCoordinate(directionIndex);
156 :
157 8 : Vector<String> dirunit=directionCoord.worldAxisUnits();
158 :
159 : Double freq;
160 : {
161 8 : Int spectralIndex=cs.findCoordinate(Coordinate::SPECTRAL);
162 8 : AlwaysAssert(spectralIndex>=0, AipsError);
163 : SpectralCoordinate
164 8 : spectralCoord=cs.spectralCoordinate(spectralIndex);
165 :
166 :
167 8 : Vector<String> units(1);
168 8 : units = "Hz";
169 8 : spectralCoord.setWorldAxisUnits(units);
170 :
171 8 : Vector<Double> spectralWorld(1);
172 8 : Vector<Double> spectralPixel(1);
173 8 : spectralPixel(0) = 0;
174 8 : spectralCoord.toWorld(spectralWorld, spectralPixel);
175 8 : freq = spectralWorld(0);
176 8 : }
177 :
178 :
179 :
180 : // maximumRadius_p: maximum radius at 1 GHz frequency
181 : //Double delta = maximumRadius_p.getValue("rad") * 1.0e+9 / freq;
182 :
183 :
184 : //Number of pix at freq
185 8 : Double numpix=maximumRadius_p.getValue(dirunit(0))/fabs(directionCoord.increment()(0))*2.0*1.0e9/freq ;
186 :
187 :
188 8 : return Int(floor(numpix));
189 :
190 :
191 8 : }
192 0 : void PBMath1D::refineSize(Vector<Float>& blc, Vector<Float>& trc, const IPosition& shape,
193 : SkyJones::SizeType sizeType)
194 : {
195 : // Round Down and Up for BLC and TRC, make them integers
196 0 : Vector<Bool> blcTrouble(blc.nelements(), false);
197 0 : Vector<Bool> trcTrouble(blc.nelements(), false);
198 0 : Vector<Float> d1(2);
199 0 : Vector<Float> d2(2);
200 :
201 0 : for (Int i=0; i<2; i++) {
202 :
203 0 : blc(i) = (Int)(blc(i));
204 0 : trc(i) = (Int)(trc(i)+0.99); // OK, its ALMOST rounding up
205 :
206 0 : if (blc(i) < 0) {
207 0 : blc(i) = 0;
208 0 : blcTrouble(i) = true;
209 : }
210 0 : if (trc(i) > shape(i)-1) {
211 0 : trc(i) = shape(i)-1;
212 0 : trcTrouble(i) = true;
213 : }
214 :
215 0 : d1(i) = trc(i) - blc(i) + 1;
216 :
217 0 : if (sizeType == SkyJones::POWEROF2) {
218 0 : d2(i) = (Int)( pow( 2.0, (Double)(Int)(log((Double)d1(i))/log(2.0) + 1.0) )+0.01);
219 0 : } else if (sizeType == SkyJones::COMPOSITE) {
220 0 : d2(i) = composite_p.nextLarger( (Int)d1(i) );
221 : } else {
222 0 : d2(i) = d1(i);
223 : }
224 :
225 : // Deal with cases:
226 :
227 0 : if (d2(i) >= shape(i)) {
228 : // requested size doesn't even fit into image:
229 : // ----- revert to image size
230 0 : blc(i) = 0; trc(i) = shape(i)-1;
231 :
232 0 : } else if (blcTrouble(i)) {
233 : // requseted size fits, but buts up against the "bottom";
234 : // ----- make full adjustment to the "top"
235 0 : blc(i) = 0; trc(i) = d2(i)-1;
236 :
237 0 : } else if (trcTrouble(i)) {
238 : // requseted size fits, but buts up against the "top";
239 : // ----- make full adjustment to the "bottom"
240 0 : trc(i) = shape(i)-1; blc(i) = shape(i) - d2(i);
241 :
242 : } else {
243 : // requested subimage does not exceed starting image
244 : // ----- do appropriate thing, based on even or odd
245 0 : Float diff = d2(i) - d1(i);
246 0 : Bool even = (Bool)( (Int)diff == 2 * (Int)(diff/2) );
247 0 : if (even) {
248 0 : blc(i) = blc(i) - diff/2;
249 0 : trc(i) = trc(i) + diff/2;
250 : } else {
251 0 : blc(i) = blc(i) - diff/2 + 0.5;
252 0 : trc(i) = trc(i) + diff/2 + 0.5;
253 : }
254 : }
255 : }
256 0 : };
257 :
258 :
259 :
260 :
261 : void
262 0 : PBMath1D::extentguts (const CoordinateSystem& coords, const MDirection& pointDir,
263 : const Float fPad, const Int iChan, Vector<Float>& blc, Vector<Float>& trc)
264 :
265 : {
266 0 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
267 0 : AlwaysAssert(directionIndex>=0, AipsError);
268 : DirectionCoordinate
269 0 : directionCoord=coords.directionCoordinate(directionIndex);
270 0 : Vector<String> units(2); units = "deg";
271 0 : directionCoord.setWorldAxisUnits(units);
272 :
273 : // convert to the EPOCH of these coords
274 0 : MDirection::Types t2 = (MDirection::Types) (pointDir.getRef().getType());
275 0 : MDirection pointDirE;
276 0 : if (t2 != directionCoord.directionType()) {
277 0 : MDirection::Convert converter;
278 0 : ObsInfo oi=coords.obsInfo();
279 0 : converter.setOut(MDirection::Ref(directionCoord.directionType(),
280 0 : MeasFrame(oi.obsDate(), oi.telescopePosition())));
281 0 : pointDirE = converter(pointDir);
282 0 : } else {
283 0 : pointDirE = pointDir;
284 : }
285 :
286 : Double freq;
287 : {
288 0 : Int spectralIndex=coords.findCoordinate(Coordinate::SPECTRAL);
289 0 : AlwaysAssert(spectralIndex>=0, AipsError);
290 : SpectralCoordinate
291 0 : spectralCoord=coords.spectralCoordinate(spectralIndex);
292 :
293 0 : units.resize(1);
294 0 : units = "Hz";
295 0 : spectralCoord.setWorldAxisUnits(units);
296 :
297 0 : Vector<Double> spectralWorld(1);
298 0 : Vector<Double> spectralPixel(1);
299 0 : spectralPixel(0) = iChan;
300 0 : spectralCoord.toWorld(spectralWorld, spectralPixel);
301 0 : freq = spectralWorld(0);
302 0 : }
303 :
304 0 : Vector<Double> edgeWorld(2);
305 0 : Vector<Double> edge1Pixel(2);
306 0 : Vector<Double> edge2Pixel(2);
307 :
308 :
309 : // maximumRadius_p: maximum radius at 1 GHz frequency
310 0 : Double delta = maximumRadius_p.getValue("rad") * 1.0e+9 / freq;
311 : {
312 0 : MDirection edgeDir( pointDirE );
313 0 : edgeDir.shift( delta, 0.0, true);
314 0 : edgeWorld(0) = edgeDir.getAngle().getValue("deg")(0);
315 0 : edgeWorld(1) = edgeDir.getAngle().getValue("deg")(1);
316 0 : directionCoord.toPixel(edge1Pixel, edgeWorld);
317 0 : }
318 : {
319 0 : MDirection edgeDir( pointDirE );
320 0 : edgeDir.shift( -delta, 0.0, true);
321 0 : edgeWorld(0) = edgeDir.getAngle().getValue("deg")(0);
322 0 : edgeWorld(1) = edgeDir.getAngle().getValue("deg")(1);
323 0 : directionCoord.toPixel(edge2Pixel, edgeWorld);
324 0 : }
325 0 : blc(0) = min( edge1Pixel(0), edge2Pixel(0) );
326 0 : trc(0) = max( edge1Pixel(0), edge2Pixel(0) );
327 0 : if (fPad > 0.1) {
328 0 : Float pad = (trc(0) - blc(0)) * (fPad - 1.0)/2;
329 0 : blc(0) = blc(0) - pad;
330 0 : trc(0) = trc(0) + pad;
331 : }
332 : {
333 0 : MDirection edgeDir( pointDirE );
334 0 : edgeDir.shift( 0.0, delta, true);
335 0 : edgeWorld(0) = edgeDir.getAngle().getValue("deg")(0);
336 0 : edgeWorld(1) = edgeDir.getAngle().getValue("deg")(1);
337 0 : directionCoord.toPixel(edge1Pixel, edgeWorld);
338 0 : }
339 : {
340 0 : MDirection edgeDir( pointDirE );
341 0 : edgeDir.shift( 0.0, -delta, true);
342 0 : edgeWorld(0) = edgeDir.getAngle().getValue("deg")(0);
343 0 : edgeWorld(1) = edgeDir.getAngle().getValue("deg")(1);
344 0 : directionCoord.toPixel(edge2Pixel, edgeWorld);
345 0 : }
346 0 : blc(1) = min( edge1Pixel(1), edge2Pixel(1) );
347 0 : trc(1) = max( edge1Pixel(1), edge2Pixel(1) );
348 0 : if (fPad > 0.1) {
349 0 : Float pad = (trc(1) - blc(1)) * (fPad - 1.0)/2;
350 0 : blc(1) = blc(1) - pad;
351 0 : trc(1) = trc(1) + pad;
352 : }
353 0 : };
354 :
355 :
356 :
357 :
358 0 : void PBMath1D::symmetrizeSquintedBeam()
359 : {
360 : // eventually we need to create the 2D squinted RR and LL
361 : // beams and average them. For now, we just return the
362 : // unsquinted beams
363 :
364 0 : if (vp_p.nelements() == 0) {
365 0 : fillPBArray();
366 : }
367 0 : esvp_p = vp_p;
368 0 : };
369 :
370 : ImageInterface<Complex>&
371 48 : PBMath1D::apply(const ImageInterface<Complex>& in,
372 : ImageInterface<Complex>& out,
373 : const MDirection& pointDir,
374 : const Quantity parAngle,
375 : const BeamSquint::SquintType doSquint,
376 : Bool inverse,
377 : Bool conjugate,
378 : Int iPower,
379 : Float cutoff,
380 : Bool forward)
381 : {
382 96 : LogIO os(LogOrigin("PBMath1D", "apply"));
383 : // Check that in and out are comparable:
384 48 : if (in.shape() != out.shape()) {
385 0 : throw(AipsError("PBMath1D::apply(ImageInterface...) - in and out images have different shapes"));
386 : }
387 48 : CoordinateSystem coords=in.coordinates();
388 48 : if (!coords.near(out.coordinates()) ) {
389 0 : throw(AipsError("PBMath1D::apply(ImageInterface...) - in and out images have different coordinates"));
390 : }
391 :
392 48 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
393 48 : AlwaysAssert(directionIndex>=0, AipsError);
394 : DirectionCoordinate
395 48 : directionCoord=coords.directionCoordinate(directionIndex);
396 48 : Vector<String> units(2); units = "deg";
397 48 : directionCoord.setWorldAxisUnits(units);
398 :
399 : // convert to the EPOCH of these coords
400 48 : MDirection::Types t2 = (MDirection::Types) (pointDir.getRef().getType());
401 48 : MDirection pointDirE;
402 :
403 :
404 48 : if (t2 != directionCoord.directionType()) {
405 0 : MDirection::Convert converter;
406 0 : ObsInfo oi=coords.obsInfo();
407 0 : converter.setOut(MDirection::Ref(directionCoord.directionType(),
408 0 : MeasFrame(oi.obsDate(), oi.telescopePosition())));
409 0 : pointDirE = converter(pointDir);
410 0 : } else {
411 48 : pointDirE = pointDir;
412 : }
413 :
414 48 : Int stokesIndex=coords.findCoordinate(Coordinate::STOKES);
415 48 : AlwaysAssert(stokesIndex>=0, AipsError);
416 : StokesCoordinate
417 48 : stokesCoord=coords.stokesCoordinate(stokesIndex);
418 :
419 48 : Int spectralIndex=coords.findCoordinate(Coordinate::SPECTRAL);
420 48 : AlwaysAssert(spectralIndex>=0, AipsError);
421 : SpectralCoordinate
422 48 : spectralCoord=coords.spectralCoordinate(spectralIndex);
423 :
424 48 : units.resize(1);
425 48 : units = "Hz";
426 48 : spectralCoord.setWorldAxisUnits(units);
427 :
428 48 : Int nchan=in.shape()(3);
429 :
430 48 : Vector<Double> pointingCenterWorld(2);
431 48 : Vector<Double> pointingCenterPixel(2);
432 48 : Vector<Double> directionPixel(2);
433 :
434 48 : pointingCenterWorld(0) = pointDirE.getAngle().getValue("deg")(0);
435 48 : pointingCenterWorld(1) = pointDirE.getAngle().getValue("deg")(1);
436 48 : directionCoord.toPixel(pointingCenterPixel, pointingCenterWorld);
437 48 : MDirection newpointDirE;
438 48 : Vector<Double> nonSquintedPointingPixel = pointingCenterPixel.copy();
439 :
440 48 : os << "pointingCenterWorld " << pointingCenterWorld << LogIO::DEBUGGING;
441 48 : os << "pointingCenterPixel " << pointingCenterPixel << LogIO::DEBUGGING;
442 :
443 : // Fill in a cache of the frequencies & squints
444 48 : Vector<Double> spectralWorld(1);
445 48 : Vector<Double> spectralPixel(1);
446 48 : Matrix<Double> xSquintPixCache(2, nchan);
447 48 : Matrix<Double> ySquintPixCache(2, nchan);
448 48 : Vector<Double> spectralCache(nchan);
449 :
450 : {
451 96 : for(Int chan=0;chan<nchan;chan++) {
452 48 : spectralPixel(0)=chan;
453 48 : if(!spectralCoord.toWorld(spectralWorld, spectralPixel)) {
454 0 : os << "Cannot find frequency for this plane" << LogIO::EXCEPTION;
455 : }
456 48 : spectralCache(chan)=spectralWorld(0);
457 :
458 :
459 48 : if (doSquint == BeamSquint::RR || doSquint == BeamSquint::GOFIGURE) {
460 40 : squint_p.getPointingDirection (pointDirE,
461 : parAngle,
462 80 : Quantity(spectralWorld(0),"Hz"),
463 : BeamSquint::RR, newpointDirE);
464 40 : pointingCenterWorld(0) = newpointDirE.getAngle().getValue("deg")(0);
465 40 : pointingCenterWorld(1) = newpointDirE.getAngle().getValue("deg")(1);
466 40 : directionCoord.toPixel(pointingCenterPixel, pointingCenterWorld);
467 40 : xSquintPixCache(0, chan) = pointingCenterPixel(0);
468 40 : ySquintPixCache(0, chan) = pointingCenterPixel(1);
469 : } else {
470 8 : xSquintPixCache(0, chan) = nonSquintedPointingPixel(0);
471 8 : ySquintPixCache(0, chan) = nonSquintedPointingPixel(1);
472 : }
473 48 : if (doSquint == BeamSquint::LL || doSquint == BeamSquint::GOFIGURE) {
474 40 : squint_p.getPointingDirection (pointDirE,
475 : parAngle,
476 80 : Quantity(spectralWorld(0),"Hz"),
477 : BeamSquint::LL, newpointDirE);
478 40 : pointingCenterWorld(0) = newpointDirE.getAngle().getValue("deg")(0);
479 40 : pointingCenterWorld(1) = newpointDirE.getAngle().getValue("deg")(1);
480 40 : directionCoord.toPixel(pointingCenterPixel, pointingCenterWorld);
481 40 : xSquintPixCache(1, chan) = pointingCenterPixel(0);
482 40 : ySquintPixCache(1, chan) = pointingCenterPixel(1);
483 : } else {
484 8 : xSquintPixCache(1, chan) = nonSquintedPointingPixel(0);
485 8 : ySquintPixCache(1, chan) = nonSquintedPointingPixel(1);
486 : }
487 : }
488 : }
489 :
490 : /*
491 : cout << "pointingCenterPixel x,y = " << nonSquintedPointingPixel << endl;
492 : cout << "squinted pointingCenterPixel x,y RR = " << xSquintPixCache(0, 0) << ", "
493 : << ySquintPixCache(0, 0) << endl;
494 : cout << "squinted pointingCenterPixel x,y LL = " << xSquintPixCache(1, 0) << ", "
495 : << ySquintPixCache(1, 0) << endl;
496 : */
497 :
498 : // Iterate through in minimum IO/Memory chunks
499 : //IPosition ncs = in.niceCursorShape();
500 48 : IPosition ncs=in.shape();
501 48 : ncs(2) = 1; ncs(3) = 1;
502 96 : RO_LatticeIterator<Complex> li(in, LatticeStepper(in.shape(), ncs, IPosition(4,0,1,2,3) ) );
503 96 : LatticeIterator<Complex> oli(out, LatticeStepper(in.shape(), ncs, IPosition(4,0,1,2,3)) );
504 :
505 : // taper no longer appears to be used
506 : //Complex taper;
507 : //Float r2=0.0;
508 : //Float r=0.0;
509 :
510 48 : Vector<Double> increment = directionCoord.increment();
511 48 : Int rrplane = -1;
512 48 : Int llplane = -1;
513 48 : stokesCoord.toPixel( rrplane, Stokes::RR );
514 48 : stokesCoord.toPixel( llplane, Stokes::LL );
515 :
516 : /*
517 : cout << "stokes types in image = " << stokesCoord.stokes() << endl;
518 : cout << "rr plane = " << rrplane << " ll plane = " << llplane << endl;
519 :
520 : */
521 : Double xPixel; Double yPixel;
522 :
523 48 : Int laststokes = -1;
524 48 : Int lastChan = -1;
525 : Int ichan;
526 : Int istokes;
527 : Int ix0, iy0;
528 : //Int indx;
529 96 : for(li.reset(),oli.reset();!li.atEnd();li++,oli++) {
530 :
531 48 : IPosition itsShape(li.matrixCursor().shape());
532 48 : IPosition loc(li.position());
533 :
534 48 : ichan = loc(3);
535 48 : istokes = loc(2);
536 48 : iy0 = loc(1);
537 48 : ix0 = loc(0);
538 :
539 : // determine the pointing: RR, LL, or Center? We make a slight mistake
540 : // here since we ignore the difference between the RR beam and the
541 : // RL beam, say. The latter is slightly smaller because of the
542 : // squint. Hence this code should be deprecated in favor of the
543 : // correct 2D version (when mosaicing in polarization)
544 48 : if ((doSquint == BeamSquint::RR) ||
545 40 : ((doSquint == BeamSquint::GOFIGURE) && (istokes == rrplane)) ) {
546 0 : xPixel = xSquintPixCache(0, ichan);
547 0 : yPixel = ySquintPixCache(0, ichan);
548 48 : } else if ((doSquint == BeamSquint::LL) ||
549 40 : ((doSquint == BeamSquint::GOFIGURE) && (istokes == llplane)) ) {
550 0 : xPixel = xSquintPixCache(1, ichan);
551 0 : yPixel = ySquintPixCache(1, ichan);
552 : } else {
553 48 : xPixel = nonSquintedPointingPixel(0);
554 48 : yPixel = nonSquintedPointingPixel(1);
555 : }
556 :
557 48 : if (istokes != laststokes) {
558 : //cerr << "Stokes = " << istokes << " pix = " << xPixel << ", " << yPixel << endl;
559 48 : laststokes = istokes;
560 : }
561 :
562 48 : Double factor = 60.0 * spectralCache(ichan)/1.0e+9 ; // arcminutes * GHz
563 48 : Double rmax2 = square( maximumRadius_p.getValue("'") / factor );
564 48 : if (wideFit_p) {
565 : // fill vp with interpolated values for current frequency
566 0 : if (ichan!=lastChan) {
567 0 : nearestVPArray(spectralCache(ichan));
568 0 : lastChan=ichan;
569 : }
570 : }
571 :
572 48 : Vector<Float> rx2(itsShape(0));
573 48 : Vector<Float> ry2(itsShape(1));
574 324568 : for(Int ix=0;ix<itsShape(0);ix++) {
575 324520 : rx2(ix) = square( increment(0)*((Double)(ix+ix0) - xPixel) );
576 : }
577 74984 : for(Int iy=0;iy<itsShape(1);iy++) {
578 74936 : ry2(iy) = square( increment(1)*((Double)(iy+iy0) - yPixel) );
579 : }
580 :
581 48 : const Matrix<Complex>& inmat = li.matrixCursor();
582 48 : Matrix<Complex>& outmat=oli.rwMatrixCursor();
583 :
584 : Bool incopy, outcopy, del;
585 48 : const Complex * inpoint = inmat.getStorage(incopy);
586 48 : Complex *outpoint =outmat.getStorage(outcopy);
587 48 : Float * rx2point = rx2.getStorage(del);
588 48 : Float * ry2point= ry2.getStorage(del);
589 48 : Complex* vppoint=vp_p.getStorage(del);
590 48 : Int nx=itsShape(0);
591 48 : Int ny=itsShape(1);
592 48 : Double inverseIncrementRadius=inverseIncrementRadius_p;
593 48 : #pragma omp parallel default(none) firstprivate(inpoint, outpoint, rx2point, ry2point, vppoint, iPower, conjugate, inverse, forward, nx, ny, rmax2, factor, inverseIncrementRadius, cutoff)
594 : {
595 : #pragma omp for
596 : for(Int iy=0;iy<ny;iy++) {
597 : Float ry2val=ry2point[iy];
598 : applyXLine(inpoint, outpoint , rx2point , vppoint , ry2val, iPower, conjugate, inverse, forward, nx, iy, rmax2,
599 : factor, inverseIncrementRadius, cutoff);
600 : /*for(Int ix=0;ix<itsShape(0);ix++) {
601 :
602 : r2 = rx2(ix) + ry2(iy);
603 :
604 : if (r2 > rmax2) {
605 : oli.rwMatrixCursor()(ix, iy) = 0.0;
606 : } else {
607 : r = sqrt(r2) * factor;
608 : indx = Int(r*inverseIncrementRadius_p);
609 : if (norm(vp_p(indx)) > 0.0) {
610 : if(iPower==2) {
611 : taper = vp_p(indx) * conj(vp_p(indx));
612 : }
613 : else {
614 : taper = vp_p(indx);
615 : }
616 : } else {
617 : taper = 0.0;
618 : }
619 : if (conjugate) {
620 : taper = conj(taper);
621 : }
622 : // Differentiate between forward (Sky->UV) and
623 : // inverse (UV->Sky) - these need different
624 : // applications of the PB
625 : if(!forward) {
626 : taper = conj(taper);
627 : }
628 : if (inverse) {
629 : if (abs(taper) < cutoff ) {
630 : oli.rwMatrixCursor()(ix, iy) = 0.0;
631 : } else {
632 : oli.rwMatrixCursor()(ix, iy) = li.matrixCursor()(ix, iy) / taper ;
633 : }
634 : } else { // not inverse!
635 : oli.rwMatrixCursor()(ix, iy) = li.matrixCursor()(ix, iy) * taper;
636 : }
637 : }
638 : }
639 : */
640 : }
641 : } //end pragma
642 48 : outmat.putStorage(outpoint, outcopy);
643 48 : inmat.freeStorage(inpoint, incopy);
644 48 : }
645 :
646 48 : return out;
647 :
648 48 : };
649 :
650 74936 : void PBMath1D::applyXLine(const Complex*& in, Complex*& out, Float*& rx2, Complex*& vp, const Float ry2, const Int ipower, const Bool conjugate, const Bool inverse, const Bool forward, const Int nx, const Int iy, const Double rmax2, const Double factor, const Double inverseIncrementRadius, const Float cutoff)
651 : {
652 : Float r;
653 : Int indx;
654 74936 : Complex taper;
655 140899296 : for(Int ix=0;ix<nx;ix++) {
656 :
657 140824360 : Float r2 = rx2[ix] + ry2;
658 :
659 140824360 : if (r2 > rmax2){
660 105274055 : out[ix+iy*nx] = 0.0;
661 : }
662 : else {
663 35550305 : r = sqrt(r2) * factor;
664 35550305 : indx = Int(r*inverseIncrementRadius);
665 35550305 : if (norm(vp[indx]) > 0.0) {
666 35550305 : if(ipower==2) {
667 35550305 : taper = vp[indx] * conj(vp[indx]);
668 : }
669 : else {
670 0 : taper = vp[indx];
671 : }
672 : } else {
673 0 : taper = 0.0;
674 : }
675 35550305 : if (conjugate) {
676 0 : taper = conj(taper);
677 : }
678 : // Differentiate between forward (Sky->UV) and
679 : // inverse (UV->Sky) - these need different
680 : // applications of the PB
681 35550305 : if(!forward) {
682 0 : taper = conj(taper);
683 : }
684 35550305 : if (inverse) {
685 0 : if (abs(taper) < cutoff ) {
686 0 : out[ix+iy*nx] = 0.0;
687 : } else {
688 0 : out[ix+iy*nx] = (in[ix+iy*nx]) / taper ;
689 : }
690 : } else { // not inverse!
691 35550305 : out[ix+iy*nx] = (in[ix+iy*nx]) * taper ;
692 : }
693 : }
694 : }
695 74936 : };
696 :
697 : ImageInterface<Float>&
698 0 : PBMath1D::apply(const ImageInterface<Float>& in,
699 : ImageInterface<Float>& out,
700 : const MDirection& pointDir,
701 : const Quantity parAngle,
702 : const BeamSquint::SquintType doSquint,
703 : Float /*cutoff*/, const Int ipower)
704 : {
705 0 : LogIO os(LogOrigin("PBMath1D", "apply"));
706 :
707 : // cout << "PBMath1D::apply: image shape: " << in.shape() << endl;
708 : // Check that in and out are comparable:
709 0 : if (in.shape() != out.shape()) {
710 0 : throw(AipsError("PBMath1D::apply(ImageInterface...) - in and out images have different shapes"));
711 :
712 : }
713 0 : CoordinateSystem coords=in.coordinates();
714 0 : if (!coords.near(out.coordinates())) {
715 0 : throw(AipsError("PBMath1D::apply(ImageInterface...) - in and out images have different coordinates"));
716 : }
717 :
718 0 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
719 0 : AlwaysAssert(directionIndex>=0, AipsError);
720 : DirectionCoordinate
721 0 : directionCoord=coords.directionCoordinate(directionIndex);
722 0 : Vector<String> units(2); units = "deg";
723 0 : directionCoord.setWorldAxisUnits(units);
724 :
725 : // convert to the EPOCH of these coords
726 0 : MDirection::Types t2 = (MDirection::Types) (pointDir.getRef().getType());
727 0 : MDirection pointDirE(pointDir);
728 :
729 0 : if (t2 != directionCoord.directionType()) {
730 0 : MDirection::Convert converter;
731 0 : ObsInfo oi=coords.obsInfo();
732 0 : converter.setOut(MDirection::Ref(directionCoord.directionType(),
733 0 : MeasFrame(oi.obsDate(), oi.telescopePosition())));
734 0 : pointDirE = converter(pointDir);
735 0 : } else {
736 0 : pointDirE = pointDir;
737 : }
738 0 : Int stokesIndex=coords.findCoordinate(Coordinate::STOKES);
739 0 : AlwaysAssert(stokesIndex>=0, AipsError);
740 : StokesCoordinate
741 0 : stokesCoord=coords.stokesCoordinate(stokesIndex);
742 :
743 0 : Int spectralIndex=coords.findCoordinate(Coordinate::SPECTRAL);
744 0 : AlwaysAssert(spectralIndex>=0, AipsError);
745 : SpectralCoordinate
746 0 : spectralCoord=coords.spectralCoordinate(spectralIndex);
747 :
748 0 : units.resize(1);
749 0 : units = "Hz";
750 0 : spectralCoord.setWorldAxisUnits(units);
751 :
752 0 : Int nchan=in.shape()(3);
753 :
754 0 : Vector<Double> pointingCenterWorld(2);
755 0 : Vector<Double> pointingCenterPixel(2);
756 0 : Vector<Double> directionPixel(2);
757 :
758 0 : pointingCenterWorld(0) = pointDirE.getAngle().getValue("deg")(0);
759 0 : pointingCenterWorld(1) = pointDirE.getAngle().getValue("deg")(1);
760 0 : directionCoord.toPixel(pointingCenterPixel, pointingCenterWorld);
761 0 : MDirection newpointDirE;
762 0 : Vector<Double> nonSquintedPointingPixel = pointingCenterPixel.copy();
763 :
764 0 : os << "pointingCenterWorld " << pointingCenterWorld << LogIO::DEBUGGING;
765 0 : os << "pointingCenterPixel " << pointingCenterPixel << LogIO::DEBUGGING;
766 :
767 : // Fill in a cache of the frequencies & squints
768 0 : Vector<Double> spectralWorld(1);
769 0 : Vector<Double> spectralPixel(1);
770 0 : Matrix<Double> xSquintPixCache(2, nchan); // kludge: prevent errors when nchan = 1
771 0 : Matrix<Double> ySquintPixCache(2, nchan);
772 0 : Vector<Double> spectralCache(nchan);
773 :
774 : {
775 0 : for(Int chan=0;chan<nchan;chan++) {
776 0 : spectralPixel(0)=chan;
777 0 : if(!spectralCoord.toWorld(spectralWorld, spectralPixel)) {
778 0 : os << "Cannot find frequency for this plane" << LogIO::EXCEPTION;
779 : }
780 0 : spectralCache(chan)=spectralWorld(0);
781 :
782 :
783 0 : if (doSquint == BeamSquint::RR || doSquint == BeamSquint::GOFIGURE) {
784 0 : squint_p.getPointingDirection (pointDirE,
785 : parAngle,
786 0 : Quantity(spectralWorld(0),"Hz"),
787 : BeamSquint::RR, newpointDirE);
788 0 : pointingCenterWorld(0) = newpointDirE.getAngle().getValue("deg")(0);
789 0 : pointingCenterWorld(1) = newpointDirE.getAngle().getValue("deg")(1);
790 0 : directionCoord.toPixel(pointingCenterPixel, pointingCenterWorld);
791 0 : xSquintPixCache(0, chan) = pointingCenterPixel(0);
792 0 : ySquintPixCache(0, chan) = pointingCenterPixel(1);
793 : } else {
794 0 : xSquintPixCache(0, chan) = nonSquintedPointingPixel(0);
795 0 : ySquintPixCache(0, chan) = nonSquintedPointingPixel(1);
796 : }
797 0 : if (doSquint == BeamSquint::LL || doSquint == BeamSquint::GOFIGURE) {
798 0 : squint_p.getPointingDirection (pointDirE,
799 : parAngle,
800 0 : Quantity(spectralWorld(0),"Hz"),
801 : BeamSquint::LL, newpointDirE);
802 0 : pointingCenterWorld(0) = newpointDirE.getAngle().getValue("deg")(0);
803 0 : pointingCenterWorld(1) = newpointDirE.getAngle().getValue("deg")(1);
804 0 : directionCoord.toPixel(pointingCenterPixel, pointingCenterWorld);
805 0 : xSquintPixCache(1, chan) = pointingCenterPixel(0);
806 0 : ySquintPixCache(1, chan) = pointingCenterPixel(1);
807 : } else {
808 0 : xSquintPixCache(1, chan) = nonSquintedPointingPixel(0);
809 0 : ySquintPixCache(1, chan) = nonSquintedPointingPixel(1);
810 : }
811 : }
812 : }
813 :
814 :
815 : // Iterate through in minimum IO/Memory chunks
816 0 : IPosition ncs = in.niceCursorShape();
817 : // IPosition ncs=in.shape();
818 0 : ncs(2) = 1; ncs(3) = 1;
819 0 : ncs(2) = 1; ncs(3) = 1;
820 0 : RO_LatticeIterator<Float> li(in, LatticeStepper(in.shape(), ncs, IPosition(4,0,1,2,3) ) );
821 0 : LatticeIterator<Float> oli(out, LatticeStepper(in.shape(), ncs, IPosition(4,0,1,2,3)) );
822 :
823 : //Vector<Float> taper(vp_p.nelements());
824 0 : Float taper=0.0;
825 0 : Float r2=0.0;
826 0 : Float r=0.0;
827 :
828 0 : Vector<Double> increment = directionCoord.increment();
829 0 : Int rrplane = -1;
830 0 : Int llplane = -1;
831 0 : stokesCoord.toPixel( rrplane, Stokes::RR );
832 0 : stokesCoord.toPixel( llplane, Stokes::LL );
833 :
834 : Double xPixel; Double yPixel;
835 :
836 0 : Int laststokes = -1;
837 0 : Int lastChan = -1;
838 : Int ichan;
839 : Int istokes;
840 : Int ix0, iy0;
841 : Int indx;
842 :
843 0 : for(li.reset(),oli.reset();!li.atEnd();li++,oli++) {
844 :
845 0 : IPosition itsShape(li.matrixCursor().shape());
846 0 : IPosition loc(li.position());
847 :
848 0 : ichan = loc(3);
849 0 : istokes = loc(2);
850 0 : iy0 = loc(1);
851 0 : ix0 = loc(0);
852 :
853 : // determine the pointing: RR, LL, or Center?
854 0 : if ((doSquint == BeamSquint::RR) ||
855 0 : ((doSquint == BeamSquint::GOFIGURE) && (istokes == rrplane)) ) {
856 0 : xPixel = xSquintPixCache(0, ichan);
857 0 : yPixel = ySquintPixCache(0, ichan);
858 0 : } else if ((doSquint == BeamSquint::LL) ||
859 0 : ((doSquint == BeamSquint::GOFIGURE) && (istokes == llplane ))) {
860 0 : xPixel = xSquintPixCache(1, ichan);
861 0 : yPixel = ySquintPixCache(1, ichan);
862 : } else {
863 0 : xPixel = nonSquintedPointingPixel(0);
864 0 : yPixel = nonSquintedPointingPixel(1);
865 : }
866 :
867 0 : if (istokes != laststokes) {
868 : // cout << "Stokes = " << istokes << " pix = " << xPixel << ", " << yPixel << endl;
869 0 : laststokes = istokes;
870 : }
871 :
872 0 : Double factor = 60.0 * spectralCache(ichan)/1.0e+9 ; // arcminutes * GHz
873 0 : Double rmax2 = square( maximumRadius_p.getValue("'") / factor );
874 0 : if (wideFit_p) {
875 : // fill vp with interpolated values for current frequency
876 0 : if (ichan!=lastChan) {
877 0 : nearestVPArray(spectralCache(ichan));
878 0 : lastChan = ichan;
879 : }
880 : }
881 :
882 : //taper.set(0.0);
883 : /*for (uInt inda=0; inda < vp_p.nelements(); ++inda){
884 : if (norm(vp_p(inda)) > 0.0) {
885 : taper[inda] = real(vp_p(inda)) * real(vp_p(inda)) + imag(vp_p(inda))*imag(vp_p(inda));
886 : if(ipower==4)
887 : taper[inda] *= taper[inda];
888 : }
889 :
890 : }
891 : */
892 :
893 0 : Vector<Float> rx2(itsShape(0));
894 0 : Vector<Float> ry2(itsShape(1));
895 0 : for(Int ix=0;ix<itsShape(0);ix++) {
896 0 : rx2(ix) = square( increment(0)*((Double)(ix+ix0) - xPixel) );
897 : }
898 0 : for(Int iy=0;iy<itsShape(1);iy++) {
899 0 : ry2(iy) = square( increment(1)*((Double)(iy+iy0) - yPixel) );
900 : }
901 :
902 0 : for(Int iy=0;iy<itsShape(1);iy++) {
903 0 : for(Int ix=0;ix<itsShape(0);ix++) {
904 :
905 0 : r2 = rx2(ix) + ry2(iy);
906 0 : if (r2 > rmax2) {
907 0 : oli.rwMatrixCursor()(ix, iy) = 0.0;
908 : } else {
909 0 : r = sqrt(r2) * factor;
910 0 : indx = Int(r*inverseIncrementRadius_p);
911 0 : if (norm(vp_p(indx)) > 0.0) {
912 0 : taper = real(vp_p(indx)) * real(vp_p(indx))+ imag(vp_p(indx)) * imag(vp_p(indx));
913 0 : if(ipower==4)
914 0 : taper *= taper;
915 : }
916 : //else {
917 : // taper = 0.0;
918 : //}
919 0 : oli.rwMatrixCursor()(ix, iy) = li.matrixCursor()(ix,iy) * taper;
920 : }
921 : }
922 : }
923 0 : }
924 0 : return out;
925 :
926 0 : };
927 :
928 : // Behavior: doSquint == RR or LL don't make sense here
929 : //
930 : //
931 :
932 : SkyComponent&
933 2491 : PBMath1D::apply(SkyComponent& in,
934 : SkyComponent& out,
935 : const MDirection& pointDir,
936 : const Quantity frequency,
937 : const Quantity parAngle,
938 : const BeamSquint::SquintType doSquint,
939 : Bool inverse,
940 : Bool conjugate,
941 : Int iPower,
942 : Float cutoff,
943 : Bool /*forward*/)
944 : {
945 : // if ( doSquint == NONE ) we can deal with any polarisation representation
946 : // if ( doSquint == GOFIGURE) an exception is thrown if polarisation is not CIRCULAR
947 : // if ( doSquint == RR || doSquint == LL ) an exception is thrown,
948 : // as it is not valid to apply the RR or LL squint to ALL polarisations
949 :
950 : // Also: we can do nothing with spectral index models
951 :
952 :
953 : // convert to the EPOCH of these coords
954 2491 : MDirection::Types t1 = (MDirection::Types) (in.shape().refDirection().getRef().getType());
955 2491 : MDirection::Types t2 = (MDirection::Types) (pointDir.getRef().getType());
956 :
957 2491 : MDirection pointDirE;
958 2491 : if ( t1 != t2) {
959 0 : MDirection::Convert converter;
960 0 : converter.setOut( t1 );
961 0 : pointDirE = converter(pointDir);
962 0 : } else {
963 2491 : pointDirE = pointDir;
964 : }
965 :
966 2491 : if (doSquint == BeamSquint::RR || doSquint == BeamSquint::LL) {
967 0 : throw(AipsError("PBMath1D::apply(SkyComponent...) - cannot force a SkyComponent to have Squint RR or LL"));
968 : }
969 2491 : if (doSquint == BeamSquint::GOFIGURE) {
970 2491 : if (in.flux().pol() != ComponentType::CIRCULAR) {
971 397 : in.flux().convertPol(ComponentType::CIRCULAR);
972 : }
973 : }
974 :
975 2491 : Vector<DComplex> compFluxIn = in.flux().value();
976 2491 : Vector<DComplex> compFlux = out.flux().value();
977 2491 : compFlux = compFluxIn.copy();
978 :
979 : // Find the direction of the component
980 2491 : MDirection compDir = in.shape().refDirection();
981 :
982 : // Now taper all polarizations appropriately
983 :
984 : // Sort out any frequency interpolation
985 : //Int ifit=0;
986 : //Float lfit=0;
987 : //Int nFreq=wFreqs_p.nelements();
988 2491 : if (wideFit_p) {
989 0 : Double freq = frequency.getValue("Hz");
990 0 : nearestVPArray(freq);
991 : /* for (ifit=0; ifit<nFreq; ifit++) {
992 : if (freq<=wFreqs_p(ifit)) break;
993 : }
994 : if (ifit>0 && ifit<nFreq) {
995 : lfit=(freq-wFreqs_p(ifit-1)) / (wFreqs_p(ifit)-wFreqs_p(ifit-1));
996 : }
997 : */
998 : }
999 :
1000 2491 : MDirection newpointDirE;
1001 12455 : for (Int pol=0;pol<4;pol++) {
1002 9964 : Stokes::StokesTypes stokes=Stokes::type(pol+5);
1003 :
1004 9964 : if (stokes == Stokes::RR && doSquint == BeamSquint::GOFIGURE) {
1005 2491 : squint_p.getPointingDirection (pointDirE, parAngle, frequency, BeamSquint::RR,
1006 : newpointDirE );
1007 7473 : } else if (stokes == Stokes::LL && doSquint == BeamSquint::GOFIGURE) {
1008 2491 : squint_p.getPointingDirection (pointDirE, parAngle, frequency, BeamSquint::LL,
1009 : newpointDirE );
1010 : } else {
1011 4982 : newpointDirE = pointDirE;
1012 : }
1013 :
1014 9964 : MVDirection mvd1( compDir.getAngle() );
1015 9964 : MVDirection mvd2( newpointDirE.getAngle() );
1016 9964 : Quantity sep = mvd1.separation(mvd2, "'");
1017 9964 : double r = sep.getValue("'") * frequency.getValue("Hz") / 1.0e+9; // arcminutes * GHz
1018 9964 : Complex taper;
1019 9964 : Int ir = Int(r*inverseIncrementRadius_p);
1020 : //vp_p is interpolated wvp_p from above
1021 9964 : Complex vpVal = ir >= Int(vp_p.nelements()) ? Complex(0) : vp_p(ir);
1022 : /*if (wideFit_p) {
1023 : if (ifit==0) {
1024 : vpVal = wbvp_p(ir,0);
1025 : } else if (ifit==nFreq) {
1026 : vpVal = wbvp_p(ir,nFreq-1);
1027 : } else {
1028 : vpVal = wbvp_p(ir,ifit-1)*(1-lfit) + wbvp_p(ir,ifit)*lfit;
1029 : }
1030 : }
1031 : */
1032 :
1033 9964 : if (r > maximumRadius_p.getValue("'")) {
1034 192 : compFlux(pol) = 0.0;
1035 : } else {
1036 9772 : if (norm(vpVal) > 0.0) {
1037 9772 : if(iPower>1){
1038 9772 : taper=vpVal*conj(vpVal);
1039 9772 : if(iPower==4)
1040 0 : taper*=taper;
1041 : }
1042 : else{
1043 0 : taper = vpVal;
1044 : //taper = pow( vp_p(Int(r*inverseIncrementRadius_p)), (Float)iPower);
1045 : }
1046 : } else {
1047 0 : taper = 0.0;
1048 : }
1049 9772 : if (conjugate) {
1050 0 : taper = conj(taper);
1051 : }
1052 9772 : if (inverse) {
1053 0 : if (abs(taper) < cutoff ) {
1054 0 : compFlux(pol) = 0.0;
1055 : } else {
1056 0 : compFlux(pol) /= taper ;
1057 : }
1058 : } else { // not inverse!
1059 9772 : compFlux(pol) *= taper;
1060 : }
1061 : }
1062 9964 : }
1063 :
1064 : // Set the output component fluxes
1065 2491 : out = in.copy();
1066 2491 : out.flux().setValue(compFlux);
1067 :
1068 2491 : return out;
1069 :
1070 2491 : };
1071 :
1072 14 : void PBMath1D::summary(Int nValues)
1073 : {
1074 14 : String name;
1075 14 : namePBClass(name);
1076 28 : LogIO os(LogOrigin("PBMath1D", "summary"));
1077 14 : os << "Using " << name << " PB Class " << LogIO::POST;
1078 14 : PBMathInterface::summary(nValues);
1079 :
1080 14 : if (nValues > 0) {
1081 0 : os << "Primary Beam Sampled Data: " << LogIO::POST;
1082 0 : os << " r['] pb[@ 1 GHz] " << LogIO::POST;
1083 0 : Vector<Float> rr;
1084 0 : Vector<Float> pb;
1085 0 : viewPB(rr, pb, nValues);
1086 0 : for (Int ii=0;ii<nValues;ii++) {
1087 0 : os << rr(ii) << " " << pb(ii) << LogIO::POST;
1088 : }
1089 0 : }
1090 : os << "Max Radius at " << refFreq_p.getValue("GHz") << " GHz: "
1091 14 : << maximumRadius_p.getValue("'") << " arcmin " << LogIO::POST;
1092 :
1093 14 : };
1094 :
1095 :
1096 453 : Bool PBMath1D::ok()
1097 : {
1098 453 : if (vp_p.nelements() == 0) {
1099 0 : return false;
1100 453 : } else if (maximumRadius_p.getValue() <= 0.0) {
1101 0 : return false;
1102 453 : } else if (refFreq_p.getValue() <= 0.0) {
1103 0 : return false;
1104 453 : } else if (inverseIncrementRadius_p <= 0.0) {
1105 0 : return false;
1106 : } else {
1107 453 : return true;
1108 : }
1109 : };
1110 :
1111 :
1112 0 : void PBMath1D::viewPB(Vector<Float>& r, Vector<Float>& pb, Int n_pb, const Double freq)
1113 : {
1114 0 : r.resize(n_pb);
1115 0 : pb.resize(n_pb);
1116 0 : if(wideFit_p)
1117 0 : nearestVPArray(freq);
1118 0 : Int nSamples= vp_p.nelements();
1119 0 : for (Int i=0; i< n_pb; i++) {
1120 0 : pb(i) = norm( vp_p( (Int) ((nSamples-1)* (((Float)i)/(n_pb-1) ) ) ) );
1121 0 : r(i) = maximumRadius_p.getValue("'") * (((Float)i)/(n_pb-1));
1122 : }
1123 :
1124 0 : };
1125 0 : void PBMath1D::nearestVPArray(Double freq, bool printINFO){
1126 0 : LogIO os(LogOrigin("PBMATH1D", "nearestVPArray"));
1127 0 : Int ifit=0;
1128 :
1129 0 : Int nFreq=wFreqs_p.nelements();
1130 0 : for (ifit=0; ifit<nFreq; ifit++) {
1131 0 : if (freq <=wFreqs_p(ifit)) break;
1132 : }
1133 0 : if(printINFO)
1134 0 : os << LogIO::NORMAL1 << "Using beam model of frequency " << ((ifit==nFreq) ? wFreqs_p(nFreq-1) : wFreqs_p(ifit)) << " MHz " << LogIO::POST;
1135 0 : if (ifit==0) {
1136 0 : vp_p = wbvp_p.column(0);
1137 0 : } else if (ifit==nFreq) {
1138 0 : vp_p = wbvp_p.column(nFreq-1);
1139 : } else {
1140 0 : Float l = (freq - wFreqs_p(ifit-1))/
1141 0 : (wFreqs_p(ifit)-wFreqs_p(ifit-1));
1142 0 : vp_p = wbvp_p.column(ifit-1)*(1-l) + wbvp_p.column(ifit)*l;
1143 : }
1144 :
1145 :
1146 0 : }
1147 :
1148 : } //# NAMESPACE CASA - END
1149 :
|
