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 73 : PBMath1D::PBMath1D(Quantity maximumRadius,
73 : Quantity refFreq,
74 : Bool isThisVP,
75 : BeamSquint squint,
76 73 : Bool useSymmetricBeam) :
77 : PBMathInterface(isThisVP, squint, useSymmetricBeam),
78 73 : wideFit_p(false),maximumRadius_p(maximumRadius),
79 73 : refFreq_p(refFreq),
80 146 : composite_p(2048)
81 : {
82 73 : fScale_p = refFreq_p.getValue("GHz"); // scale is ratio of refFreq_p to 1GHz
83 73 : refFreq_p = Quantity( 1.0, "GHz"); // internal Ref Freq is now 1GHz
84 : // convert instantiation parameters to GHz*arcmin reference
85 73 : maximumRadius_p = maximumRadius_p * fScale_p;
86 73 : scale_p = 1.0/(C::arcmin * C::giga);
87 73 : };
88 :
89 73 : PBMath1D::~PBMath1D()
90 : {
91 73 : };
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 74 : 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 148 : LogIO os(LogOrigin("PBMath1D", "apply"));
383 : // Check that in and out are comparable:
384 74 : if (in.shape() != out.shape()) {
385 0 : throw(AipsError("PBMath1D::apply(ImageInterface...) - in and out images have different shapes"));
386 : }
387 74 : CoordinateSystem coords=in.coordinates();
388 74 : if (!coords.near(out.coordinates()) ) {
389 0 : throw(AipsError("PBMath1D::apply(ImageInterface...) - in and out images have different coordinates"));
390 : }
391 :
392 74 : Int directionIndex=coords.findCoordinate(Coordinate::DIRECTION);
393 74 : AlwaysAssert(directionIndex>=0, AipsError);
394 : DirectionCoordinate
395 74 : directionCoord=coords.directionCoordinate(directionIndex);
396 74 : Vector<String> units(2); units = "deg";
397 74 : directionCoord.setWorldAxisUnits(units);
398 :
399 : // convert to the EPOCH of these coords
400 74 : MDirection::Types t2 = (MDirection::Types) (pointDir.getRef().getType());
401 74 : MDirection pointDirE;
402 :
403 :
404 74 : 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 74 : pointDirE = pointDir;
412 : }
413 :
414 74 : Int stokesIndex=coords.findCoordinate(Coordinate::STOKES);
415 74 : AlwaysAssert(stokesIndex>=0, AipsError);
416 : StokesCoordinate
417 74 : stokesCoord=coords.stokesCoordinate(stokesIndex);
418 :
419 74 : Int spectralIndex=coords.findCoordinate(Coordinate::SPECTRAL);
420 74 : AlwaysAssert(spectralIndex>=0, AipsError);
421 : SpectralCoordinate
422 74 : spectralCoord=coords.spectralCoordinate(spectralIndex);
423 :
424 74 : units.resize(1);
425 74 : units = "Hz";
426 74 : spectralCoord.setWorldAxisUnits(units);
427 :
428 74 : Int nchan=in.shape()(3);
429 :
430 74 : Vector<Double> pointingCenterWorld(2);
431 74 : Vector<Double> pointingCenterPixel(2);
432 74 : Vector<Double> directionPixel(2);
433 :
434 74 : pointingCenterWorld(0) = pointDirE.getAngle().getValue("deg")(0);
435 74 : pointingCenterWorld(1) = pointDirE.getAngle().getValue("deg")(1);
436 74 : directionCoord.toPixel(pointingCenterPixel, pointingCenterWorld);
437 74 : MDirection newpointDirE;
438 74 : Vector<Double> nonSquintedPointingPixel = pointingCenterPixel.copy();
439 :
440 74 : os << "pointingCenterWorld " << pointingCenterWorld << LogIO::DEBUGGING;
441 74 : os << "pointingCenterPixel " << pointingCenterPixel << LogIO::DEBUGGING;
442 :
443 : // Fill in a cache of the frequencies & squints
444 74 : Vector<Double> spectralWorld(1);
445 74 : Vector<Double> spectralPixel(1);
446 74 : Matrix<Double> xSquintPixCache(2, nchan);
447 74 : Matrix<Double> ySquintPixCache(2, nchan);
448 74 : Vector<Double> spectralCache(nchan);
449 :
450 : {
451 148 : for(Int chan=0;chan<nchan;chan++) {
452 74 : spectralPixel(0)=chan;
453 74 : if(!spectralCoord.toWorld(spectralWorld, spectralPixel)) {
454 0 : os << "Cannot find frequency for this plane" << LogIO::EXCEPTION;
455 : }
456 74 : spectralCache(chan)=spectralWorld(0);
457 :
458 :
459 74 : if (doSquint == BeamSquint::RR || doSquint == BeamSquint::GOFIGURE) {
460 0 : squint_p.getPointingDirection (pointDirE,
461 : parAngle,
462 0 : Quantity(spectralWorld(0),"Hz"),
463 : BeamSquint::RR, newpointDirE);
464 0 : pointingCenterWorld(0) = newpointDirE.getAngle().getValue("deg")(0);
465 0 : pointingCenterWorld(1) = newpointDirE.getAngle().getValue("deg")(1);
466 0 : directionCoord.toPixel(pointingCenterPixel, pointingCenterWorld);
467 0 : xSquintPixCache(0, chan) = pointingCenterPixel(0);
468 0 : ySquintPixCache(0, chan) = pointingCenterPixel(1);
469 : } else {
470 74 : xSquintPixCache(0, chan) = nonSquintedPointingPixel(0);
471 74 : ySquintPixCache(0, chan) = nonSquintedPointingPixel(1);
472 : }
473 74 : if (doSquint == BeamSquint::LL || doSquint == BeamSquint::GOFIGURE) {
474 0 : squint_p.getPointingDirection (pointDirE,
475 : parAngle,
476 0 : Quantity(spectralWorld(0),"Hz"),
477 : BeamSquint::LL, newpointDirE);
478 0 : pointingCenterWorld(0) = newpointDirE.getAngle().getValue("deg")(0);
479 0 : pointingCenterWorld(1) = newpointDirE.getAngle().getValue("deg")(1);
480 0 : directionCoord.toPixel(pointingCenterPixel, pointingCenterWorld);
481 0 : xSquintPixCache(1, chan) = pointingCenterPixel(0);
482 0 : ySquintPixCache(1, chan) = pointingCenterPixel(1);
483 : } else {
484 74 : xSquintPixCache(1, chan) = nonSquintedPointingPixel(0);
485 74 : 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 74 : IPosition ncs=in.shape();
501 74 : ncs(2) = 1; ncs(3) = 1;
502 148 : RO_LatticeIterator<Complex> li(in, LatticeStepper(in.shape(), ncs, IPosition(4,0,1,2,3) ) );
503 148 : 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 74 : Vector<Double> increment = directionCoord.increment();
511 74 : Int rrplane = -1;
512 74 : Int llplane = -1;
513 74 : stokesCoord.toPixel( rrplane, Stokes::RR );
514 74 : 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 74 : Int laststokes = -1;
524 74 : Int lastChan = -1;
525 : Int ichan;
526 : Int istokes;
527 : Int ix0, iy0;
528 : //Int indx;
529 148 : for(li.reset(),oli.reset();!li.atEnd();li++,oli++) {
530 :
531 74 : IPosition itsShape(li.matrixCursor().shape());
532 74 : IPosition loc(li.position());
533 :
534 74 : ichan = loc(3);
535 74 : istokes = loc(2);
536 74 : iy0 = loc(1);
537 74 : 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 74 : if ((doSquint == BeamSquint::RR) ||
545 0 : ((doSquint == BeamSquint::GOFIGURE) && (istokes == rrplane)) ) {
546 0 : xPixel = xSquintPixCache(0, ichan);
547 0 : yPixel = ySquintPixCache(0, ichan);
548 74 : } else if ((doSquint == BeamSquint::LL) ||
549 0 : ((doSquint == BeamSquint::GOFIGURE) && (istokes == llplane)) ) {
550 0 : xPixel = xSquintPixCache(1, ichan);
551 0 : yPixel = ySquintPixCache(1, ichan);
552 : } else {
553 74 : xPixel = nonSquintedPointingPixel(0);
554 74 : yPixel = nonSquintedPointingPixel(1);
555 : }
556 :
557 74 : if (istokes != laststokes) {
558 : //cerr << "Stokes = " << istokes << " pix = " << xPixel << ", " << yPixel << endl;
559 74 : laststokes = istokes;
560 : }
561 :
562 74 : Double factor = 60.0 * spectralCache(ichan)/1.0e+9 ; // arcminutes * GHz
563 74 : Double rmax2 = square( maximumRadius_p.getValue("'") / factor );
564 74 : 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 74 : Vector<Float> rx2(itsShape(0));
573 74 : Vector<Float> ry2(itsShape(1));
574 26186 : for(Int ix=0;ix<itsShape(0);ix++) {
575 26112 : rx2(ix) = square( increment(0)*((Double)(ix+ix0) - xPixel) );
576 : }
577 26186 : for(Int iy=0;iy<itsShape(1);iy++) {
578 26112 : ry2(iy) = square( increment(1)*((Double)(iy+iy0) - yPixel) );
579 : }
580 :
581 74 : const Matrix<Complex>& inmat = li.matrixCursor();
582 74 : Matrix<Complex>& outmat=oli.rwMatrixCursor();
583 :
584 : Bool incopy, outcopy, del;
585 74 : const Complex * inpoint = inmat.getStorage(incopy);
586 74 : Complex *outpoint =outmat.getStorage(outcopy);
587 74 : Float * rx2point = rx2.getStorage(del);
588 74 : Float * ry2point= ry2.getStorage(del);
589 74 : Complex* vppoint=vp_p.getStorage(del);
590 74 : Int nx=itsShape(0);
591 74 : Int ny=itsShape(1);
592 74 : Double inverseIncrementRadius=inverseIncrementRadius_p;
593 74 : #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 74 : outmat.putStorage(outpoint, outcopy);
643 74 : inmat.freeStorage(inpoint, incopy);
644 74 : }
645 :
646 74 : return out;
647 :
648 74 : };
649 :
650 26112 : 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 26112 : Complex taper;
655 21390848 : for(Int ix=0;ix<nx;ix++) {
656 :
657 21364736 : Float r2 = rx2[ix] + ry2;
658 :
659 21364736 : if (r2 > rmax2){
660 14793424 : out[ix+iy*nx] = 0.0;
661 : }
662 : else {
663 6571312 : r = sqrt(r2) * factor;
664 6571312 : indx = Int(r*inverseIncrementRadius);
665 6571312 : if (norm(vp[indx]) > 0.0) {
666 6571312 : if(ipower==2) {
667 6571312 : taper = vp[indx] * conj(vp[indx]);
668 : }
669 : else {
670 0 : taper = vp[indx];
671 : }
672 : } else {
673 0 : taper = 0.0;
674 : }
675 6571312 : 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 6571312 : if(!forward) {
682 0 : taper = conj(taper);
683 : }
684 6571312 : 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 6571312 : out[ix+iy*nx] = (in[ix+iy*nx]) * taper ;
692 : }
693 : }
694 : }
695 26112 : };
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 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 0 : Bool PBMath1D::ok()
1097 : {
1098 0 : if (vp_p.nelements() == 0) {
1099 0 : return false;
1100 0 : } else if (maximumRadius_p.getValue() <= 0.0) {
1101 0 : return false;
1102 0 : } else if (refFreq_p.getValue() <= 0.0) {
1103 0 : return false;
1104 0 : } else if (inverseIncrementRadius_p <= 0.0) {
1105 0 : return false;
1106 : } else {
1107 0 : 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 :
|
