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 4 : PBMath1D::PBMath1D()
67 4 : : composite_p(2048)
68 : {
69 4 : };
70 :
71 :
72 29 : PBMath1D::PBMath1D(Quantity maximumRadius,
73 : Quantity refFreq,
74 : Bool isThisVP,
75 : BeamSquint squint,
76 29 : Bool useSymmetricBeam) :
77 : PBMathInterface(isThisVP, squint, useSymmetricBeam),
78 29 : wideFit_p(false),maximumRadius_p(maximumRadius),
79 29 : refFreq_p(refFreq),
80 58 : composite_p(2048)
81 : {
82 29 : fScale_p = refFreq_p.getValue("GHz"); // scale is ratio of refFreq_p to 1GHz
83 29 : refFreq_p = Quantity( 1.0, "GHz"); // internal Ref Freq is now 1GHz
84 : // convert instantiation parameters to GHz*arcmin reference
85 29 : maximumRadius_p = maximumRadius_p * fScale_p;
86 29 : scale_p = 1.0/(C::arcmin * C::giga);
87 29 : };
88 :
89 33 : PBMath1D::~PBMath1D()
90 : {
91 33 : };
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 0 : Int PBMath1D::support(const CoordinateSystem& cs){
152 0 : Int directionIndex=cs.findCoordinate(Coordinate::DIRECTION);
153 0 : AlwaysAssert(directionIndex>=0, AipsError);
154 : DirectionCoordinate
155 0 : directionCoord=cs.directionCoordinate(directionIndex);
156 :
157 0 : Vector<String> dirunit=directionCoord.worldAxisUnits();
158 :
159 : Double freq;
160 : {
161 0 : Int spectralIndex=cs.findCoordinate(Coordinate::SPECTRAL);
162 0 : AlwaysAssert(spectralIndex>=0, AipsError);
163 : SpectralCoordinate
164 0 : spectralCoord=cs.spectralCoordinate(spectralIndex);
165 :
166 :
167 0 : Vector<String> units(1);
168 0 : units = "Hz";
169 0 : spectralCoord.setWorldAxisUnits(units);
170 :
171 0 : Vector<Double> spectralWorld(1);
172 0 : Vector<Double> spectralPixel(1);
173 0 : spectralPixel(0) = 0;
174 0 : spectralCoord.toWorld(spectralWorld, spectralPixel);
175 0 : freq = spectralWorld(0);
176 0 : }
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 0 : Double numpix=maximumRadius_p.getValue(dirunit(0))/fabs(directionCoord.increment()(0))*2.0*1.0e9/freq ;
186 :
187 :
188 0 : return Int(floor(numpix));
189 :
190 :
191 0 : }
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 16 : 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 32 : LogIO os(LogOrigin("PBMath1D", "apply"));
383 : // Check that in and out are comparable:
384 16 : if (in.shape() != out.shape()) {
385 0 : throw(AipsError("PBMath1D::apply(ImageInterface...) - in and out images have different shapes"));
386 : }
387 16 : CoordinateSystem coords=in.coordinates();
388 16 : if (!coords.near(out.coordinates()) ) {
389 0 : throw(AipsError("PBMath1D::apply(ImageInterface...) - in and out images have different coordinates"));
390 : }
391 :
392 16 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
393 16 : AlwaysAssert(directionIndex>=0, AipsError);
394 : DirectionCoordinate
395 16 : directionCoord=coords.directionCoordinate(directionIndex);
396 16 : Vector<String> units(2); units = "deg";
397 16 : directionCoord.setWorldAxisUnits(units);
398 :
399 : // convert to the EPOCH of these coords
400 16 : MDirection::Types t2 = (MDirection::Types) (pointDir.getRef().getType());
401 16 : MDirection pointDirE;
402 :
403 :
404 16 : 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 16 : pointDirE = pointDir;
412 : }
413 :
414 16 : Int stokesIndex=coords.findCoordinate(Coordinate::STOKES);
415 16 : AlwaysAssert(stokesIndex>=0, AipsError);
416 : StokesCoordinate
417 16 : stokesCoord=coords.stokesCoordinate(stokesIndex);
418 :
419 16 : Int spectralIndex=coords.findCoordinate(Coordinate::SPECTRAL);
420 16 : AlwaysAssert(spectralIndex>=0, AipsError);
421 : SpectralCoordinate
422 16 : spectralCoord=coords.spectralCoordinate(spectralIndex);
423 :
424 16 : units.resize(1);
425 16 : units = "Hz";
426 16 : spectralCoord.setWorldAxisUnits(units);
427 :
428 16 : Int nchan=in.shape()(3);
429 :
430 16 : Vector<Double> pointingCenterWorld(2);
431 16 : Vector<Double> pointingCenterPixel(2);
432 16 : Vector<Double> directionPixel(2);
433 :
434 16 : pointingCenterWorld(0) = pointDirE.getAngle().getValue("deg")(0);
435 16 : pointingCenterWorld(1) = pointDirE.getAngle().getValue("deg")(1);
436 16 : directionCoord.toPixel(pointingCenterPixel, pointingCenterWorld);
437 16 : MDirection newpointDirE;
438 16 : Vector<Double> nonSquintedPointingPixel = pointingCenterPixel.copy();
439 :
440 16 : os << "pointingCenterWorld " << pointingCenterWorld << LogIO::DEBUGGING;
441 16 : os << "pointingCenterPixel " << pointingCenterPixel << LogIO::DEBUGGING;
442 :
443 : // Fill in a cache of the frequencies & squints
444 16 : Vector<Double> spectralWorld(1);
445 16 : Vector<Double> spectralPixel(1);
446 16 : Matrix<Double> xSquintPixCache(2, nchan);
447 16 : Matrix<Double> ySquintPixCache(2, nchan);
448 16 : Vector<Double> spectralCache(nchan);
449 :
450 : {
451 38 : for(Int chan=0;chan<nchan;chan++) {
452 22 : spectralPixel(0)=chan;
453 22 : if(!spectralCoord.toWorld(spectralWorld, spectralPixel)) {
454 0 : os << "Cannot find frequency for this plane" << LogIO::EXCEPTION;
455 : }
456 22 : spectralCache(chan)=spectralWorld(0);
457 :
458 :
459 22 : if (doSquint == BeamSquint::RR || doSquint == BeamSquint::GOFIGURE) {
460 12 : squint_p.getPointingDirection (pointDirE,
461 : parAngle,
462 24 : Quantity(spectralWorld(0),"Hz"),
463 : BeamSquint::RR, newpointDirE);
464 12 : pointingCenterWorld(0) = newpointDirE.getAngle().getValue("deg")(0);
465 12 : pointingCenterWorld(1) = newpointDirE.getAngle().getValue("deg")(1);
466 12 : directionCoord.toPixel(pointingCenterPixel, pointingCenterWorld);
467 12 : xSquintPixCache(0, chan) = pointingCenterPixel(0);
468 12 : ySquintPixCache(0, chan) = pointingCenterPixel(1);
469 : } else {
470 10 : xSquintPixCache(0, chan) = nonSquintedPointingPixel(0);
471 10 : ySquintPixCache(0, chan) = nonSquintedPointingPixel(1);
472 : }
473 22 : if (doSquint == BeamSquint::LL || doSquint == BeamSquint::GOFIGURE) {
474 12 : squint_p.getPointingDirection (pointDirE,
475 : parAngle,
476 24 : Quantity(spectralWorld(0),"Hz"),
477 : BeamSquint::LL, newpointDirE);
478 12 : pointingCenterWorld(0) = newpointDirE.getAngle().getValue("deg")(0);
479 12 : pointingCenterWorld(1) = newpointDirE.getAngle().getValue("deg")(1);
480 12 : directionCoord.toPixel(pointingCenterPixel, pointingCenterWorld);
481 12 : xSquintPixCache(1, chan) = pointingCenterPixel(0);
482 12 : ySquintPixCache(1, chan) = pointingCenterPixel(1);
483 : } else {
484 10 : xSquintPixCache(1, chan) = nonSquintedPointingPixel(0);
485 10 : 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 16 : IPosition ncs=in.shape();
501 16 : ncs(2) = 1; ncs(3) = 1;
502 32 : RO_LatticeIterator<Complex> li(in, LatticeStepper(in.shape(), ncs, IPosition(4,0,1,2,3) ) );
503 32 : 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 16 : Vector<Double> increment = directionCoord.increment();
511 16 : Int rrplane = -1;
512 16 : Int llplane = -1;
513 16 : stokesCoord.toPixel( rrplane, Stokes::RR );
514 16 : 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 16 : Int laststokes = -1;
524 16 : Int lastChan = -1;
525 : Int ichan;
526 : Int istokes;
527 : Int ix0, iy0;
528 : //Int indx;
529 38 : for(li.reset(),oli.reset();!li.atEnd();li++,oli++) {
530 :
531 22 : IPosition itsShape(li.matrixCursor().shape());
532 22 : IPosition loc(li.position());
533 :
534 22 : ichan = loc(3);
535 22 : istokes = loc(2);
536 22 : iy0 = loc(1);
537 22 : 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 22 : if ((doSquint == BeamSquint::RR) ||
545 12 : ((doSquint == BeamSquint::GOFIGURE) && (istokes == rrplane)) ) {
546 0 : xPixel = xSquintPixCache(0, ichan);
547 0 : yPixel = ySquintPixCache(0, ichan);
548 22 : } else if ((doSquint == BeamSquint::LL) ||
549 12 : ((doSquint == BeamSquint::GOFIGURE) && (istokes == llplane)) ) {
550 0 : xPixel = xSquintPixCache(1, ichan);
551 0 : yPixel = ySquintPixCache(1, ichan);
552 : } else {
553 22 : xPixel = nonSquintedPointingPixel(0);
554 22 : yPixel = nonSquintedPointingPixel(1);
555 : }
556 :
557 22 : if (istokes != laststokes) {
558 : //cerr << "Stokes = " << istokes << " pix = " << xPixel << ", " << yPixel << endl;
559 16 : laststokes = istokes;
560 : }
561 :
562 22 : Double factor = 60.0 * spectralCache(ichan)/1.0e+9 ; // arcminutes * GHz
563 22 : Double rmax2 = square( maximumRadius_p.getValue("'") / factor );
564 22 : if (wideFit_p) {
565 : // fill vp with interpolated values for current frequency
566 10 : if (ichan!=lastChan) {
567 10 : nearestVPArray(spectralCache(ichan));
568 10 : lastChan=ichan;
569 : }
570 : }
571 :
572 22 : Vector<Float> rx2(itsShape(0));
573 22 : Vector<Float> ry2(itsShape(1));
574 56622 : for(Int ix=0;ix<itsShape(0);ix++) {
575 56600 : rx2(ix) = square( increment(0)*((Double)(ix+ix0) - xPixel) );
576 : }
577 56622 : for(Int iy=0;iy<itsShape(1);iy++) {
578 56600 : ry2(iy) = square( increment(1)*((Double)(iy+iy0) - yPixel) );
579 : }
580 :
581 22 : const Matrix<Complex>& inmat = li.matrixCursor();
582 22 : Matrix<Complex>& outmat=oli.rwMatrixCursor();
583 :
584 : Bool incopy, outcopy, del;
585 22 : const Complex * inpoint = inmat.getStorage(incopy);
586 22 : Complex *outpoint =outmat.getStorage(outcopy);
587 22 : Float * rx2point = rx2.getStorage(del);
588 22 : Float * ry2point= ry2.getStorage(del);
589 22 : Complex* vppoint=vp_p.getStorage(del);
590 22 : Int nx=itsShape(0);
591 22 : Int ny=itsShape(1);
592 22 : Double inverseIncrementRadius=inverseIncrementRadius_p;
593 22 : #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 22 : outmat.putStorage(outpoint, outcopy);
643 22 : inmat.freeStorage(inpoint, incopy);
644 22 : }
645 :
646 16 : return out;
647 :
648 16 : };
649 :
650 56600 : 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 56600 : Complex taper;
655 203286600 : for(Int ix=0;ix<nx;ix++) {
656 :
657 203230000 : Float r2 = rx2[ix] + ry2;
658 :
659 203230000 : if (r2 > rmax2){
660 199197946 : out[ix+iy*nx] = 0.0;
661 : }
662 : else {
663 4032054 : r = sqrt(r2) * factor;
664 4032054 : indx = Int(r*inverseIncrementRadius);
665 4032054 : if (norm(vp[indx]) > 0.0) {
666 2752806 : if(ipower==2) {
667 2752806 : taper = vp[indx] * conj(vp[indx]);
668 : }
669 : else {
670 0 : taper = vp[indx];
671 : }
672 : } else {
673 1279248 : taper = 0.0;
674 : }
675 4032054 : 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 4032054 : if(!forward) {
682 0 : taper = conj(taper);
683 : }
684 4032054 : 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 4032054 : out[ix+iy*nx] = (in[ix+iy*nx]) * taper ;
692 : }
693 : }
694 : }
695 56600 : };
696 :
697 : ImageInterface<Float>&
698 24 : 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 48 : LogIO os(LogOrigin("PBMath1D", "apply"));
706 :
707 : // cout << "PBMath1D::apply: image shape: " << in.shape() << endl;
708 : // Check that in and out are comparable:
709 24 : if (in.shape() != out.shape()) {
710 0 : throw(AipsError("PBMath1D::apply(ImageInterface...) - in and out images have different shapes"));
711 :
712 : }
713 24 : CoordinateSystem coords=in.coordinates();
714 24 : if (!coords.near(out.coordinates())) {
715 0 : throw(AipsError("PBMath1D::apply(ImageInterface...) - in and out images have different coordinates"));
716 : }
717 :
718 24 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
719 24 : AlwaysAssert(directionIndex>=0, AipsError);
720 : DirectionCoordinate
721 24 : directionCoord=coords.directionCoordinate(directionIndex);
722 24 : Vector<String> units(2); units = "deg";
723 24 : directionCoord.setWorldAxisUnits(units);
724 :
725 : // convert to the EPOCH of these coords
726 24 : MDirection::Types t2 = (MDirection::Types) (pointDir.getRef().getType());
727 24 : MDirection pointDirE(pointDir);
728 :
729 24 : 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 24 : pointDirE = pointDir;
737 : }
738 24 : Int stokesIndex=coords.findCoordinate(Coordinate::STOKES);
739 24 : AlwaysAssert(stokesIndex>=0, AipsError);
740 : StokesCoordinate
741 24 : stokesCoord=coords.stokesCoordinate(stokesIndex);
742 :
743 24 : Int spectralIndex=coords.findCoordinate(Coordinate::SPECTRAL);
744 24 : AlwaysAssert(spectralIndex>=0, AipsError);
745 : SpectralCoordinate
746 24 : spectralCoord=coords.spectralCoordinate(spectralIndex);
747 :
748 24 : units.resize(1);
749 24 : units = "Hz";
750 24 : spectralCoord.setWorldAxisUnits(units);
751 :
752 24 : Int nchan=in.shape()(3);
753 :
754 24 : Vector<Double> pointingCenterWorld(2);
755 24 : Vector<Double> pointingCenterPixel(2);
756 24 : Vector<Double> directionPixel(2);
757 :
758 24 : pointingCenterWorld(0) = pointDirE.getAngle().getValue("deg")(0);
759 24 : pointingCenterWorld(1) = pointDirE.getAngle().getValue("deg")(1);
760 24 : directionCoord.toPixel(pointingCenterPixel, pointingCenterWorld);
761 24 : MDirection newpointDirE;
762 24 : Vector<Double> nonSquintedPointingPixel = pointingCenterPixel.copy();
763 :
764 24 : os << "pointingCenterWorld " << pointingCenterWorld << LogIO::DEBUGGING;
765 24 : os << "pointingCenterPixel " << pointingCenterPixel << LogIO::DEBUGGING;
766 :
767 : // Fill in a cache of the frequencies & squints
768 24 : Vector<Double> spectralWorld(1);
769 24 : Vector<Double> spectralPixel(1);
770 24 : Matrix<Double> xSquintPixCache(2, nchan); // kludge: prevent errors when nchan = 1
771 24 : Matrix<Double> ySquintPixCache(2, nchan);
772 24 : Vector<Double> spectralCache(nchan);
773 :
774 : {
775 60 : for(Int chan=0;chan<nchan;chan++) {
776 36 : spectralPixel(0)=chan;
777 36 : if(!spectralCoord.toWorld(spectralWorld, spectralPixel)) {
778 0 : os << "Cannot find frequency for this plane" << LogIO::EXCEPTION;
779 : }
780 36 : spectralCache(chan)=spectralWorld(0);
781 :
782 :
783 36 : if (doSquint == BeamSquint::RR || doSquint == BeamSquint::GOFIGURE) {
784 36 : squint_p.getPointingDirection (pointDirE,
785 : parAngle,
786 72 : Quantity(spectralWorld(0),"Hz"),
787 : BeamSquint::RR, newpointDirE);
788 36 : pointingCenterWorld(0) = newpointDirE.getAngle().getValue("deg")(0);
789 36 : pointingCenterWorld(1) = newpointDirE.getAngle().getValue("deg")(1);
790 36 : directionCoord.toPixel(pointingCenterPixel, pointingCenterWorld);
791 36 : xSquintPixCache(0, chan) = pointingCenterPixel(0);
792 36 : ySquintPixCache(0, chan) = pointingCenterPixel(1);
793 : } else {
794 0 : xSquintPixCache(0, chan) = nonSquintedPointingPixel(0);
795 0 : ySquintPixCache(0, chan) = nonSquintedPointingPixel(1);
796 : }
797 36 : if (doSquint == BeamSquint::LL || doSquint == BeamSquint::GOFIGURE) {
798 36 : squint_p.getPointingDirection (pointDirE,
799 : parAngle,
800 72 : Quantity(spectralWorld(0),"Hz"),
801 : BeamSquint::LL, newpointDirE);
802 36 : pointingCenterWorld(0) = newpointDirE.getAngle().getValue("deg")(0);
803 36 : pointingCenterWorld(1) = newpointDirE.getAngle().getValue("deg")(1);
804 36 : directionCoord.toPixel(pointingCenterPixel, pointingCenterWorld);
805 36 : xSquintPixCache(1, chan) = pointingCenterPixel(0);
806 36 : 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 24 : IPosition ncs = in.niceCursorShape();
817 : // IPosition ncs=in.shape();
818 24 : ncs(2) = 1; ncs(3) = 1;
819 24 : ncs(2) = 1; ncs(3) = 1;
820 48 : RO_LatticeIterator<Float> li(in, LatticeStepper(in.shape(), ncs, IPosition(4,0,1,2,3) ) );
821 48 : LatticeIterator<Float> oli(out, LatticeStepper(in.shape(), ncs, IPosition(4,0,1,2,3)) );
822 :
823 : //Vector<Float> taper(vp_p.nelements());
824 24 : Float taper=0.0;
825 24 : Float r2=0.0;
826 24 : Float r=0.0;
827 :
828 24 : Vector<Double> increment = directionCoord.increment();
829 24 : Int rrplane = -1;
830 24 : Int llplane = -1;
831 24 : stokesCoord.toPixel( rrplane, Stokes::RR );
832 24 : stokesCoord.toPixel( llplane, Stokes::LL );
833 :
834 : Double xPixel; Double yPixel;
835 :
836 24 : Int laststokes = -1;
837 24 : Int lastChan = -1;
838 : Int ichan;
839 : Int istokes;
840 : Int ix0, iy0;
841 : Int indx;
842 :
843 60 : for(li.reset(),oli.reset();!li.atEnd();li++,oli++) {
844 :
845 36 : IPosition itsShape(li.matrixCursor().shape());
846 36 : IPosition loc(li.position());
847 :
848 36 : ichan = loc(3);
849 36 : istokes = loc(2);
850 36 : iy0 = loc(1);
851 36 : ix0 = loc(0);
852 :
853 : // determine the pointing: RR, LL, or Center?
854 36 : if ((doSquint == BeamSquint::RR) ||
855 36 : ((doSquint == BeamSquint::GOFIGURE) && (istokes == rrplane)) ) {
856 0 : xPixel = xSquintPixCache(0, ichan);
857 0 : yPixel = ySquintPixCache(0, ichan);
858 36 : } else if ((doSquint == BeamSquint::LL) ||
859 36 : ((doSquint == BeamSquint::GOFIGURE) && (istokes == llplane ))) {
860 0 : xPixel = xSquintPixCache(1, ichan);
861 0 : yPixel = ySquintPixCache(1, ichan);
862 : } else {
863 36 : xPixel = nonSquintedPointingPixel(0);
864 36 : yPixel = nonSquintedPointingPixel(1);
865 : }
866 :
867 36 : if (istokes != laststokes) {
868 : // cout << "Stokes = " << istokes << " pix = " << xPixel << ", " << yPixel << endl;
869 24 : laststokes = istokes;
870 : }
871 :
872 36 : Double factor = 60.0 * spectralCache(ichan)/1.0e+9 ; // arcminutes * GHz
873 36 : Double rmax2 = square( maximumRadius_p.getValue("'") / factor );
874 36 : 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 36 : Vector<Float> rx2(itsShape(0));
894 36 : Vector<Float> ry2(itsShape(1));
895 49116 : for(Int ix=0;ix<itsShape(0);ix++) {
896 49080 : rx2(ix) = square( increment(0)*((Double)(ix+ix0) - xPixel) );
897 : }
898 49116 : for(Int iy=0;iy<itsShape(1);iy++) {
899 49080 : ry2(iy) = square( increment(1)*((Double)(iy+iy0) - yPixel) );
900 : }
901 :
902 49116 : for(Int iy=0;iy<itsShape(1);iy++) {
903 196888680 : for(Int ix=0;ix<itsShape(0);ix++) {
904 :
905 196839600 : r2 = rx2(ix) + ry2(iy);
906 196839600 : if (r2 > rmax2) {
907 195395076 : oli.rwMatrixCursor()(ix, iy) = 0.0;
908 : } else {
909 1444524 : r = sqrt(r2) * factor;
910 1444524 : indx = Int(r*inverseIncrementRadius_p);
911 1444524 : if (norm(vp_p(indx)) > 0.0) {
912 1444524 : taper = real(vp_p(indx)) * real(vp_p(indx))+ imag(vp_p(indx)) * imag(vp_p(indx));
913 1444524 : if(ipower==4)
914 1444524 : taper *= taper;
915 : }
916 : //else {
917 : // taper = 0.0;
918 : //}
919 1444524 : oli.rwMatrixCursor()(ix, iy) = li.matrixCursor()(ix,iy) * taper;
920 : }
921 : }
922 : }
923 36 : }
924 24 : return out;
925 :
926 24 : };
927 :
928 : // Behavior: doSquint == RR or LL don't make sense here
929 : //
930 : //
931 :
932 : SkyComponent&
933 0 : 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 0 : MDirection::Types t1 = (MDirection::Types) (in.shape().refDirection().getRef().getType());
955 0 : MDirection::Types t2 = (MDirection::Types) (pointDir.getRef().getType());
956 :
957 0 : MDirection pointDirE;
958 0 : if ( t1 != t2) {
959 0 : MDirection::Convert converter;
960 0 : converter.setOut( t1 );
961 0 : pointDirE = converter(pointDir);
962 0 : } else {
963 0 : pointDirE = pointDir;
964 : }
965 :
966 0 : 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 0 : if (doSquint == BeamSquint::GOFIGURE) {
970 0 : if (in.flux().pol() != ComponentType::CIRCULAR) {
971 0 : in.flux().convertPol(ComponentType::CIRCULAR);
972 : }
973 : }
974 :
975 0 : Vector<DComplex> compFluxIn = in.flux().value();
976 0 : Vector<DComplex> compFlux = out.flux().value();
977 0 : compFlux = compFluxIn.copy();
978 :
979 : // Find the direction of the component
980 0 : 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 0 : 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 0 : MDirection newpointDirE;
1001 0 : for (Int pol=0;pol<4;pol++) {
1002 0 : Stokes::StokesTypes stokes=Stokes::type(pol+5);
1003 :
1004 0 : if (stokes == Stokes::RR && doSquint == BeamSquint::GOFIGURE) {
1005 0 : squint_p.getPointingDirection (pointDirE, parAngle, frequency, BeamSquint::RR,
1006 : newpointDirE );
1007 0 : } else if (stokes == Stokes::LL && doSquint == BeamSquint::GOFIGURE) {
1008 0 : squint_p.getPointingDirection (pointDirE, parAngle, frequency, BeamSquint::LL,
1009 : newpointDirE );
1010 : } else {
1011 0 : newpointDirE = pointDirE;
1012 : }
1013 :
1014 0 : MVDirection mvd1( compDir.getAngle() );
1015 0 : MVDirection mvd2( newpointDirE.getAngle() );
1016 0 : Quantity sep = mvd1.separation(mvd2, "'");
1017 0 : double r = sep.getValue("'") * frequency.getValue("Hz") / 1.0e+9; // arcminutes * GHz
1018 0 : Complex taper;
1019 0 : Int ir = Int(r*inverseIncrementRadius_p);
1020 : //vp_p is interpolated wvp_p from above
1021 0 : 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 0 : if (r > maximumRadius_p.getValue("'")) {
1034 0 : compFlux(pol) = 0.0;
1035 : } else {
1036 0 : if (norm(vpVal) > 0.0) {
1037 0 : if(iPower>1){
1038 0 : taper=vpVal*conj(vpVal);
1039 0 : 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 0 : if (conjugate) {
1050 0 : taper = conj(taper);
1051 : }
1052 0 : 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 0 : compFlux(pol) *= taper;
1060 : }
1061 : }
1062 0 : }
1063 :
1064 : // Set the output component fluxes
1065 0 : out = in.copy();
1066 0 : out.flux().setValue(compFlux);
1067 :
1068 0 : return out;
1069 :
1070 0 : };
1071 :
1072 0 : void PBMath1D::summary(Int nValues)
1073 : {
1074 0 : String name;
1075 0 : namePBClass(name);
1076 0 : LogIO os(LogOrigin("PBMath1D", "summary"));
1077 0 : os << "Using " << name << " PB Class " << LogIO::POST;
1078 0 : PBMathInterface::summary(nValues);
1079 :
1080 0 : 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 0 : << maximumRadius_p.getValue("'") << " arcmin " << LogIO::POST;
1092 :
1093 0 : };
1094 :
1095 :
1096 36 : Bool PBMath1D::ok()
1097 : {
1098 36 : if (vp_p.nelements() == 0) {
1099 0 : return false;
1100 36 : } else if (maximumRadius_p.getValue() <= 0.0) {
1101 0 : return false;
1102 36 : } else if (refFreq_p.getValue() <= 0.0) {
1103 0 : return false;
1104 36 : } else if (inverseIncrementRadius_p <= 0.0) {
1105 0 : return false;
1106 : } else {
1107 36 : 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 :
|
