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 286 : PBMath1D::PBMath1D()
67 286 : : composite_p(2048)
68 : {
69 286 : };
70 :
71 :
72 2389 : PBMath1D::PBMath1D(Quantity maximumRadius,
73 : Quantity refFreq,
74 : Bool isThisVP,
75 : BeamSquint squint,
76 2389 : Bool useSymmetricBeam) :
77 : PBMathInterface(isThisVP, squint, useSymmetricBeam),
78 2389 : wideFit_p(false),maximumRadius_p(maximumRadius),
79 2389 : refFreq_p(refFreq),
80 4778 : composite_p(2048)
81 : {
82 2389 : fScale_p = refFreq_p.getValue("GHz"); // scale is ratio of refFreq_p to 1GHz
83 2389 : refFreq_p = Quantity( 1.0, "GHz"); // internal Ref Freq is now 1GHz
84 : // convert instantiation parameters to GHz*arcmin reference
85 2389 : maximumRadius_p = maximumRadius_p * fScale_p;
86 2389 : scale_p = 1.0/(C::arcmin * C::giga);
87 2389 : };
88 :
89 2653 : PBMath1D::~PBMath1D()
90 : {
91 2653 : };
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 377 : Int PBMath1D::support(const CoordinateSystem& cs){
152 377 : Int directionIndex=cs.findCoordinate(Coordinate::DIRECTION);
153 377 : AlwaysAssert(directionIndex>=0, AipsError);
154 : DirectionCoordinate
155 377 : directionCoord=cs.directionCoordinate(directionIndex);
156 :
157 377 : Vector<String> dirunit=directionCoord.worldAxisUnits();
158 :
159 : Double freq;
160 : {
161 377 : Int spectralIndex=cs.findCoordinate(Coordinate::SPECTRAL);
162 377 : AlwaysAssert(spectralIndex>=0, AipsError);
163 : SpectralCoordinate
164 377 : spectralCoord=cs.spectralCoordinate(spectralIndex);
165 :
166 :
167 377 : Vector<String> units(1);
168 377 : units = "Hz";
169 377 : spectralCoord.setWorldAxisUnits(units);
170 :
171 377 : Vector<Double> spectralWorld(1);
172 377 : Vector<Double> spectralPixel(1);
173 377 : spectralPixel(0) = 0;
174 377 : spectralCoord.toWorld(spectralWorld, spectralPixel);
175 377 : freq = spectralWorld(0);
176 377 : }
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 377 : Double numpix=maximumRadius_p.getValue(dirunit(0))/fabs(directionCoord.increment()(0))*2.0*1.0e9/freq ;
186 :
187 :
188 377 : return Int(floor(numpix));
189 :
190 :
191 377 : }
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 1174 : 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 1174 : if (vp_p.nelements() == 0) {
365 0 : fillPBArray();
366 : }
367 1174 : esvp_p = vp_p;
368 1174 : };
369 :
370 : ImageInterface<Complex>&
371 3362 : 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 6724 : LogIO os(LogOrigin("PBMath1D", "apply"));
383 : // Check that in and out are comparable:
384 3362 : if (in.shape() != out.shape()) {
385 0 : throw(AipsError("PBMath1D::apply(ImageInterface...) - in and out images have different shapes"));
386 : }
387 3362 : CoordinateSystem coords=in.coordinates();
388 3362 : if (!coords.near(out.coordinates()) ) {
389 0 : throw(AipsError("PBMath1D::apply(ImageInterface...) - in and out images have different coordinates"));
390 : }
391 :
392 3362 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
393 3362 : AlwaysAssert(directionIndex>=0, AipsError);
394 : DirectionCoordinate
395 3362 : directionCoord=coords.directionCoordinate(directionIndex);
396 3362 : Vector<String> units(2); units = "deg";
397 3362 : directionCoord.setWorldAxisUnits(units);
398 :
399 : // convert to the EPOCH of these coords
400 3362 : MDirection::Types t2 = (MDirection::Types) (pointDir.getRef().getType());
401 3362 : MDirection pointDirE;
402 :
403 :
404 3362 : 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 3362 : pointDirE = pointDir;
412 : }
413 :
414 3362 : Int stokesIndex=coords.findCoordinate(Coordinate::STOKES);
415 3362 : AlwaysAssert(stokesIndex>=0, AipsError);
416 : StokesCoordinate
417 3362 : stokesCoord=coords.stokesCoordinate(stokesIndex);
418 :
419 3362 : Int spectralIndex=coords.findCoordinate(Coordinate::SPECTRAL);
420 3362 : AlwaysAssert(spectralIndex>=0, AipsError);
421 : SpectralCoordinate
422 3362 : spectralCoord=coords.spectralCoordinate(spectralIndex);
423 :
424 3362 : units.resize(1);
425 3362 : units = "Hz";
426 3362 : spectralCoord.setWorldAxisUnits(units);
427 :
428 3362 : Int nchan=in.shape()(3);
429 :
430 3362 : Vector<Double> pointingCenterWorld(2);
431 3362 : Vector<Double> pointingCenterPixel(2);
432 3362 : Vector<Double> directionPixel(2);
433 :
434 3362 : pointingCenterWorld(0) = pointDirE.getAngle().getValue("deg")(0);
435 3362 : pointingCenterWorld(1) = pointDirE.getAngle().getValue("deg")(1);
436 3362 : directionCoord.toPixel(pointingCenterPixel, pointingCenterWorld);
437 3362 : MDirection newpointDirE;
438 3362 : Vector<Double> nonSquintedPointingPixel = pointingCenterPixel.copy();
439 :
440 3362 : os << "pointingCenterWorld " << pointingCenterWorld << LogIO::DEBUGGING;
441 3362 : os << "pointingCenterPixel " << pointingCenterPixel << LogIO::DEBUGGING;
442 :
443 : // Fill in a cache of the frequencies & squints
444 3362 : Vector<Double> spectralWorld(1);
445 3362 : Vector<Double> spectralPixel(1);
446 3362 : Matrix<Double> xSquintPixCache(2, nchan);
447 3362 : Matrix<Double> ySquintPixCache(2, nchan);
448 3362 : Vector<Double> spectralCache(nchan);
449 :
450 : {
451 8680 : for(Int chan=0;chan<nchan;chan++) {
452 5318 : spectralPixel(0)=chan;
453 5318 : if(!spectralCoord.toWorld(spectralWorld, spectralPixel)) {
454 0 : os << "Cannot find frequency for this plane" << LogIO::EXCEPTION;
455 : }
456 5318 : spectralCache(chan)=spectralWorld(0);
457 :
458 :
459 5318 : if (doSquint == BeamSquint::RR || doSquint == BeamSquint::GOFIGURE) {
460 77 : squint_p.getPointingDirection (pointDirE,
461 : parAngle,
462 154 : Quantity(spectralWorld(0),"Hz"),
463 : BeamSquint::RR, newpointDirE);
464 77 : pointingCenterWorld(0) = newpointDirE.getAngle().getValue("deg")(0);
465 77 : pointingCenterWorld(1) = newpointDirE.getAngle().getValue("deg")(1);
466 77 : directionCoord.toPixel(pointingCenterPixel, pointingCenterWorld);
467 77 : xSquintPixCache(0, chan) = pointingCenterPixel(0);
468 77 : ySquintPixCache(0, chan) = pointingCenterPixel(1);
469 : } else {
470 5241 : xSquintPixCache(0, chan) = nonSquintedPointingPixel(0);
471 5241 : ySquintPixCache(0, chan) = nonSquintedPointingPixel(1);
472 : }
473 5318 : if (doSquint == BeamSquint::LL || doSquint == BeamSquint::GOFIGURE) {
474 63 : squint_p.getPointingDirection (pointDirE,
475 : parAngle,
476 126 : Quantity(spectralWorld(0),"Hz"),
477 : BeamSquint::LL, newpointDirE);
478 63 : pointingCenterWorld(0) = newpointDirE.getAngle().getValue("deg")(0);
479 63 : pointingCenterWorld(1) = newpointDirE.getAngle().getValue("deg")(1);
480 63 : directionCoord.toPixel(pointingCenterPixel, pointingCenterWorld);
481 63 : xSquintPixCache(1, chan) = pointingCenterPixel(0);
482 63 : ySquintPixCache(1, chan) = pointingCenterPixel(1);
483 : } else {
484 5255 : xSquintPixCache(1, chan) = nonSquintedPointingPixel(0);
485 5255 : 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 3362 : IPosition ncs=in.shape();
501 3362 : ncs(2) = 1; ncs(3) = 1;
502 6724 : RO_LatticeIterator<Complex> li(in, LatticeStepper(in.shape(), ncs, IPosition(4,0,1,2,3) ) );
503 6724 : 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 3362 : Vector<Double> increment = directionCoord.increment();
511 3362 : Int rrplane = -1;
512 3362 : Int llplane = -1;
513 3362 : stokesCoord.toPixel( rrplane, Stokes::RR );
514 3362 : 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 3362 : Int laststokes = -1;
524 3362 : Int lastChan = -1;
525 : Int ichan;
526 : Int istokes;
527 : Int ix0, iy0;
528 : //Int indx;
529 8913 : for(li.reset(),oli.reset();!li.atEnd();li++,oli++) {
530 :
531 5551 : IPosition itsShape(li.matrixCursor().shape());
532 5551 : IPosition loc(li.position());
533 :
534 5551 : ichan = loc(3);
535 5551 : istokes = loc(2);
536 5551 : iy0 = loc(1);
537 5551 : 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 5551 : if ((doSquint == BeamSquint::RR) ||
545 63 : ((doSquint == BeamSquint::GOFIGURE) && (istokes == rrplane)) ) {
546 14 : xPixel = xSquintPixCache(0, ichan);
547 14 : yPixel = ySquintPixCache(0, ichan);
548 5537 : } else if ((doSquint == BeamSquint::LL) ||
549 63 : ((doSquint == BeamSquint::GOFIGURE) && (istokes == llplane)) ) {
550 0 : xPixel = xSquintPixCache(1, ichan);
551 0 : yPixel = ySquintPixCache(1, ichan);
552 : } else {
553 5537 : xPixel = nonSquintedPointingPixel(0);
554 5537 : yPixel = nonSquintedPointingPixel(1);
555 : }
556 :
557 5551 : if (istokes != laststokes) {
558 : //cerr << "Stokes = " << istokes << " pix = " << xPixel << ", " << yPixel << endl;
559 3635 : laststokes = istokes;
560 : }
561 :
562 5551 : Double factor = 60.0 * spectralCache(ichan)/1.0e+9 ; // arcminutes * GHz
563 5551 : Double rmax2 = square( maximumRadius_p.getValue("'") / factor );
564 5551 : if (wideFit_p) {
565 : // fill vp with interpolated values for current frequency
566 1340 : if (ichan!=lastChan) {
567 1340 : nearestVPArray(spectralCache(ichan));
568 1340 : lastChan=ichan;
569 : }
570 : }
571 :
572 5551 : Vector<Float> rx2(itsShape(0));
573 5551 : Vector<Float> ry2(itsShape(1));
574 3892346 : for(Int ix=0;ix<itsShape(0);ix++) {
575 3886795 : rx2(ix) = square( increment(0)*((Double)(ix+ix0) - xPixel) );
576 : }
577 2772064 : for(Int iy=0;iy<itsShape(1);iy++) {
578 2766513 : ry2(iy) = square( increment(1)*((Double)(iy+iy0) - yPixel) );
579 : }
580 :
581 5551 : const Matrix<Complex>& inmat = li.matrixCursor();
582 5551 : Matrix<Complex>& outmat=oli.rwMatrixCursor();
583 :
584 : Bool incopy, outcopy, del;
585 5551 : const Complex * inpoint = inmat.getStorage(incopy);
586 5551 : Complex *outpoint =outmat.getStorage(outcopy);
587 5551 : Float * rx2point = rx2.getStorage(del);
588 5551 : Float * ry2point= ry2.getStorage(del);
589 5551 : Complex* vppoint=vp_p.getStorage(del);
590 5551 : Int nx=itsShape(0);
591 5551 : Int ny=itsShape(1);
592 5551 : Double inverseIncrementRadius=inverseIncrementRadius_p;
593 5551 : #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 5551 : outmat.putStorage(outpoint, outcopy);
643 5551 : inmat.freeStorage(inpoint, incopy);
644 5551 : }
645 :
646 3362 : return out;
647 :
648 3362 : };
649 :
650 2766513 : 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 2766513 : Complex taper;
655 3207568346 : for(Int ix=0;ix<nx;ix++) {
656 :
657 3204801833 : Float r2 = rx2[ix] + ry2;
658 :
659 3204801833 : if (r2 > rmax2){
660 2053817738 : out[ix+iy*nx] = 0.0;
661 : }
662 : else {
663 1150984095 : r = sqrt(r2) * factor;
664 1150984095 : indx = Int(r*inverseIncrementRadius);
665 1150984095 : if (norm(vp[indx]) > 0.0) {
666 1092646399 : if(ipower==2) {
667 597319256 : taper = vp[indx] * conj(vp[indx]);
668 : }
669 : else {
670 495327143 : taper = vp[indx];
671 : }
672 : } else {
673 58337696 : taper = 0.0;
674 : }
675 1150984095 : 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 1150984095 : if(!forward) {
682 0 : taper = conj(taper);
683 : }
684 1150984095 : 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 1150984095 : out[ix+iy*nx] = (in[ix+iy*nx]) * taper ;
692 : }
693 : }
694 : }
695 2766513 : };
696 :
697 : ImageInterface<Float>&
698 40 : 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 80 : LogIO os(LogOrigin("PBMath1D", "apply"));
706 :
707 : // cout << "PBMath1D::apply: image shape: " << in.shape() << endl;
708 : // Check that in and out are comparable:
709 40 : if (in.shape() != out.shape()) {
710 0 : throw(AipsError("PBMath1D::apply(ImageInterface...) - in and out images have different shapes"));
711 :
712 : }
713 40 : CoordinateSystem coords=in.coordinates();
714 40 : if (!coords.near(out.coordinates())) {
715 0 : throw(AipsError("PBMath1D::apply(ImageInterface...) - in and out images have different coordinates"));
716 : }
717 :
718 40 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
719 40 : AlwaysAssert(directionIndex>=0, AipsError);
720 : DirectionCoordinate
721 40 : directionCoord=coords.directionCoordinate(directionIndex);
722 40 : Vector<String> units(2); units = "deg";
723 40 : directionCoord.setWorldAxisUnits(units);
724 :
725 : // convert to the EPOCH of these coords
726 40 : MDirection::Types t2 = (MDirection::Types) (pointDir.getRef().getType());
727 40 : MDirection pointDirE(pointDir);
728 :
729 40 : 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 40 : pointDirE = pointDir;
737 : }
738 40 : Int stokesIndex=coords.findCoordinate(Coordinate::STOKES);
739 40 : AlwaysAssert(stokesIndex>=0, AipsError);
740 : StokesCoordinate
741 40 : stokesCoord=coords.stokesCoordinate(stokesIndex);
742 :
743 40 : Int spectralIndex=coords.findCoordinate(Coordinate::SPECTRAL);
744 40 : AlwaysAssert(spectralIndex>=0, AipsError);
745 : SpectralCoordinate
746 40 : spectralCoord=coords.spectralCoordinate(spectralIndex);
747 :
748 40 : units.resize(1);
749 40 : units = "Hz";
750 40 : spectralCoord.setWorldAxisUnits(units);
751 :
752 40 : Int nchan=in.shape()(3);
753 :
754 40 : Vector<Double> pointingCenterWorld(2);
755 40 : Vector<Double> pointingCenterPixel(2);
756 40 : Vector<Double> directionPixel(2);
757 :
758 40 : pointingCenterWorld(0) = pointDirE.getAngle().getValue("deg")(0);
759 40 : pointingCenterWorld(1) = pointDirE.getAngle().getValue("deg")(1);
760 40 : directionCoord.toPixel(pointingCenterPixel, pointingCenterWorld);
761 40 : MDirection newpointDirE;
762 40 : Vector<Double> nonSquintedPointingPixel = pointingCenterPixel.copy();
763 :
764 40 : os << "pointingCenterWorld " << pointingCenterWorld << LogIO::DEBUGGING;
765 40 : os << "pointingCenterPixel " << pointingCenterPixel << LogIO::DEBUGGING;
766 :
767 : // Fill in a cache of the frequencies & squints
768 40 : Vector<Double> spectralWorld(1);
769 40 : Vector<Double> spectralPixel(1);
770 40 : Matrix<Double> xSquintPixCache(2, nchan); // kludge: prevent errors when nchan = 1
771 40 : Matrix<Double> ySquintPixCache(2, nchan);
772 40 : Vector<Double> spectralCache(nchan);
773 :
774 : {
775 108 : for(Int chan=0;chan<nchan;chan++) {
776 68 : spectralPixel(0)=chan;
777 68 : if(!spectralCoord.toWorld(spectralWorld, spectralPixel)) {
778 0 : os << "Cannot find frequency for this plane" << LogIO::EXCEPTION;
779 : }
780 68 : spectralCache(chan)=spectralWorld(0);
781 :
782 :
783 68 : if (doSquint == BeamSquint::RR || doSquint == BeamSquint::GOFIGURE) {
784 68 : squint_p.getPointingDirection (pointDirE,
785 : parAngle,
786 136 : Quantity(spectralWorld(0),"Hz"),
787 : BeamSquint::RR, newpointDirE);
788 68 : pointingCenterWorld(0) = newpointDirE.getAngle().getValue("deg")(0);
789 68 : pointingCenterWorld(1) = newpointDirE.getAngle().getValue("deg")(1);
790 68 : directionCoord.toPixel(pointingCenterPixel, pointingCenterWorld);
791 68 : xSquintPixCache(0, chan) = pointingCenterPixel(0);
792 68 : ySquintPixCache(0, chan) = pointingCenterPixel(1);
793 : } else {
794 0 : xSquintPixCache(0, chan) = nonSquintedPointingPixel(0);
795 0 : ySquintPixCache(0, chan) = nonSquintedPointingPixel(1);
796 : }
797 68 : if (doSquint == BeamSquint::LL || doSquint == BeamSquint::GOFIGURE) {
798 68 : squint_p.getPointingDirection (pointDirE,
799 : parAngle,
800 136 : Quantity(spectralWorld(0),"Hz"),
801 : BeamSquint::LL, newpointDirE);
802 68 : pointingCenterWorld(0) = newpointDirE.getAngle().getValue("deg")(0);
803 68 : pointingCenterWorld(1) = newpointDirE.getAngle().getValue("deg")(1);
804 68 : directionCoord.toPixel(pointingCenterPixel, pointingCenterWorld);
805 68 : xSquintPixCache(1, chan) = pointingCenterPixel(0);
806 68 : 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 40 : IPosition ncs = in.niceCursorShape();
817 : // IPosition ncs=in.shape();
818 40 : ncs(2) = 1; ncs(3) = 1;
819 40 : ncs(2) = 1; ncs(3) = 1;
820 80 : RO_LatticeIterator<Float> li(in, LatticeStepper(in.shape(), ncs, IPosition(4,0,1,2,3) ) );
821 80 : LatticeIterator<Float> oli(out, LatticeStepper(in.shape(), ncs, IPosition(4,0,1,2,3)) );
822 :
823 : //Vector<Float> taper(vp_p.nelements());
824 40 : Float taper=0.0;
825 40 : Float r2=0.0;
826 40 : Float r=0.0;
827 :
828 40 : Vector<Double> increment = directionCoord.increment();
829 40 : Int rrplane = -1;
830 40 : Int llplane = -1;
831 40 : stokesCoord.toPixel( rrplane, Stokes::RR );
832 40 : stokesCoord.toPixel( llplane, Stokes::LL );
833 :
834 : Double xPixel; Double yPixel;
835 :
836 40 : Int laststokes = -1;
837 40 : Int lastChan = -1;
838 : Int ichan;
839 : Int istokes;
840 : Int ix0, iy0;
841 : Int indx;
842 :
843 108 : for(li.reset(),oli.reset();!li.atEnd();li++,oli++) {
844 :
845 68 : IPosition itsShape(li.matrixCursor().shape());
846 68 : IPosition loc(li.position());
847 :
848 68 : ichan = loc(3);
849 68 : istokes = loc(2);
850 68 : iy0 = loc(1);
851 68 : ix0 = loc(0);
852 :
853 : // determine the pointing: RR, LL, or Center?
854 68 : if ((doSquint == BeamSquint::RR) ||
855 68 : ((doSquint == BeamSquint::GOFIGURE) && (istokes == rrplane)) ) {
856 0 : xPixel = xSquintPixCache(0, ichan);
857 0 : yPixel = ySquintPixCache(0, ichan);
858 68 : } else if ((doSquint == BeamSquint::LL) ||
859 68 : ((doSquint == BeamSquint::GOFIGURE) && (istokes == llplane ))) {
860 0 : xPixel = xSquintPixCache(1, ichan);
861 0 : yPixel = ySquintPixCache(1, ichan);
862 : } else {
863 68 : xPixel = nonSquintedPointingPixel(0);
864 68 : yPixel = nonSquintedPointingPixel(1);
865 : }
866 :
867 68 : if (istokes != laststokes) {
868 : // cout << "Stokes = " << istokes << " pix = " << xPixel << ", " << yPixel << endl;
869 40 : laststokes = istokes;
870 : }
871 :
872 68 : Double factor = 60.0 * spectralCache(ichan)/1.0e+9 ; // arcminutes * GHz
873 68 : Double rmax2 = square( maximumRadius_p.getValue("'") / factor );
874 68 : 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 68 : Vector<Float> rx2(itsShape(0));
894 68 : Vector<Float> ry2(itsShape(1));
895 135868 : for(Int ix=0;ix<itsShape(0);ix++) {
896 135800 : rx2(ix) = square( increment(0)*((Double)(ix+ix0) - xPixel) );
897 : }
898 135868 : for(Int iy=0;iy<itsShape(1);iy++) {
899 135800 : ry2(iy) = square( increment(1)*((Double)(iy+iy0) - yPixel) );
900 : }
901 :
902 135868 : for(Int iy=0;iy<itsShape(1);iy++) {
903 1036796200 : for(Int ix=0;ix<itsShape(0);ix++) {
904 :
905 1036660400 : r2 = rx2(ix) + ry2(iy);
906 1036660400 : if (r2 > rmax2) {
907 1032186012 : oli.rwMatrixCursor()(ix, iy) = 0.0;
908 : } else {
909 4474388 : r = sqrt(r2) * factor;
910 4474388 : indx = Int(r*inverseIncrementRadius_p);
911 4474388 : if (norm(vp_p(indx)) > 0.0) {
912 4474388 : taper = real(vp_p(indx)) * real(vp_p(indx))+ imag(vp_p(indx)) * imag(vp_p(indx));
913 4474388 : if(ipower==4)
914 4474388 : taper *= taper;
915 : }
916 : //else {
917 : // taper = 0.0;
918 : //}
919 4474388 : oli.rwMatrixCursor()(ix, iy) = li.matrixCursor()(ix,iy) * taper;
920 : }
921 : }
922 : }
923 68 : }
924 40 : return out;
925 :
926 40 : };
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 17 : void PBMath1D::summary(Int nValues)
1073 : {
1074 17 : String name;
1075 17 : namePBClass(name);
1076 34 : LogIO os(LogOrigin("PBMath1D", "summary"));
1077 17 : os << "Using " << name << " PB Class " << LogIO::POST;
1078 17 : PBMathInterface::summary(nValues);
1079 :
1080 17 : 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 17 : << maximumRadius_p.getValue("'") << " arcmin " << LogIO::POST;
1092 :
1093 17 : };
1094 :
1095 :
1096 636 : Bool PBMath1D::ok()
1097 : {
1098 636 : if (vp_p.nelements() == 0) {
1099 0 : return false;
1100 636 : } else if (maximumRadius_p.getValue() <= 0.0) {
1101 0 : return false;
1102 636 : } else if (refFreq_p.getValue() <= 0.0) {
1103 0 : return false;
1104 636 : } else if (inverseIncrementRadius_p <= 0.0) {
1105 0 : return false;
1106 : } else {
1107 636 : 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 :
|
