Line data Source code
1 : //# TwoSidedShape.cc:
2 : //# Copyright (C) 1999,2000,2001,2002
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: TwoSidedShape.cc 21292 2012-11-28 14:58:19Z gervandiepen $
27 :
28 : #include <components/ComponentModels/TwoSidedShape.h>
29 : #include <iomanip>
30 : #include <casacore/casa/Arrays/Vector.h>
31 : #include <casacore/casa/Arrays/ArrayLogical.h>
32 : #include <casacore/casa/Containers/Record.h>
33 : #include <casacore/casa/Containers/RecordFieldId.h>
34 : #include <casacore/casa/Containers/RecordInterface.h>
35 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
36 : #include <casacore/casa/Exceptions/Error.h>
37 : #include <casacore/casa/Logging/LogIO.h>
38 : #include <casacore/casa/Logging/LogOrigin.h>
39 : #include <casacore/casa/BasicSL/Constants.h>
40 : #include <casacore/casa/BasicMath/Math.h>
41 : #include <casacore/casa/Quanta/Quantum.h>
42 : #include <casacore/casa/Quanta/QuantumHolder.h>
43 : #include <casacore/casa/Quanta/MVAngle.h>
44 : #include <casacore/casa/Utilities/Assert.h>
45 : #include <casacore/casa/Utilities/Precision.h>
46 : #include <casacore/casa/BasicSL/String.h>
47 :
48 : using namespace casacore;
49 : namespace casa { //# NAMESPACE CASA - BEGIN
50 :
51 14382 : TwoSidedShape::~TwoSidedShape() {
52 14382 : DebugAssert(ok(), AipsError);
53 14382 : }
54 :
55 0 : TwoSidedShape& TwoSidedShape::operator=(const TwoSidedShape& other) {
56 0 : if (this != &other) {
57 0 : ComponentShape::operator=(other);
58 0 : itsMajUnit = other.itsMajUnit;
59 0 : itsMinUnit = other.itsMinUnit;
60 0 : itsPaUnit = other.itsPaUnit;
61 0 : itsMajErr = other.itsMajErr;
62 0 : itsMinErr = other.itsMinErr;
63 0 : itsPaErr = other.itsPaErr;
64 : }
65 0 : DebugAssert(ok(), AipsError);
66 0 : return *this;
67 : }
68 :
69 6350 : void TwoSidedShape::setWidth(const Quantity& majorAxis,
70 : const Quantity& minorAxis,
71 : const Quantity& positionAngle) {
72 6350 : itsMajUnit = majorAxis.getFullUnit();
73 6350 : itsMinUnit = minorAxis.getFullUnit();
74 6350 : itsPaUnit = positionAngle.getFullUnit();
75 6350 : const Unit rad("rad");
76 6350 : setWidthInRad(majorAxis.getValue(rad), minorAxis.getValue(rad),
77 : positionAngle.getValue(rad));
78 6350 : DebugAssert(ok(), AipsError);
79 6350 : }
80 :
81 0 : void TwoSidedShape::setWidth(const Quantum<Double>& majorAxis,
82 : const Double axialRatio,
83 : const Quantum<Double>& positionAngle) {
84 0 : itsMinUnit = itsMajUnit = majorAxis.getFullUnit();
85 0 : itsPaUnit = positionAngle.getFullUnit();
86 0 : const Unit rad("rad");
87 0 : const Double majWidth = majorAxis.getValue(rad);
88 0 : setWidthInRad(majWidth, majWidth*axialRatio, positionAngle.getValue(rad));
89 0 : DebugAssert(ok(), AipsError);
90 0 : }
91 :
92 274 : Quantum<Double> TwoSidedShape::majorAxis() const {
93 274 : Quantum<Double> retVal(majorAxisInRad(), Unit("rad"));
94 274 : retVal.convert(itsMajUnit);
95 274 : return retVal;
96 0 : }
97 :
98 274 : Quantum<Double> TwoSidedShape::minorAxis() const {
99 274 : Quantum<Double> retVal(minorAxisInRad(), Unit("rad"));
100 274 : retVal.convert(itsMinUnit);
101 274 : return retVal;
102 0 : }
103 :
104 274 : Quantum<Double> TwoSidedShape::positionAngle() const {
105 274 : Quantum<Double> retVal(positionAngleInRad(), Unit("rad"));
106 274 : retVal.convert(itsPaUnit);
107 274 : return retVal;
108 0 : }
109 :
110 0 : Double TwoSidedShape::axialRatio() const {
111 0 : return minorAxisInRad()/majorAxisInRad();
112 : }
113 :
114 6350 : void TwoSidedShape::setErrors(const Quantum<Double>& majorAxisError,
115 : const Quantum<Double>& minorAxisError,
116 : const Quantum<Double>& positionAngleError) {
117 6350 : if (ComponentShape::badError(majorAxisError) ||
118 12700 : ComponentShape::badError(minorAxisError) ||
119 6350 : ComponentShape::badError(positionAngleError)) {
120 0 : LogIO logErr(LogOrigin("TwoSidedShape", "setErrors(...)"));
121 : logErr << "The errors must be non-negative angular quantities."
122 0 : << LogIO::EXCEPTION;
123 0 : }
124 6350 : itsMajErr = majorAxisError;
125 6350 : itsMinErr = minorAxisError;
126 6350 : itsPaErr = positionAngleError;
127 6350 : }
128 :
129 274 : const Quantum<Double>& TwoSidedShape::majorAxisError() const {
130 274 : return itsMajErr;
131 : }
132 :
133 274 : const Quantum<Double>& TwoSidedShape::minorAxisError() const {
134 274 : return itsMinErr;
135 : }
136 :
137 274 : const Quantum<Double>& TwoSidedShape::positionAngleError() const {
138 274 : return itsPaErr;
139 : }
140 :
141 0 : Double TwoSidedShape::axialRatioError() const {
142 0 : const Unit rad("rad");
143 0 : const Double relErr = itsMajErr.getValue(rad)/majorAxisInRad() +
144 0 : itsMinErr.getValue(rad)/minorAxisInRad();
145 0 : return axialRatio() * relErr;
146 0 : }
147 :
148 0 : void TwoSidedShape::sample(Vector<Double>& scale,
149 : const Vector<MDirection::MVType>& directions,
150 : const MDirection::Ref& refFrame,
151 : const MVAngle& pixelLatSize,
152 : const MVAngle& pixelLongSize) const {
153 0 : ComponentShape::sample(scale, directions, refFrame, pixelLatSize,
154 : pixelLongSize);
155 0 : }
156 :
157 0 : void TwoSidedShape::visibility(Vector<DComplex>& scale,
158 : const Matrix<Double>& uvw,
159 : const Double& frequency) const {
160 0 : ComponentShape::visibility(scale, uvw, frequency);
161 0 : }
162 :
163 0 : void TwoSidedShape::visibility(Matrix<DComplex>& scale,
164 : const Matrix<Double>& uvw,
165 : const Vector<Double>& frequency) const {
166 0 : ComponentShape::visibility(scale, uvw, frequency);
167 0 : }
168 :
169 0 : Bool TwoSidedShape::isSymmetric() const {
170 0 : DebugAssert(ok(), AipsError);
171 0 : return true;
172 : }
173 :
174 308 : uInt TwoSidedShape::nParameters() const {
175 308 : DebugAssert(ok(), AipsError);
176 308 : return 3;
177 : }
178 :
179 154 : void TwoSidedShape::setParameters(const Vector<Double>& newParms) {
180 154 : DebugAssert(newParms.nelements() == nParameters(), AipsError);
181 154 : DebugAssert(newParms(0) >= newParms(1), AipsError);
182 154 : DebugAssert(abs(newParms(2)) <= C::_2pi, AipsError);
183 154 : setWidthInRad(newParms(0), newParms(1), newParms(2));
184 154 : DebugAssert(ok(), AipsError);
185 154 : }
186 :
187 126 : Vector<Double> TwoSidedShape::parameters() const {
188 126 : DebugAssert(ok(), AipsError);
189 126 : Vector<Double> compParms(3);
190 126 : compParms(0) = majorAxisInRad();
191 126 : compParms(1) = minorAxisInRad();
192 126 : compParms(2) = positionAngleInRad();
193 126 : return compParms;
194 0 : }
195 :
196 154 : void TwoSidedShape::setErrors(const Vector<Double>& newErrors) {
197 154 : DebugAssert(newErrors.nelements() == nParameters(), AipsError);
198 154 : DebugAssert(allGE(newErrors, 0.0), AipsError);
199 154 : const Unit rad("rad");
200 154 : itsMajErr.setValue(newErrors(0));
201 154 : itsMajErr.setUnit(rad);
202 154 : itsMinErr.setValue(newErrors(1));
203 154 : itsMinErr.setUnit(rad);
204 154 : itsPaErr.setValue(newErrors(2));
205 154 : itsPaErr.setUnit(rad);
206 154 : DebugAssert(ok(), AipsError);
207 154 : }
208 :
209 126 : Vector<Double> TwoSidedShape::errors() const {
210 126 : DebugAssert(ok(), AipsError);
211 126 : Vector<Double> compErrors(3);
212 126 : compErrors(0) = itsMajErr.getBaseValue();
213 126 : compErrors(1) = itsMinErr.getBaseValue();
214 126 : compErrors(2) = itsPaErr.getBaseValue();
215 126 : return compErrors;
216 0 : }
217 :
218 0 : Vector<Double> TwoSidedShape::optParameters() const {
219 0 : DebugAssert(ok(), AipsError);
220 0 : return Vector<Double>(0);
221 : }
222 :
223 0 : void TwoSidedShape::setOptParameters(const Vector<Double>& newOptParms){
224 0 : DebugAssert(ok(), AipsError);
225 : // squash compiler warning, maybe just get rid of DebugAssert statement
226 0 : if (newOptParms.empty()) {};
227 0 : }
228 :
229 6350 : Bool TwoSidedShape::fromRecord(String& errorMessage,
230 : const RecordInterface& record) {
231 6350 : if (!ComponentShape::fromRecord(errorMessage, record)) return false;
232 6350 : Quantum<Double> majorAxis, minorAxis, pa;
233 6350 : if (!fromAngQRecord(majorAxis, errorMessage, "majoraxis", record) ||
234 19050 : !fromAngQRecord(minorAxis, errorMessage, "minoraxis", record) ||
235 12700 : !fromAngQRecord(pa, errorMessage, "positionangle", record)) {
236 0 : errorMessage += "Shape not changed\n";
237 0 : return false;
238 : }
239 6350 : const Unit rad("rad");
240 6350 : const Double majorAxisInRad = majorAxis.getValue(rad);
241 6350 : const Double minorAxisInRad = minorAxis.getValue(rad);
242 : // // The near function is necessary for Intel processors (and doesn't hurt for
243 : // // other architectures) because of the extra precision that floating point
244 : // // variables have when returned in floating point registers. See
245 : // // http://aips2.nrao.edu/mail/aips2-lib/1101 for a discussion of this. The
246 : // // near function was added here and in the setMinorAxis function to fix
247 : // // defect AOCso00071
248 6350 : if (majorAxisInRad < minorAxisInRad &&
249 0 : !near(minorAxisInRad, minorAxisInRad, 2*C::dbl_epsilon)) {
250 0 : errorMessage += "The major axis cannot be smaller than the minor axis\n";
251 0 : return false;
252 : }
253 6350 : setWidth(majorAxis, minorAxis, pa);
254 6350 : if (!fromAngQRecord(majorAxis, errorMessage, "majoraxiserror", record) ||
255 19050 : !fromAngQRecord(minorAxis, errorMessage, "minoraxiserror", record) ||
256 12700 : !fromAngQRecord(pa, errorMessage, "positionangleerror", record)) {
257 0 : errorMessage += "Shape errors not changed\n";
258 0 : return false;
259 : }
260 6350 : setErrors(majorAxis, minorAxis, pa);
261 6350 : DebugAssert(ok(), AipsError);
262 6350 : return true;
263 6350 : }
264 :
265 274 : Bool TwoSidedShape::toRecord(String& errorMessage,
266 : RecordInterface& record) const {
267 274 : DebugAssert(ok(), AipsError);
268 274 : if (!ComponentShape::toRecord(errorMessage, record)) return false;
269 : {
270 274 : const QuantumHolder qHolder(majorAxis());
271 274 : Record qRecord;
272 274 : if (!qHolder.toRecord(errorMessage, qRecord)) {
273 0 : errorMessage += "Cannot convert the major axis to a record\n";
274 0 : return false;
275 : }
276 274 : record.defineRecord(RecordFieldId("majoraxis"), qRecord);
277 274 : }
278 : {
279 274 : const QuantumHolder qHolder(minorAxis());
280 274 : Record qRecord;
281 274 : if (!qHolder.toRecord(errorMessage, qRecord)) {
282 0 : errorMessage += "Cannot convert the minor axis to a record\n";
283 0 : return false;
284 : }
285 274 : record.defineRecord(RecordFieldId("minoraxis"), qRecord);
286 274 : }
287 : {
288 274 : const QuantumHolder qHolder(positionAngle());
289 274 : Record qRecord;
290 274 : if (!qHolder.toRecord(errorMessage, qRecord)) {
291 0 : errorMessage += "Cannot convert the position angle to a record\n";
292 0 : return false;
293 : }
294 274 : record.defineRecord(RecordFieldId("positionangle"), qRecord);
295 274 : }
296 : {
297 274 : const QuantumHolder qHolder(majorAxisError());
298 274 : Record qRecord;
299 274 : if (!qHolder.toRecord(errorMessage, qRecord)) {
300 0 : errorMessage += "Cannot convert the major axis error to a record\n";
301 0 : return false;
302 : }
303 274 : record.defineRecord(RecordFieldId("majoraxiserror"), qRecord);
304 274 : }
305 : {
306 274 : const QuantumHolder qHolder(minorAxisError());
307 274 : Record qRecord;
308 274 : if (!qHolder.toRecord(errorMessage, qRecord)) {
309 0 : errorMessage += "Cannot convert the minor axis error to a record\n";
310 0 : return false;
311 : }
312 274 : record.defineRecord(RecordFieldId("minoraxiserror"), qRecord);
313 274 : }
314 : {
315 274 : const QuantumHolder qHolder(positionAngleError());
316 274 : Record qRecord;
317 274 : if (!qHolder.toRecord(errorMessage, qRecord)) {
318 0 : errorMessage += "Cannot convert the position angle error to a record\n";
319 0 : return false;
320 : }
321 274 : record.defineRecord(RecordFieldId("positionangleerror"), qRecord);
322 274 : }
323 274 : return true;
324 : }
325 :
326 0 : Bool TwoSidedShape::convertUnit(String& errorMessage,
327 : const RecordInterface& record) {
328 0 : const Unit deg("deg");
329 : {
330 0 : const String fieldString("majoraxis");
331 0 : if (!record.isDefined(fieldString)) {
332 0 : errorMessage += "The 'majoraxis' field does not exist\n";
333 0 : return false;
334 : }
335 0 : const RecordFieldId field(fieldString);
336 0 : if (!((record.dataType(field) == TpString) &&
337 0 : (record.shape(field) == IPosition(1,1)))) {
338 0 : errorMessage += "The 'majoraxis' field must be a string\n";
339 0 : errorMessage += "(but not a vector of strings)\n";
340 0 : return false;
341 : }
342 0 : const Unit unit = Unit(record.asString(field));
343 0 : if (unit != deg) {
344 : errorMessage +=
345 0 : "Cannot convert the major axis width to a non angular unit";
346 0 : return false;
347 : }
348 0 : itsMajUnit = unit;
349 0 : }
350 : {
351 0 : const String fieldString("minoraxis");
352 0 : if (!record.isDefined(fieldString)) {
353 0 : errorMessage += "The 'minoraxis' field does not exist\n";
354 0 : return false;
355 : }
356 0 : const RecordFieldId field(fieldString);
357 0 : if (!((record.dataType(field) == TpString) &&
358 0 : (record.shape(field) == IPosition(1,1)))) {
359 0 : errorMessage += "The 'minoraxis' field must be a string\n";
360 0 : errorMessage += "(but not a vector of strings)\n";
361 0 : return false;
362 : }
363 0 : const Unit unit = Unit(record.asString(field));
364 0 : if (unit != deg) {
365 : errorMessage +=
366 0 : "Cannot convert the minor axis width to a non angular unit";
367 0 : return false;
368 : }
369 0 : itsMinUnit = unit;
370 0 : }
371 : {
372 0 : const String fieldString("positionangle");
373 0 : if (!record.isDefined(fieldString)) {
374 0 : errorMessage += "The 'positionangle' field does not exist\n";
375 0 : return false;
376 : }
377 0 : const RecordFieldId field(fieldString);
378 0 : if (!((record.dataType(field) == TpString) &&
379 0 : (record.shape(field) == IPosition(1,1)))) {
380 0 : errorMessage += "The 'positionangle' field must be a string\n";
381 0 : errorMessage += "(but not a vector of strings)\n";
382 0 : return false;
383 : }
384 0 : const Unit unit = Unit(record.asString(field));
385 0 : if (unit != deg) {
386 : errorMessage +=
387 0 : "Cannot convert the position angle to a non angular unit";
388 0 : return false;
389 : }
390 0 : itsPaUnit = unit;
391 0 : }
392 0 : DebugAssert(ok(), AipsError);
393 0 : return true;
394 0 : }
395 :
396 694106 : Bool TwoSidedShape::ok() const {
397 : // The LogIO class is only constructed if an error is detected for
398 : // performance reasons. Both function static and file static variables
399 : // where considered and rejected for this purpose.
400 694106 : if (!ComponentShape::ok()) return false;
401 694106 : const Unit deg("deg");
402 694106 : if (itsMajUnit != deg) {
403 0 : LogIO logErr(LogOrigin("TwoSidedCompRep", "ok()"));
404 : logErr << LogIO::SEVERE << "The major axis does not have angular units."
405 0 : << LogIO::POST;
406 0 : return false;
407 0 : }
408 694106 : if (itsMinUnit != deg) {
409 0 : LogIO logErr(LogOrigin("TwoSidedCompRep", "ok()"));
410 : logErr << LogIO::SEVERE << "The minor axis does not have angular units."
411 0 : << LogIO::POST;
412 0 : return false;
413 0 : }
414 694106 : if (itsPaUnit != deg) {
415 0 : LogIO logErr(LogOrigin("TwoSidedCompRep", "ok()"));
416 : logErr << LogIO::SEVERE <<"The position angle does not have angular units."
417 0 : << LogIO::POST;
418 0 : return false;
419 0 : }
420 694106 : return true;
421 694106 : }
422 :
423 6504 : TwoSidedShape::TwoSidedShape()
424 : :ComponentShape(),
425 6504 : itsMajUnit("arcmin"),
426 6504 : itsMinUnit("arcmin"),
427 6504 : itsPaUnit("deg"),
428 6504 : itsMajErr(0, "arcmin"),
429 6504 : itsMinErr(0, "arcmin"),
430 13008 : itsPaErr(0, "deg")
431 : {
432 6504 : DebugAssert(ok(), AipsError);
433 6504 : }
434 :
435 0 : TwoSidedShape::TwoSidedShape(const MDirection& direction,
436 : const Unit& majorAxisUnit,
437 : const Unit& minorAxisUnit,
438 0 : const Unit& paUnit)
439 : :ComponentShape(direction),
440 0 : itsMajUnit(majorAxisUnit),
441 0 : itsMinUnit(minorAxisUnit),
442 0 : itsPaUnit(paUnit),
443 0 : itsMajErr(0, "arcmin"),
444 0 : itsMinErr(0, "arcmin"),
445 0 : itsPaErr(0, "deg")
446 : {
447 0 : }
448 :
449 7878 : TwoSidedShape::TwoSidedShape(const TwoSidedShape& other)
450 : :ComponentShape(other),
451 7878 : itsMajUnit(other.itsMajUnit),
452 7878 : itsMinUnit(other.itsMinUnit),
453 7878 : itsPaUnit(other.itsPaUnit),
454 7878 : itsMajErr(other.itsMajErr),
455 7878 : itsMinErr(other.itsMinErr),
456 15756 : itsPaErr(other.itsPaErr)
457 : {
458 7878 : DebugAssert(ok(), AipsError);
459 7878 : }
460 :
461 :
462 0 : Vector<Double> TwoSidedShape::toPixel (const DirectionCoordinate& dirCoord) const
463 : //
464 : // pars(0) = long cen abs pix
465 : // pars(1) = lat cen abs pix
466 : // pars(2) = major pix
467 : // pars(3) = minor pix
468 : // pars(4) = pa radians; pos +x (long) -> +y (lat)
469 : //
470 : {
471 0 : LogIO os(LogOrigin("TwoSidedShape", "toPixel"));
472 0 : Vector<Double> parameters(5);
473 :
474 : // Do locations
475 :
476 0 : Vector<Double> pixelCen = ComponentShape::toPixel (dirCoord);
477 0 : parameters(0) = pixelCen(0);
478 0 : parameters(1) = pixelCen(1);
479 :
480 : // Now convert the tip of the major axis to x/y pixel coordinates
481 :
482 0 : const MDirection dirRef = refDirection();
483 0 : Quantum<Double> majorWorld = majorAxis();
484 0 : Quantum<Double> paMajor = positionAngle();
485 0 : majorWorld.scale(0.5);
486 0 : Vector<Double> majorCart = widthToCartesian (majorWorld, paMajor, dirRef, dirCoord, pixelCen);
487 :
488 : // Position angle of major axis.
489 : // atan2 gives pos +x (long) -> +y (lat). put in range +/- pi
490 :
491 0 : MVAngle pa(atan2(majorCart(1), majorCart(0)));
492 0 : pa();
493 :
494 : // I cannot just add 90deg to the world position angle. It is 90deg in the
495 : // pixel coordinate frame, not the world frame. So I have to work
496 : // my way along the minor axis in pixel coordinates and locate
497 : // the tip of the minor axis iteratively. The algorithm
498 : // below could be much smarter/faster with a binary search.
499 :
500 0 : Quantum<Double> minorWorld = minorAxis();
501 0 : Quantum<Double> paMinor = paMajor + Quantum<Double>(C::pi/2.0, Unit("rad"));
502 0 : minorWorld.scale(0.5);
503 : //
504 0 : Double dX = sin(pa.radian());
505 0 : Double dY = cos(pa.radian());
506 : //
507 0 : Vector<Double> posPix = pixelCen.copy();
508 0 : MDirection posWorld;
509 0 : MVDirection mvdRef = dirRef.getValue();
510 0 : Vector<Double> prevPosPix(2);
511 : //
512 0 : Double minorWorldRad = minorWorld.getValue(Unit("rad"));
513 0 : Double sep = 0.0;
514 0 : Double prevSep = 0.0;
515 0 : Bool more = true;
516 0 : while (more) {
517 0 : dirCoord.toWorld(posWorld, posPix);
518 0 : MVDirection mvd = posWorld.getValue();
519 0 : sep = mvdRef.separation(mvd);
520 0 : if (sep > minorWorldRad) break;
521 : //
522 0 : prevPosPix = posPix;
523 0 : prevSep = sep;
524 : //
525 0 : posPix(0) += dX;
526 0 : posPix(1) += dY;
527 0 : }
528 0 : Double frac = (minorWorldRad - prevSep) / (sep - prevSep);
529 0 : Double fracX = dX * frac;
530 0 : Double fracY = dY * frac;
531 : //
532 0 : Vector<Double> minorCart(2);
533 0 : minorCart(0) = prevPosPix(0) + fracX - pixelCen(0);
534 0 : minorCart(1) = prevPosPix(1) + fracY - pixelCen(1);
535 : //
536 0 : Double tmp1 = 2.0 * hypot(majorCart(0), majorCart(1));
537 0 : Double tmp2 = 2.0 * hypot(minorCart(0), minorCart(1));
538 : //
539 0 : parameters(2) = max(tmp1,tmp2);
540 0 : parameters(3) = min(tmp1,tmp2);
541 0 : parameters(4) = pa.radian();
542 : //
543 0 : return parameters;
544 0 : }
545 :
546 0 : Bool TwoSidedShape::fromPixel (const Vector<Double>& parameters,
547 : const DirectionCoordinate& dirCoord)
548 : //
549 : // pars(0) = long cen abs pix
550 : // pars(1) = lat cen abs pix
551 : // pars(2) = major pix
552 : // pars(3) = minor pix
553 : // pars(4) = pa radians; pos +x (long) -> +y (lat)
554 : //
555 : {
556 : // Direction first
557 :
558 0 : Vector<Double> pixelCen(2);
559 0 : pixelCen(0) = parameters(0);
560 0 : pixelCen(1) = parameters(1);
561 0 : ComponentShape::fromPixel (pixelCen, dirCoord);
562 : // Shape. First put x/y p.a. into +y -> -x system
563 :
564 0 : Double pa0 = parameters(4) - C::pi_2;
565 0 : MDirection tipMajor = directionFromCartesian (parameters(2), pa0, dirCoord, pixelCen);
566 : //
567 0 : pa0 += C::pi_2; // minor axis position angle
568 0 : MDirection tipMinor = directionFromCartesian (parameters(3), pa0, dirCoord, pixelCen);
569 :
570 : // Find tip directions
571 0 : const MDirection& directionRef = refDirection();
572 0 : MVDirection mvdRef = directionRef.getValue();
573 0 : MVDirection mvdMajor = tipMajor.getValue();
574 0 : MVDirection mvdMinor = tipMinor.getValue();
575 :
576 : // Separations
577 :
578 0 : Double tmp1 = 2 * mvdRef.separation(mvdMajor) * 3600 * 180.0 / C::pi;
579 0 : Double tmp2 = 2 * mvdRef.separation(mvdMinor) * 3600 * 180.0 / C::pi;
580 :
581 0 : Quantity majorAxis(max(tmp1,tmp2), Unit("arcsec"));
582 0 : Quantity minorAxis(min(tmp1,tmp2), Unit("arcsec"));
583 0 : Bool flipped = tmp2 > tmp1;
584 0 : Quantity pa;
585 0 : if (!flipped) {
586 0 : pa = mvdRef.positionAngle(mvdMajor, Unit("deg"));
587 : } else {
588 0 : pa = mvdRef.positionAngle(mvdMinor, Unit("deg"));
589 : }
590 0 : setWidth (majorAxis, minorAxis, pa);
591 0 : return flipped;
592 0 : }
593 :
594 0 : Vector<Double> TwoSidedShape::widthToCartesian (const Quantum<Double>& width,
595 : const Quantum<Double>& pa,
596 : const MDirection& dirRef,
597 : const DirectionCoordinate& dirCoord,
598 : const Vector<Double>& pixelCen) const
599 : {
600 :
601 : // Find MDirection of tip of axis
602 :
603 0 : MDirection dirTip = dirRef;
604 0 : dirTip.shiftAngle(width, pa);
605 :
606 : // Convert to pixel
607 :
608 0 : Vector<Double> pixelTip(2);
609 0 : if (!dirCoord.toPixel(pixelTip, dirTip)) {
610 0 : LogIO os(LogOrigin("TwoSidedShape", "widthToCartesian"));
611 : os << "DirectionCoordinate conversion to pixel failed because "
612 0 : << dirCoord.errorMessage() << LogIO::EXCEPTION;
613 0 : }
614 :
615 : // Find offset cartesian components
616 :
617 0 : Vector<Double> cart(2);
618 0 : cart(0) = pixelTip(0) - pixelCen(0);
619 0 : cart(1) = pixelTip(1) - pixelCen(1);
620 0 : return cart;
621 0 : }
622 :
623 0 : MDirection TwoSidedShape::directionFromCartesian (Double width, Double pa,
624 : const DirectionCoordinate& dirCoord,
625 : const Vector<Double>& pixelCen) const
626 : {
627 :
628 : // Now find tips of major and minor axes in pixel coordinates
629 : // and convert to world
630 :
631 0 : Double z = width / 2.0;
632 0 : Double x = -z * sin(pa);
633 0 : Double y = z * cos(pa);
634 0 : MDirection dir;
635 0 : Vector<Double> pixelTip(2);
636 0 : pixelTip(0) = pixelCen(0) + x;
637 0 : pixelTip(1) = pixelCen(1) + y;
638 0 : ThrowIf(
639 : ! dirCoord.toWorld(dir, pixelTip),
640 : "DirectionCoordinate conversion failed because "
641 : + dirCoord.errorMessage()
642 : );
643 0 : return dir;
644 0 : }
645 :
646 0 : String TwoSidedShape::sizeToString(
647 : Quantity major, Quantity minor, Quantity posangle,
648 : Bool includeUncertainties, Quantity majorErr,
649 : Quantity minorErr, Quantity posanErr
650 : ) {
651 : // Inputs all as angle quanta
652 0 : Vector<String> angUnits(5);
653 0 : angUnits[0] = "deg";
654 0 : angUnits[1] = "arcmin";
655 0 : angUnits[2] = "arcsec";
656 0 : angUnits[3] = "marcsec";
657 0 : angUnits[4] = "uarcsec";
658 : // First force position angle to be between 0 and 180 deg
659 0 : if(posangle.getValue() < 0) {
660 0 : posangle += Quantity(180, "deg");
661 : }
662 :
663 0 : String prefUnits;
664 0 : Quantity vmax(max(fabs(major.getValue("arcsec")), fabs(minor.getValue("arcsec"))), "arcsec");
665 :
666 0 : for (uInt i=0; i<angUnits.size(); i++) {
667 0 : prefUnits = angUnits[i];
668 0 : if(vmax.getValue(prefUnits) > 1) {
669 0 : break;
670 : }
671 : }
672 0 : major.convert(prefUnits);
673 0 : minor.convert(prefUnits);
674 0 : majorErr.convert(prefUnits);
675 0 : minorErr.convert(prefUnits);
676 :
677 0 : Double vmaj = major.getValue();
678 0 : Double vmin = minor.getValue();
679 :
680 : // Formatting as "value +/- err" for symmetric errors
681 :
682 0 : Double dmaj = majorErr.getValue();
683 0 : Double dmin = minorErr.getValue();
684 : // position angle is always in degrees cuz users like that
685 0 : Double pa = posangle.getValue("deg");
686 0 : Double dpa = posanErr.getValue("deg");
687 :
688 0 : Vector<Double> majVec(2), minVec(2), paVec(2);
689 0 : majVec[0] = vmaj;
690 0 : majVec[1] = dmaj;
691 0 : minVec[0] = vmin;
692 0 : minVec[1] = dmin;
693 0 : paVec[0] = pa;
694 0 : paVec[1] = dpa;
695 0 : uInt precision1 = precisionForValueErrorPairs(majVec, minVec);
696 0 : uInt precision2 = precisionForValueErrorPairs(paVec, Vector<Double>(0));
697 :
698 0 : ostringstream summary;
699 0 : summary << std::fixed << std::setprecision(precision1);
700 0 : summary << " --- major axis FWHM: " << major.getValue();
701 0 : if (includeUncertainties) {
702 0 : if (majorErr.getValue() == 0) {
703 0 : summary << " " << prefUnits << " (fixed)" << endl;
704 : }
705 : else {
706 0 : summary << " +/- " << majorErr.getValue()
707 0 : << " " << prefUnits << endl;
708 : }
709 : }
710 : else {
711 0 : summary << " " << prefUnits << endl;
712 : }
713 0 : summary << " --- minor axis FWHM: " << minor.getValue();
714 0 : if (includeUncertainties) {
715 0 : if (minorErr.getValue() == 0) {
716 0 : summary << " " << prefUnits << " (fixed)" << endl;
717 : }
718 : else {
719 0 : summary << " +/- " << minorErr.getValue()
720 0 : << " " << prefUnits << endl;
721 : }
722 : }
723 : else {
724 0 : summary << " " << prefUnits << endl;
725 : }
726 0 : summary << std::setprecision(precision2);
727 0 : summary << " --- position angle: " << pa;
728 0 : if (includeUncertainties) {
729 0 : if (dpa == 0) {
730 0 : summary << "deg (fixed)" << endl;
731 : }
732 : else {
733 0 : summary << " +/- " << dpa << " deg" << endl;
734 : }
735 : }
736 : else {
737 0 : summary << " deg" << endl;
738 : }
739 0 : return summary.str();
740 0 : }
741 :
742 :
743 : // Local Variables:
744 : // compile-command: "gmake OPTLIB=1 TwoSidedShape"
745 : // End:
746 :
747 :
748 : } //# NAMESPACE CASA - END
749 :
|
