Line data Source code
1 : //# DiskShape.cc:
2 : //# Copyright (C) 1998,1999,2000
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 addressed 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 : //# $Id: DiskShape.cc 21130 2011-10-18 07:39:05Z gervandiepen $
27 :
28 : #include <components/ComponentModels/DiskShape.h>
29 : #include <components/ComponentModels/Flux.h>
30 : #include <casacore/casa/Arrays/Matrix.h>
31 : #include <casacore/casa/Arrays/Vector.h>
32 : #include <casacore/casa/Arrays/VectorIter.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/casa/Arrays/Array.h>
35 : #include <casacore/casa/Exceptions/Error.h>
36 : #include <casacore/casa/Logging/LogIO.h>
37 : #include <casacore/casa/Logging/LogOrigin.h>
38 : #include <casacore/casa/BasicSL/Constants.h>
39 : #include <casacore/casa/BasicMath/Math.h>
40 : #include <casacore/measures/Measures/MCDirection.h>
41 : #include <casacore/measures/Measures/MDirection.h>
42 : #include <casacore/measures/Measures/MeasConvert.h>
43 : #include <casacore/measures/Measures/MeasRef.h>
44 : #include <casacore/casa/Quanta/MVAngle.h>
45 : #include <casacore/casa/Quanta/MVDirection.h>
46 : #include <casacore/casa/Quanta/Quantum.h>
47 : #include <casacore/casa/Utilities/Assert.h>
48 : #include <casacore/casa/BasicSL/String.h>
49 : #include <cmath>
50 :
51 : using namespace casacore;
52 : namespace casa { //# NAMESPACE CASA - BEGIN
53 :
54 13597 : DiskShape::DiskShape()
55 : :TwoSidedShape(),
56 13597 : _majorAxis(Quantity(1,"'").getValue("rad")),
57 13597 : _minorAxis(Quantity(1,"'").getValue("rad")),
58 13597 : _pa(Quantity(0,"deg").getValue("rad")),
59 27194 : _recipArea(1.0/_area())
60 : {
61 13597 : DebugAssert(ok(), AipsError);
62 13597 : }
63 :
64 0 : DiskShape::DiskShape(const MDirection& direction,
65 : const Quantum<Double>& majorAxis,
66 : const Quantum<Double>& minorAxis,
67 0 : const Quantum<Double>& positionAngle)
68 0 : :TwoSidedShape(direction, majorAxis.getFullUnit(),
69 0 : minorAxis.getFullUnit(), positionAngle.getFullUnit()),
70 0 : _majorAxis(majorAxis.getValue("rad")),
71 0 : _minorAxis(minorAxis.getValue("rad")),
72 0 : _pa(positionAngle.getValue("rad")),
73 0 : _recipArea(1.0/_area())
74 : {
75 0 : DebugAssert(ok(), AipsError);
76 0 : }
77 :
78 0 : DiskShape::DiskShape(const MDirection& direction,
79 : const Quantum<Double>& width,
80 : const Double axialRatio,
81 0 : const Quantum<Double>& positionAngle)
82 0 : :TwoSidedShape(direction, width.getFullUnit(),
83 0 : width.getFullUnit(), positionAngle.getFullUnit()),
84 0 : _majorAxis(width.getValue("rad")),
85 0 : _minorAxis(_majorAxis*axialRatio),
86 0 : _pa(positionAngle.getValue("rad")),
87 0 : _recipArea(1.0/_area())
88 : {
89 0 : DebugAssert(ok(), AipsError);
90 0 : }
91 :
92 14028 : DiskShape::DiskShape(const DiskShape& other)
93 : :TwoSidedShape(other),
94 14028 : _majorAxis(other._majorAxis),
95 14028 : _minorAxis(other._minorAxis),
96 14028 : _pa(other._pa),
97 14028 : _recipArea(other._recipArea)
98 : {
99 14028 : DebugAssert(ok(), AipsError);
100 14028 : }
101 :
102 55236 : DiskShape::~DiskShape() {
103 27625 : DebugAssert(ok(), AipsError);
104 55236 : }
105 :
106 0 : DiskShape& DiskShape::operator=(const DiskShape& other) {
107 0 : if (this != &other) {
108 0 : TwoSidedShape::operator=(other);
109 0 : _majorAxis = other._majorAxis;
110 0 : _minorAxis = other._minorAxis;
111 0 : _pa = other._pa;
112 0 : _recipArea = other._recipArea;
113 : }
114 0 : DebugAssert(ok(), AipsError);
115 0 : return *this;
116 : }
117 :
118 68136 : ComponentType::Shape DiskShape::type() const {
119 68136 : DebugAssert(ok(), AipsError);
120 68136 : return ComponentType::DISK;
121 : }
122 :
123 13607 : void DiskShape::setWidthInRad(const Double majorAxis,
124 : const Double minorAxis,
125 : const Double positionAngle) {
126 : static const double epsilon = 1.0e-17;
127 13607 : _majorAxis = majorAxis;
128 13607 : _minorAxis = minorAxis;
129 13607 : _minorAxis = std::abs(majorAxis-minorAxis) < epsilon ? majorAxis : minorAxis;
130 13607 : _pa = positionAngle;
131 13607 : AlwaysAssert(_majorAxis > 0 && _minorAxis > 0 && _majorAxis >=_minorAxis,
132 : AipsError);
133 13607 : _recipArea = 1.0/_area();
134 13607 : DebugAssert(ok(), AipsError);
135 13607 : }
136 :
137 667 : Double DiskShape::majorAxisInRad() const {
138 667 : DebugAssert(ok(), AipsError);
139 667 : return _majorAxis;
140 : }
141 :
142 667 : Double DiskShape::minorAxisInRad() const {
143 667 : DebugAssert(ok(), AipsError);
144 667 : return _minorAxis;
145 : }
146 :
147 667 : Double DiskShape::positionAngleInRad() const {
148 667 : DebugAssert(ok(), AipsError);
149 667 : return _pa;
150 : }
151 :
152 0 : Double DiskShape::sample(const MDirection& direction,
153 : const MVAngle& pixelLatSize,
154 : const MVAngle& pixelLongSize) const {
155 0 : DebugAssert(ok(), AipsError);
156 0 : const MDirection& compDir(refDirection());
157 0 : const MDirection::Ref& compDirFrame(compDir.getRef());
158 0 : const MDirection::MVType* compDirValue = &(compDir.getValue());
159 0 : Bool deleteValue = false;
160 : // Convert direction to the same frame as the reference direction
161 0 : if (ComponentShape::differentRefs(direction.getRef(), compDirFrame)) {
162 0 : compDirValue = new MDirection::MVType
163 0 : (MDirection::Convert(compDir, direction.getRef())().getValue());
164 0 : deleteValue = true;
165 : }
166 0 : Double retVal = _calcSample(*compDirValue, direction.getValue(),
167 0 : _majorAxis/2.0, _minorAxis/2.0,
168 0 : _recipArea*pixelLatSize.radian()*
169 0 : pixelLongSize.radian());
170 0 : if (deleteValue) delete compDirValue;
171 0 : return retVal;
172 0 : }
173 :
174 0 : void DiskShape::sample(Vector<Double>& scale,
175 : const Vector<MVDirection>& directions,
176 : const MDirection::Ref& refFrame,
177 : const MVAngle& pixelLatSize,
178 : const MVAngle& pixelLongSize) const {
179 0 : DebugAssert(ok(), AipsError);
180 0 : const uInt nSamples = directions.nelements();
181 0 : DebugAssert(scale.nelements() == nSamples, AipsError);
182 :
183 0 : const MDirection& compDir(refDirection());
184 0 : const MDirection::Ref& compDirFrame(compDir.getRef());
185 0 : const MDirection::MVType* compDirValue = &(compDir.getValue());
186 0 : Bool deleteValue = false;
187 : // Convert direction to the same frame as the reference direction
188 0 : if (refFrame != compDirFrame) {
189 0 : compDirValue = new MDirection::MVType
190 0 : (MDirection::Convert(compDir, refFrame)().getValue());
191 0 : deleteValue = true;
192 : }
193 0 : const Double majRad = _majorAxis/2.0;
194 0 : const Double minRad = _minorAxis/2.0;
195 0 : const Double pixValue = _recipArea *
196 0 : pixelLatSize.radian() * pixelLongSize.radian();
197 0 : for (uInt i = 0; i < nSamples; i++) {
198 0 : scale(i) = _calcSample(*compDirValue, directions(i),
199 : majRad, minRad, pixValue);
200 : }
201 0 : if (deleteValue) delete compDirValue;
202 0 : }
203 :
204 67375 : DComplex DiskShape::visibility(const Vector<Double>& uvw,
205 : const Double& frequency) const {
206 67375 : DebugAssert(uvw.nelements() == 3, AipsError);
207 67375 : DebugAssert(frequency > 0, AipsError);
208 67375 : DebugAssert(ok(), AipsError);
209 67375 : Double u = uvw(0);
210 67375 : Double v = uvw(1);
211 67375 : if (near(u + v, 0.0)) return DComplex(1.0, 0.0);
212 67375 : if (!nearAbs(_pa, 0.0)) {
213 67375 : _rotateVis(u, v, cos(_pa), sin(_pa));
214 : }
215 67375 : return DComplex(_calcVis(u, v, C::pi * frequency/C::c), 0.0);
216 : }
217 :
218 13460 : void DiskShape::visibility(Matrix<DComplex>& scale,
219 : const Matrix<Double>& uvw,
220 : const Vector<Double>& frequency) const {
221 :
222 13460 : scale.resize(uvw.ncolumn(), frequency.nelements());
223 :
224 13460 : VectorIterator<DComplex> iter(scale, 0);
225 19597828 : for ( uInt k =0 ; k < frequency.nelements() ; ++k){
226 19584368 : visibility(iter.vector(), uvw, frequency(k));
227 19584368 : iter.next();
228 : }
229 13460 : }
230 :
231 19584368 : void DiskShape::visibility(Vector<DComplex>& scale,
232 : const Matrix<Double>& uvw,
233 : const Double& frequency) const {
234 19584368 : DebugAssert(ok(), AipsError);
235 19584368 : const uInt nSamples = scale.nelements();
236 19584368 : DebugAssert(uvw.ncolumn() == nSamples, AipsError);
237 19584368 : DebugAssert(uvw.nrow() == 3, AipsError);
238 19584368 : DebugAssert(frequency > 0, AipsError);
239 :
240 19584368 : Bool doRotation = false;
241 19584368 : Double cpa = 1.0, spa = 0.0;
242 19584368 : if (!nearAbs(_pa, 0.0)) {
243 11842232 : doRotation = true;
244 11842232 : cpa = cos(_pa);
245 11842232 : spa = sin(_pa);
246 : }
247 :
248 19584368 : const Double factor = C::pi * frequency/C::c;
249 : Double u, v;
250 52622240 : for (uInt i = 0; i < nSamples; i++) {
251 33037872 : u = uvw(0, i);
252 33037872 : v = uvw(1, i);
253 : // DComplex& thisVis = scale(i);
254 : /// thisVis.imag() = 0.0;
255 33037872 : if (near(u + v, 0.0)) {
256 : /// thisVis.real() = 1.0; // avoids dividing by zero in _calcVis(...)
257 7782400 : scale[i] = DComplex(1.0, 0.0); // avoids dividing by zero
258 : // in _calcVis(...)
259 : } else {
260 25255472 : if (doRotation) _rotateVis(u, v, cpa, spa);
261 : /// thisVis.real() = _calcVis(u, v, factor);
262 25255472 : scale[i] = DComplex(_calcVis(u, v, factor), 0.0);
263 : }
264 : }
265 19584368 : }
266 :
267 14028 : ComponentShape* DiskShape::clone() const {
268 14028 : DebugAssert(ok(), AipsError);
269 14028 : ComponentShape* tmpPtr = new DiskShape(*this);
270 14028 : AlwaysAssert(tmpPtr != 0, AipsError);
271 14028 : return tmpPtr;
272 : }
273 :
274 19833450 : Bool DiskShape::ok() const {
275 : // The LogIO class is only constructed if an error is detected for
276 : // performance reasons. Both function static and file static variables
277 : // where considered and rejected for this purpose.
278 19833450 : if (!TwoSidedShape::ok()) return false;
279 19833450 : if (_majorAxis <= 0) {
280 0 : LogIO logErr(LogOrigin("DiskCompRep", "ok()"));
281 : logErr << LogIO::SEVERE << "The major axis width is zero or negative"
282 0 : << LogIO::POST;
283 0 : return false;
284 0 : }
285 19833450 : if (_minorAxis <= 0) {
286 0 : LogIO logErr(LogOrigin("DiskCompRep", "ok()"));
287 : logErr << LogIO::SEVERE << "The minor axis width is zero or negative"
288 0 : << LogIO::POST;
289 0 : return false;
290 0 : }
291 19833450 : if (_minorAxis > _majorAxis) {
292 0 : LogIO logErr(LogOrigin("DiskCompRep", "ok()"));
293 : logErr << LogIO::SEVERE << "The minor axis width is larger than "
294 : << "the major axis width"
295 0 : << LogIO::POST;
296 0 : return false;
297 0 : }
298 19833450 : if (!near(_recipArea, 1.0/_area(), 2*C::dbl_epsilon)) {
299 0 : LogIO logErr(LogOrigin("DiskCompRep", "ok()"));
300 : logErr << LogIO::SEVERE << "The disk shape does not have"
301 : << " unit area"
302 0 : << LogIO::POST;
303 0 : return false;
304 0 : }
305 19833450 : return true;
306 : }
307 :
308 0 : const ComponentShape* DiskShape::getPtr() const {
309 0 : return this;
310 : }
311 :
312 25322847 : Double DiskShape::_calcVis(Double u, Double v, const Double factor) const {
313 25322847 : u *= _minorAxis;
314 25322847 : v *= _majorAxis;
315 25322847 : const Double r = hypot(u, v) * factor;
316 25322847 : return 2.0 * j1(r)/r;
317 : }
318 :
319 15368391 : void DiskShape::_rotateVis(Double& u, Double& v,
320 : const Double cpa, const Double spa) {
321 15368391 : const Double utemp = u;
322 15368391 : u = u * cpa - v * spa;
323 15368391 : v = utemp * spa + v * cpa;
324 15368391 : }
325 :
326 19860654 : Double DiskShape::_area() const {
327 19860654 : return C::pi_4 * _majorAxis * _minorAxis;
328 : }
329 :
330 0 : Double DiskShape::_calcSample(const MDirection::MVType& compDirValue,
331 : const MDirection::MVType& dirVal,
332 : const Double majRad, const Double minRad,
333 : const Double pixValue) const {
334 0 : const Double separation = compDirValue.separation(dirVal);
335 0 : if (separation <= majRad) {
336 0 : const Double pa = compDirValue.positionAngle(dirVal) - _pa;
337 0 : const Double x = abs(separation*cos(pa));
338 0 : const Double y = abs(separation*sin(pa));
339 0 : if ((x <= majRad) &&
340 0 : (y <= minRad) &&
341 0 : (y <= minRad * sqrt(1 - square(x/majRad)))) {
342 0 : return pixValue;
343 : }
344 : }
345 0 : return 0.0;
346 : }
347 :
348 0 : String DiskShape::sizeToString() const {
349 : return TwoSidedShape::sizeToString(
350 0 : Quantity(_majorAxis, "rad"),
351 0 : Quantity(_minorAxis, "rad"),
352 0 : Quantity(_pa, "rad"),
353 : false
354 0 : );
355 : }
356 :
357 : }
358 :
|
