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: aips2-request@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 1261 : GaussianShape::GaussianShape()
50 : :TwoSidedShape(),
51 2522 : itsShape(1.0, 0.0, 0.0, Quantity(1,"'").getValue("rad"), 1.0, 0.0),
52 3783 : itsFT(itsShape)
53 : {
54 1261 : itsShape.setFlux(1.0);
55 1261 : updateFT();
56 1261 : DebugAssert(ok(), AipsError);
57 1261 : }
58 :
59 722 : GaussianShape::GaussianShape(const MDirection& direction,
60 : const Quantum<Double>& majorAxis,
61 : const Quantum<Double>& minorAxis,
62 722 : const Quantum<Double>& positionAngle)
63 722 : :TwoSidedShape(direction, majorAxis.getFullUnit(),
64 1444 : minorAxis.getFullUnit(), positionAngle.getFullUnit()),
65 724 : itsShape(1.0, 0.0, 0.0, majorAxis.getValue("rad"),
66 1446 : minorAxis.getValue("rad")/majorAxis.getValue("rad"),
67 1445 : positionAngle.getValue("rad")),
68 2166 : itsFT(itsShape)
69 : {
70 : // Adjust the flux of the Gaussian now that the width is correctly set
71 721 : itsShape.setFlux(1.0);
72 721 : updateFT();
73 721 : DebugAssert(ok(), AipsError);
74 721 : }
75 :
76 1 : GaussianShape::GaussianShape(const MDirection& direction,
77 : const Quantum<Double>& width,
78 : const Double axialRatio,
79 1 : const Quantum<Double>& positionAngle)
80 1 : :TwoSidedShape(direction, width.getFullUnit(),
81 2 : width.getFullUnit(), positionAngle.getFullUnit()),
82 1 : itsShape(1.0, 0.0, 0.0, width.getValue("rad"), axialRatio,
83 2 : positionAngle.getValue("rad")),
84 3 : itsFT(itsShape)
85 : {
86 : // Adjust the flux of the Gaussian now that the width is correctly set
87 1 : itsShape.setFlux(1.0);
88 1 : updateFT();
89 1 : DebugAssert(ok(), AipsError);
90 1 : }
91 :
92 3774 : GaussianShape::GaussianShape(const GaussianShape& other)
93 : :TwoSidedShape(other),
94 3774 : itsShape(other.itsShape),
95 3774 : itsFT(other.itsFT)
96 : {
97 3774 : DebugAssert(ok(), AipsError);
98 3774 : }
99 :
100 10266 : GaussianShape::~GaussianShape() {
101 5757 : DebugAssert(ok(), AipsError);
102 10266 : }
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 72035 : ComponentType::Shape GaussianShape::type() const {
115 72035 : DebugAssert(ok(), AipsError);
116 72035 : return ComponentType::GAUSSIAN;
117 : }
118 :
119 1425 : void GaussianShape::setWidthInRad(const Double majorAxis,
120 : const Double minorAxis,
121 : const Double positionAngle) {
122 2850 : Vector<Double> angle(2);
123 1425 : angle(0) = majorAxis;
124 1425 : angle(1) = minorAxis;
125 1425 : itsShape.setWidth(angle);
126 1425 : 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 1425 : itsShape.setFlux(1.0);
130 1425 : updateFT();
131 1425 : DebugAssert(ok(), AipsError);
132 1425 : }
133 :
134 3586 : Double GaussianShape::majorAxisInRad() const {
135 3586 : DebugAssert(ok(), AipsError);
136 3586 : return itsShape.majorAxis();
137 : }
138 :
139 3586 : Double GaussianShape::minorAxisInRad() const {
140 3586 : DebugAssert(ok(), AipsError);
141 3586 : return itsShape.minorAxis();
142 : }
143 :
144 0 : Double GaussianShape::axialRatio() const {
145 0 : DebugAssert(ok(), AipsError);
146 0 : return itsShape.axialRatio();
147 : }
148 :
149 3528 : Double GaussianShape::positionAngleInRad() const {
150 3528 : DebugAssert(ok(), AipsError);
151 3528 : 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 : }
179 :
180 4438 : 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 4438 : DebugAssert(ok(), AipsError);
186 4438 : const uInt nSamples = directions.nelements();
187 4438 : DebugAssert(scale.nelements() == nSamples, AipsError);
188 :
189 4438 : const MDirection& compDir(refDirection());
190 8876 : const MDirection::Ref& compDirFrame(compDir.getRef());
191 4438 : const MDirection::MVType* compDirValue = &(compDir.getValue());
192 4438 : Bool deleteValue = false;
193 : // Convert direction to the same frame as the reference direction
194 4438 : if (refFrame != compDirFrame) {
195 4438 : compDirValue = new MDirection::MVType
196 4438 : (MDirection::Convert(compDir, refFrame)().getValue());
197 4438 : deleteValue = true;
198 : }
199 4438 : const Double pixArea = pixelLatSize.radian() * pixelLongSize.radian();
200 4438 : const Double maxSep = 4.0 * itsShape.majorAxis();
201 : Double separation, pa;
202 135765490 : for (uInt i = 0; i < nSamples; i++) {
203 135761052 : const MDirection::MVType& dirVal = directions(i);
204 135761052 : separation = compDirValue->separation(dirVal);
205 135761052 : if (separation > maxSep) {
206 130264456 : scale(i) = 0.0;
207 : } else {
208 5496596 : pa = - compDirValue->positionAngle(dirVal);
209 5496596 : scale(i) = pixArea * itsShape(separation*sin(pa), separation*cos(pa));
210 : }
211 : }
212 4438 : if (deleteValue) delete compDirValue;
213 4438 : }
214 :
215 667625 : DComplex GaussianShape::visibility(const Vector<Double>& uvw,
216 : const Double& frequency) const {
217 667625 : DebugAssert(uvw.nelements() == 3, AipsError);
218 667625 : DebugAssert(frequency > 0, AipsError);
219 667625 : DebugAssert(ok(), AipsError);
220 667625 : const Double wavenumber = frequency/C::c;
221 667625 : 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 196 : void GaussianShape::visibility(Matrix<DComplex>& scale,
231 : const Matrix<Double>& uvw,
232 : const Vector<Double>& frequency) const {
233 196 : ComponentShape::visibility(scale, uvw, frequency);
234 196 : }
235 :
236 3774 : ComponentShape* GaussianShape::clone() const {
237 3774 : DebugAssert(ok(), AipsError);
238 3774 : ComponentShape* tmpPtr = new GaussianShape(*this);
239 3774 : AlwaysAssert(tmpPtr != 0, AipsError);
240 3774 : return tmpPtr;
241 : }
242 :
243 777680 : 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 777680 : if (!TwoSidedShape::ok()) return false;
248 777680 : 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 : }
255 777680 : 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 : }
263 777680 : return true;
264 : }
265 :
266 1492 : const ComponentShape* GaussianShape::getPtr() const {
267 1492 : return this;
268 : }
269 :
270 344 : Quantity GaussianShape::getArea() const {
271 344 : Double majorAxis = itsShape.majorAxis();
272 344 : Double minorAxis = itsShape.minorAxis();
273 :
274 344 : Quantity area(C::pi/(4*C::ln2) * majorAxis * minorAxis, "sr");
275 344 : return area;
276 : }
277 :
278 3408 : void GaussianShape::updateFT() {
279 3408 : const Double factor = 4.0*C::ln2/C::pi;
280 3408 : Vector<Double> width(2);
281 3408 : width(0) = factor/itsShape.minorAxis();
282 3408 : width(1) = factor/itsShape.majorAxis();
283 3408 : itsFT.setWidth(width);
284 3408 : itsFT.setPA(itsShape.PA() + C::pi_2);
285 3408 : }
286 :
287 393 : String GaussianShape::sizeToString() const {
288 : return TwoSidedShape::sizeToString(
289 786 : Quantity(itsShape.majorAxis(), "rad"),
290 786 : Quantity(itsShape.minorAxis(), "rad"),
291 786 : Quantity(itsShape.PA(), "rad"), true,
292 : majorAxisError(), minorAxisError(),
293 : positionAngleError()
294 1179 : );
295 : }
296 :
297 :
298 :
299 :
300 :
301 : // Local Variables:
302 : // compile-command: "gmake GaussianShape"
303 : // End:
304 :
305 : } //# NAMESPACE CASA - END
306 :
|
