Line data Source code
1 : //# GaussianShape.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: GaussianShape.cc 21130 2011-10-18 07:39:05Z gervandiepen $
27 :
28 : #include <components/ComponentModels/GaussianShape.h>
29 : #include <components/ComponentModels/Flux.h>
30 : #include <casacore/casa/Arrays/Vector.h>
31 : #include <casacore/casa/Exceptions/Error.h>
32 : #include <casacore/casa/Logging/LogIO.h>
33 : #include <casacore/casa/Logging/LogOrigin.h>
34 : #include <casacore/casa/BasicSL/Constants.h>
35 : #include <casacore/casa/BasicMath/Math.h>
36 : #include <casacore/measures/Measures/MCDirection.h>
37 : #include <casacore/measures/Measures/MDirection.h>
38 : #include <casacore/measures/Measures/MeasConvert.h>
39 : #include <casacore/measures/Measures/MeasRef.h>
40 : #include <casacore/casa/Quanta/MVAngle.h>
41 : #include <casacore/casa/Quanta/MVDirection.h>
42 : #include <casacore/casa/Quanta/Quantum.h>
43 : #include <casacore/casa/Utilities/Assert.h>
44 : #include <casacore/casa/BasicSL/String.h>
45 :
46 : using namespace casacore;
47 : namespace casa { //# NAMESPACE CASA - BEGIN
48 :
49 6504 : GaussianShape::GaussianShape()
50 : :TwoSidedShape(),
51 6504 : itsShape(1.0, 0.0, 0.0, Quantity(1,"'").getValue("rad"), 1.0, 0.0),
52 13008 : itsFT(itsShape)
53 : {
54 6504 : itsShape.setFlux(1.0);
55 6504 : updateFT();
56 6504 : DebugAssert(ok(), AipsError);
57 6504 : }
58 :
59 0 : GaussianShape::GaussianShape(const MDirection& direction,
60 : const Quantum<Double>& majorAxis,
61 : const Quantum<Double>& minorAxis,
62 0 : const Quantum<Double>& positionAngle)
63 0 : :TwoSidedShape(direction, majorAxis.getFullUnit(),
64 0 : minorAxis.getFullUnit(), positionAngle.getFullUnit()),
65 0 : itsShape(1.0, 0.0, 0.0, majorAxis.getValue("rad"),
66 0 : minorAxis.getValue("rad")/majorAxis.getValue("rad"),
67 0 : positionAngle.getValue("rad")),
68 0 : itsFT(itsShape)
69 : {
70 : // Adjust the flux of the Gaussian now that the width is correctly set
71 0 : itsShape.setFlux(1.0);
72 0 : updateFT();
73 0 : DebugAssert(ok(), AipsError);
74 0 : }
75 :
76 0 : GaussianShape::GaussianShape(const MDirection& direction,
77 : const Quantum<Double>& width,
78 : const Double axialRatio,
79 0 : const Quantum<Double>& positionAngle)
80 0 : :TwoSidedShape(direction, width.getFullUnit(),
81 0 : width.getFullUnit(), positionAngle.getFullUnit()),
82 0 : itsShape(1.0, 0.0, 0.0, width.getValue("rad"), axialRatio,
83 0 : positionAngle.getValue("rad")),
84 0 : itsFT(itsShape)
85 : {
86 : // Adjust the flux of the Gaussian now that the width is correctly set
87 0 : itsShape.setFlux(1.0);
88 0 : updateFT();
89 0 : DebugAssert(ok(), AipsError);
90 0 : }
91 :
92 7878 : GaussianShape::GaussianShape(const GaussianShape& other)
93 : :TwoSidedShape(other),
94 7878 : itsShape(other.itsShape),
95 7878 : itsFT(other.itsFT)
96 : {
97 7878 : DebugAssert(ok(), AipsError);
98 7878 : }
99 :
100 28764 : GaussianShape::~GaussianShape() {
101 14382 : DebugAssert(ok(), AipsError);
102 28764 : }
103 :
104 0 : GaussianShape& GaussianShape::operator=(const GaussianShape& other) {
105 0 : if (this != &other) {
106 0 : TwoSidedShape::operator=(other);
107 0 : itsShape = other.itsShape;
108 0 : itsFT = other.itsFT;
109 : }
110 0 : DebugAssert(ok(), AipsError);
111 0 : return *this;
112 : }
113 :
114 944 : ComponentType::Shape GaussianShape::type() const {
115 944 : DebugAssert(ok(), AipsError);
116 944 : return ComponentType::GAUSSIAN;
117 : }
118 :
119 6504 : void GaussianShape::setWidthInRad(const Double majorAxis,
120 : const Double minorAxis,
121 : const Double positionAngle) {
122 6504 : Vector<Double> angle(2);
123 6504 : angle(0) = majorAxis;
124 6504 : angle(1) = minorAxis;
125 6504 : itsShape.setWidth(angle);
126 6504 : itsShape.setPA(positionAngle);
127 : // Adjusting the width normally keeps the height constant and modifies the
128 : // flux. Modify this behaviour by restoring the flux
129 6504 : itsShape.setFlux(1.0);
130 6504 : updateFT();
131 6504 : DebugAssert(ok(), AipsError);
132 6504 : }
133 :
134 400 : Double GaussianShape::majorAxisInRad() const {
135 400 : DebugAssert(ok(), AipsError);
136 400 : return itsShape.majorAxis();
137 : }
138 :
139 400 : Double GaussianShape::minorAxisInRad() const {
140 400 : DebugAssert(ok(), AipsError);
141 400 : return itsShape.minorAxis();
142 : }
143 :
144 0 : Double GaussianShape::axialRatio() const {
145 0 : DebugAssert(ok(), AipsError);
146 0 : return itsShape.axialRatio();
147 : }
148 :
149 400 : Double GaussianShape::positionAngleInRad() const {
150 400 : DebugAssert(ok(), AipsError);
151 400 : return itsShape.PA();
152 : }
153 :
154 0 : Double GaussianShape::sample(const MDirection& direction,
155 : const MVAngle& pixelLatSize,
156 : const MVAngle& pixelLongSize) const {
157 0 : DebugAssert(ok(), AipsError);
158 0 : const MDirection& compDir(refDirection());
159 0 : const MDirection::Ref& compDirFrame(compDir.getRef());
160 0 : const MDirection::MVType* compDirValue = &(compDir.getValue());
161 0 : Bool deleteValue = false;
162 : // Convert direction to the same frame as the reference direction
163 0 : if (ComponentShape::differentRefs(direction.getRef(), compDirFrame)) {
164 0 : compDirValue = new MDirection::MVType
165 0 : (MDirection::Convert(compDir, direction.getRef())().getValue());
166 0 : deleteValue = true;
167 : }
168 0 : const MDirection::MVType& dirValue = direction.getValue();
169 0 : const Double separation = compDirValue->separation(dirValue);
170 0 : Double retVal = 0.0;
171 0 : if (separation < 4 * itsShape.majorAxis()) {
172 0 : const Double pa = - compDirValue->positionAngle(dirValue);
173 0 : retVal = pixelLatSize.radian() * pixelLongSize.radian() *
174 0 : itsShape(separation*sin(pa), separation*cos(pa));
175 : }
176 0 : if (deleteValue) delete compDirValue;
177 0 : return retVal;
178 0 : }
179 :
180 4278 : void GaussianShape::sample(Vector<Double>& scale,
181 : const Vector<MDirection::MVType>& directions,
182 : const MDirection::Ref& refFrame,
183 : const MVAngle& pixelLatSize,
184 : const MVAngle& pixelLongSize) const {
185 4278 : DebugAssert(ok(), AipsError);
186 4278 : const uInt nSamples = directions.nelements();
187 4278 : DebugAssert(scale.nelements() == nSamples, AipsError);
188 :
189 4278 : const MDirection& compDir(refDirection());
190 4278 : const MDirection::Ref& compDirFrame(compDir.getRef());
191 4278 : const MDirection::MVType* compDirValue = &(compDir.getValue());
192 4278 : Bool deleteValue = false;
193 : // Convert direction to the same frame as the reference direction
194 4278 : if (refFrame != compDirFrame) {
195 4278 : compDirValue = new MDirection::MVType
196 4278 : (MDirection::Convert(compDir, refFrame)().getValue());
197 4278 : deleteValue = true;
198 : }
199 4278 : const Double pixArea = pixelLatSize.radian() * pixelLongSize.radian();
200 4278 : const Double maxSep = 4.0 * itsShape.majorAxis();
201 : Double separation, pa;
202 131595330 : for (uInt i = 0; i < nSamples; i++) {
203 131591052 : const MDirection::MVType& dirVal = directions(i);
204 131591052 : separation = compDirValue->separation(dirVal);
205 131591052 : if (separation > maxSep) {
206 126735604 : scale(i) = 0.0;
207 : } else {
208 4855448 : pa = - compDirValue->positionAngle(dirVal);
209 4855448 : scale(i) = pixArea * itsShape(separation*sin(pa), separation*cos(pa));
210 : }
211 : }
212 4278 : if (deleteValue) delete compDirValue;
213 4278 : }
214 :
215 600250 : DComplex GaussianShape::visibility(const Vector<Double>& uvw,
216 : const Double& frequency) const {
217 600250 : DebugAssert(uvw.nelements() == 3, AipsError);
218 600250 : DebugAssert(frequency > 0, AipsError);
219 600250 : DebugAssert(ok(), AipsError);
220 600250 : const Double wavenumber = frequency/C::c;
221 600250 : return DComplex(itsFT(-uvw(0)*wavenumber, uvw(1)*wavenumber), 0.0);
222 : }
223 :
224 0 : void GaussianShape::visibility(Vector<DComplex>& scale,
225 : const Matrix<Double>& uvw,
226 : const Double& frequency) const {
227 0 : ComponentShape::visibility(scale, uvw, frequency);
228 0 : }
229 :
230 6272 : void GaussianShape::visibility(Matrix<DComplex>& scale,
231 : const Matrix<Double>& uvw,
232 : const Vector<Double>& frequency) const {
233 6272 : ComponentShape::visibility(scale, uvw, frequency);
234 6272 : }
235 :
236 7878 : ComponentShape* GaussianShape::clone() const {
237 7878 : DebugAssert(ok(), AipsError);
238 7878 : ComponentShape* tmpPtr = new GaussianShape(*this);
239 7878 : AlwaysAssert(tmpPtr != 0, AipsError);
240 7878 : return tmpPtr;
241 : }
242 :
243 665342 : Bool GaussianShape::ok() const {
244 : // The LogIO class is only constructed if an error is detected for
245 : // performance reasons. Both function static and file static variables
246 : // where considered and rejected for this purpose.
247 665342 : if (!TwoSidedShape::ok()) return false;
248 665342 : if (!near(itsShape.flux(), Double(1.0), 2*C::dbl_epsilon)) {
249 0 : LogIO logErr(LogOrigin("GaussianCompRep", "ok()"));
250 : logErr << LogIO::SEVERE << "The internal Gaussian shape does not have"
251 : << " unit area"
252 0 : << LogIO::POST;
253 0 : return false;
254 0 : }
255 665342 : if (!near(itsFT.height(), 1.0, 2*C::dbl_epsilon)) {
256 0 : LogIO logErr(LogOrigin("GaussianCompRep", "ok()"));
257 : logErr << LogIO::SEVERE << "The cached Fourier Transform of"
258 : << " the internal Gaussian shape does not have"
259 : << " unit height"
260 0 : << LogIO::POST;
261 0 : return false;
262 0 : }
263 665342 : return true;
264 : }
265 :
266 78 : const ComponentShape* GaussianShape::getPtr() const {
267 78 : return this;
268 : }
269 :
270 0 : Quantity GaussianShape::getArea() const {
271 0 : Double majorAxis = itsShape.majorAxis();
272 0 : Double minorAxis = itsShape.minorAxis();
273 :
274 0 : Quantity area(C::pi/(4*C::ln2) * majorAxis * minorAxis, "sr");
275 0 : return area;
276 : }
277 :
278 13008 : void GaussianShape::updateFT() {
279 13008 : const Double factor = 4.0*C::ln2/C::pi;
280 13008 : Vector<Double> width(2);
281 13008 : width(0) = factor/itsShape.minorAxis();
282 13008 : width(1) = factor/itsShape.majorAxis();
283 13008 : itsFT.setWidth(width);
284 13008 : itsFT.setPA(itsShape.PA() + C::pi_2);
285 13008 : }
286 :
287 0 : String GaussianShape::sizeToString() const {
288 : return TwoSidedShape::sizeToString(
289 0 : Quantity(itsShape.majorAxis(), "rad"),
290 0 : Quantity(itsShape.minorAxis(), "rad"),
291 0 : Quantity(itsShape.PA(), "rad"), true,
292 : majorAxisError(), minorAxisError(),
293 : positionAngleError()
294 0 : );
295 : }
296 :
297 :
298 :
299 :
300 :
301 : // Local Variables:
302 : // compile-command: "gmake GaussianShape"
303 : // End:
304 :
305 : } //# NAMESPACE CASA - END
306 :
|
