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 34951 : TwoSidedShape::~TwoSidedShape() {
52 34951 : DebugAssert(ok(), AipsError);
53 34951 : }
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 15534 : void TwoSidedShape::setWidth(const Quantity& majorAxis,
70 : const Quantity& minorAxis,
71 : const Quantity& positionAngle) {
72 15534 : itsMajUnit = majorAxis.getFullUnit();
73 15534 : itsMinUnit = minorAxis.getFullUnit();
74 15534 : itsPaUnit = positionAngle.getFullUnit();
75 15534 : const Unit rad("rad");
76 15534 : setWidthInRad(majorAxis.getValue(rad), minorAxis.getValue(rad),
77 : positionAngle.getValue(rad));
78 15534 : DebugAssert(ok(), AipsError);
79 15534 : }
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 3929 : Quantum<Double> TwoSidedShape::majorAxis() const {
93 3929 : Quantum<Double> retVal(majorAxisInRad(), Unit("rad"));
94 3929 : retVal.convert(itsMajUnit);
95 3929 : return retVal;
96 0 : }
97 :
98 3929 : Quantum<Double> TwoSidedShape::minorAxis() const {
99 3929 : Quantum<Double> retVal(minorAxisInRad(), Unit("rad"));
100 3929 : retVal.convert(itsMinUnit);
101 3929 : return retVal;
102 0 : }
103 :
104 3871 : Quantum<Double> TwoSidedShape::positionAngle() const {
105 3871 : Quantum<Double> retVal(positionAngleInRad(), Unit("rad"));
106 3871 : retVal.convert(itsPaUnit);
107 3871 : return retVal;
108 0 : }
109 :
110 0 : Double TwoSidedShape::axialRatio() const {
111 0 : return minorAxisInRad()/majorAxisInRad();
112 : }
113 :
114 15398 : void TwoSidedShape::setErrors(const Quantum<Double>& majorAxisError,
115 : const Quantum<Double>& minorAxisError,
116 : const Quantum<Double>& positionAngleError) {
117 15398 : if (ComponentShape::badError(majorAxisError) ||
118 30796 : ComponentShape::badError(minorAxisError) ||
119 15398 : 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 15398 : itsMajErr = majorAxisError;
125 15398 : itsMinErr = minorAxisError;
126 15398 : itsPaErr = positionAngleError;
127 15398 : }
128 :
129 2214 : const Quantum<Double>& TwoSidedShape::majorAxisError() const {
130 2214 : return itsMajErr;
131 : }
132 :
133 2214 : const Quantum<Double>& TwoSidedShape::minorAxisError() const {
134 2214 : return itsMinErr;
135 : }
136 :
137 2214 : const Quantum<Double>& TwoSidedShape::positionAngleError() const {
138 2214 : 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 546 : uInt TwoSidedShape::nParameters() const {
175 546 : DebugAssert(ok(), AipsError);
176 546 : return 3;
177 : }
178 :
179 268 : void TwoSidedShape::setParameters(const Vector<Double>& newParms) {
180 268 : DebugAssert(newParms.nelements() == nParameters(), AipsError);
181 268 : DebugAssert(newParms(0) >= newParms(1), AipsError);
182 268 : DebugAssert(abs(newParms(2)) <= C::_2pi, AipsError);
183 268 : setWidthInRad(newParms(0), newParms(1), newParms(2));
184 268 : DebugAssert(ok(), AipsError);
185 268 : }
186 :
187 324 : Vector<Double> TwoSidedShape::parameters() const {
188 324 : DebugAssert(ok(), AipsError);
189 324 : Vector<Double> compParms(3);
190 324 : compParms(0) = majorAxisInRad();
191 324 : compParms(1) = minorAxisInRad();
192 324 : compParms(2) = positionAngleInRad();
193 324 : return compParms;
194 0 : }
195 :
196 246 : void TwoSidedShape::setErrors(const Vector<Double>& newErrors) {
197 246 : DebugAssert(newErrors.nelements() == nParameters(), AipsError);
198 246 : DebugAssert(allGE(newErrors, 0.0), AipsError);
199 246 : const Unit rad("rad");
200 246 : itsMajErr.setValue(newErrors(0));
201 246 : itsMajErr.setUnit(rad);
202 246 : itsMinErr.setValue(newErrors(1));
203 246 : itsMinErr.setUnit(rad);
204 246 : itsPaErr.setValue(newErrors(2));
205 246 : itsPaErr.setUnit(rad);
206 246 : DebugAssert(ok(), AipsError);
207 246 : }
208 :
209 292 : Vector<Double> TwoSidedShape::errors() const {
210 292 : DebugAssert(ok(), AipsError);
211 292 : Vector<Double> compErrors(3);
212 292 : compErrors(0) = itsMajErr.getBaseValue();
213 292 : compErrors(1) = itsMinErr.getBaseValue();
214 292 : compErrors(2) = itsPaErr.getBaseValue();
215 292 : 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 14854 : Bool TwoSidedShape::fromRecord(String& errorMessage,
230 : const RecordInterface& record) {
231 14854 : if (!ComponentShape::fromRecord(errorMessage, record)) return false;
232 14854 : Quantum<Double> majorAxis, minorAxis, pa;
233 14854 : if (!fromAngQRecord(majorAxis, errorMessage, "majoraxis", record) ||
234 44562 : !fromAngQRecord(minorAxis, errorMessage, "minoraxis", record) ||
235 29708 : !fromAngQRecord(pa, errorMessage, "positionangle", record)) {
236 0 : errorMessage += "Shape not changed\n";
237 0 : return false;
238 : }
239 14854 : const Unit rad("rad");
240 14854 : const Double majorAxisInRad = majorAxis.getValue(rad);
241 14854 : 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 14854 : 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 14854 : setWidth(majorAxis, minorAxis, pa);
254 14854 : if (!fromAngQRecord(majorAxis, errorMessage, "majoraxiserror", record) ||
255 44562 : !fromAngQRecord(minorAxis, errorMessage, "minoraxiserror", record) ||
256 29708 : !fromAngQRecord(pa, errorMessage, "positionangleerror", record)) {
257 0 : errorMessage += "Shape errors not changed\n";
258 0 : return false;
259 : }
260 14854 : setErrors(majorAxis, minorAxis, pa);
261 14854 : DebugAssert(ok(), AipsError);
262 14854 : return true;
263 14854 : }
264 :
265 1805 : Bool TwoSidedShape::toRecord(String& errorMessage,
266 : RecordInterface& record) const {
267 1805 : DebugAssert(ok(), AipsError);
268 1805 : if (!ComponentShape::toRecord(errorMessage, record)) return false;
269 : {
270 1805 : const QuantumHolder qHolder(majorAxis());
271 1805 : Record qRecord;
272 1805 : if (!qHolder.toRecord(errorMessage, qRecord)) {
273 0 : errorMessage += "Cannot convert the major axis to a record\n";
274 0 : return false;
275 : }
276 1805 : record.defineRecord(RecordFieldId("majoraxis"), qRecord);
277 1805 : }
278 : {
279 1805 : const QuantumHolder qHolder(minorAxis());
280 1805 : Record qRecord;
281 1805 : if (!qHolder.toRecord(errorMessage, qRecord)) {
282 0 : errorMessage += "Cannot convert the minor axis to a record\n";
283 0 : return false;
284 : }
285 1805 : record.defineRecord(RecordFieldId("minoraxis"), qRecord);
286 1805 : }
287 : {
288 1805 : const QuantumHolder qHolder(positionAngle());
289 1805 : Record qRecord;
290 1805 : if (!qHolder.toRecord(errorMessage, qRecord)) {
291 0 : errorMessage += "Cannot convert the position angle to a record\n";
292 0 : return false;
293 : }
294 1805 : record.defineRecord(RecordFieldId("positionangle"), qRecord);
295 1805 : }
296 : {
297 1805 : const QuantumHolder qHolder(majorAxisError());
298 1805 : Record qRecord;
299 1805 : if (!qHolder.toRecord(errorMessage, qRecord)) {
300 0 : errorMessage += "Cannot convert the major axis error to a record\n";
301 0 : return false;
302 : }
303 1805 : record.defineRecord(RecordFieldId("majoraxiserror"), qRecord);
304 1805 : }
305 : {
306 1805 : const QuantumHolder qHolder(minorAxisError());
307 1805 : Record qRecord;
308 1805 : if (!qHolder.toRecord(errorMessage, qRecord)) {
309 0 : errorMessage += "Cannot convert the minor axis error to a record\n";
310 0 : return false;
311 : }
312 1805 : record.defineRecord(RecordFieldId("minoraxiserror"), qRecord);
313 1805 : }
314 : {
315 1805 : const QuantumHolder qHolder(positionAngleError());
316 1805 : Record qRecord;
317 1805 : if (!qHolder.toRecord(errorMessage, qRecord)) {
318 0 : errorMessage += "Cannot convert the position angle error to a record\n";
319 0 : return false;
320 : }
321 1805 : record.defineRecord(RecordFieldId("positionangleerror"), qRecord);
322 1805 : }
323 1805 : 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 20686581 : 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 20686581 : if (!ComponentShape::ok()) return false;
401 20686581 : const Unit deg("deg");
402 20686581 : 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 20686581 : 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 20686581 : 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 20686581 : return true;
421 20686581 : }
422 :
423 15642 : TwoSidedShape::TwoSidedShape()
424 : :ComponentShape(),
425 15642 : itsMajUnit("arcmin"),
426 15642 : itsMinUnit("arcmin"),
427 15642 : itsPaUnit("deg"),
428 15642 : itsMajErr(0, "arcmin"),
429 15642 : itsMinErr(0, "arcmin"),
430 31284 : itsPaErr(0, "deg")
431 : {
432 15642 : DebugAssert(ok(), AipsError);
433 15642 : }
434 :
435 723 : TwoSidedShape::TwoSidedShape(const MDirection& direction,
436 : const Unit& majorAxisUnit,
437 : const Unit& minorAxisUnit,
438 723 : const Unit& paUnit)
439 : :ComponentShape(direction),
440 723 : itsMajUnit(majorAxisUnit),
441 723 : itsMinUnit(minorAxisUnit),
442 723 : itsPaUnit(paUnit),
443 723 : itsMajErr(0, "arcmin"),
444 723 : itsMinErr(0, "arcmin"),
445 1446 : itsPaErr(0, "deg")
446 : {
447 723 : }
448 :
449 18586 : TwoSidedShape::TwoSidedShape(const TwoSidedShape& other)
450 : :ComponentShape(other),
451 18586 : itsMajUnit(other.itsMajUnit),
452 18586 : itsMinUnit(other.itsMinUnit),
453 18586 : itsPaUnit(other.itsPaUnit),
454 18586 : itsMajErr(other.itsMajErr),
455 18586 : itsMinErr(other.itsMinErr),
456 37172 : itsPaErr(other.itsPaErr)
457 : {
458 18586 : DebugAssert(ok(), AipsError);
459 18586 : }
460 :
461 :
462 738 : 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 1476 : LogIO os(LogOrigin("TwoSidedShape", "toPixel"));
472 738 : Vector<Double> parameters(5);
473 :
474 : // Do locations
475 :
476 738 : Vector<Double> pixelCen = ComponentShape::toPixel (dirCoord);
477 738 : parameters(0) = pixelCen(0);
478 738 : parameters(1) = pixelCen(1);
479 :
480 : // Now convert the tip of the major axis to x/y pixel coordinates
481 :
482 738 : const MDirection dirRef = refDirection();
483 738 : Quantum<Double> majorWorld = majorAxis();
484 738 : Quantum<Double> paMajor = positionAngle();
485 738 : majorWorld.scale(0.5);
486 738 : 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 738 : MVAngle pa(atan2(majorCart(1), majorCart(0)));
492 738 : 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 738 : Quantum<Double> minorWorld = minorAxis();
501 1476 : Quantum<Double> paMinor = paMajor + Quantum<Double>(C::pi/2.0, Unit("rad"));
502 738 : minorWorld.scale(0.5);
503 : //
504 738 : Double dX = sin(pa.radian());
505 738 : Double dY = cos(pa.radian());
506 : //
507 738 : Vector<Double> posPix = pixelCen.copy();
508 738 : MDirection posWorld;
509 738 : MVDirection mvdRef = dirRef.getValue();
510 738 : Vector<Double> prevPosPix(2);
511 : //
512 738 : Double minorWorldRad = minorWorld.getValue(Unit("rad"));
513 738 : Double sep = 0.0;
514 738 : Double prevSep = 0.0;
515 738 : Bool more = true;
516 4135 : while (more) {
517 4135 : dirCoord.toWorld(posWorld, posPix);
518 4135 : MVDirection mvd = posWorld.getValue();
519 4135 : sep = mvdRef.separation(mvd);
520 4135 : if (sep > minorWorldRad) break;
521 : //
522 3397 : prevPosPix = posPix;
523 3397 : prevSep = sep;
524 : //
525 3397 : posPix(0) += dX;
526 3397 : posPix(1) += dY;
527 4135 : }
528 738 : Double frac = (minorWorldRad - prevSep) / (sep - prevSep);
529 738 : Double fracX = dX * frac;
530 738 : Double fracY = dY * frac;
531 : //
532 738 : Vector<Double> minorCart(2);
533 738 : minorCart(0) = prevPosPix(0) + fracX - pixelCen(0);
534 738 : minorCart(1) = prevPosPix(1) + fracY - pixelCen(1);
535 : //
536 738 : Double tmp1 = 2.0 * hypot(majorCart(0), majorCart(1));
537 738 : Double tmp2 = 2.0 * hypot(minorCart(0), minorCart(1));
538 : //
539 738 : parameters(2) = max(tmp1,tmp2);
540 738 : parameters(3) = min(tmp1,tmp2);
541 738 : parameters(4) = pa.radian();
542 : //
543 1476 : return parameters;
544 738 : }
545 :
546 526 : 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 526 : Vector<Double> pixelCen(2);
559 526 : pixelCen(0) = parameters(0);
560 526 : pixelCen(1) = parameters(1);
561 526 : ComponentShape::fromPixel (pixelCen, dirCoord);
562 : // Shape. First put x/y p.a. into +y -> -x system
563 :
564 526 : Double pa0 = parameters(4) - C::pi_2;
565 526 : MDirection tipMajor = directionFromCartesian (parameters(2), pa0, dirCoord, pixelCen);
566 : //
567 526 : pa0 += C::pi_2; // minor axis position angle
568 526 : MDirection tipMinor = directionFromCartesian (parameters(3), pa0, dirCoord, pixelCen);
569 :
570 : // Find tip directions
571 526 : const MDirection& directionRef = refDirection();
572 526 : MVDirection mvdRef = directionRef.getValue();
573 526 : MVDirection mvdMajor = tipMajor.getValue();
574 526 : MVDirection mvdMinor = tipMinor.getValue();
575 :
576 : // Separations
577 :
578 526 : Double tmp1 = 2 * mvdRef.separation(mvdMajor) * 3600 * 180.0 / C::pi;
579 526 : Double tmp2 = 2 * mvdRef.separation(mvdMinor) * 3600 * 180.0 / C::pi;
580 :
581 526 : Quantity majorAxis(max(tmp1,tmp2), Unit("arcsec"));
582 526 : Quantity minorAxis(min(tmp1,tmp2), Unit("arcsec"));
583 526 : Bool flipped = tmp2 > tmp1;
584 526 : Quantity pa;
585 526 : if (!flipped) {
586 526 : pa = mvdRef.positionAngle(mvdMajor, Unit("deg"));
587 : } else {
588 0 : pa = mvdRef.positionAngle(mvdMinor, Unit("deg"));
589 : }
590 526 : setWidth (majorAxis, minorAxis, pa);
591 526 : return flipped;
592 526 : }
593 :
594 738 : 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 738 : MDirection dirTip = dirRef;
604 738 : dirTip.shiftAngle(width, pa);
605 :
606 : // Convert to pixel
607 :
608 738 : Vector<Double> pixelTip(2);
609 738 : 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 738 : Vector<Double> cart(2);
618 738 : cart(0) = pixelTip(0) - pixelCen(0);
619 738 : cart(1) = pixelTip(1) - pixelCen(1);
620 1476 : return cart;
621 738 : }
622 :
623 1052 : 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 1052 : Double z = width / 2.0;
632 1052 : Double x = -z * sin(pa);
633 1052 : Double y = z * cos(pa);
634 1052 : MDirection dir;
635 1052 : Vector<Double> pixelTip(2);
636 1052 : pixelTip(0) = pixelCen(0) + x;
637 1052 : pixelTip(1) = pixelCen(1) + y;
638 1052 : ThrowIf(
639 : ! dirCoord.toWorld(dir, pixelTip),
640 : "DirectionCoordinate conversion failed because "
641 : + dirCoord.errorMessage()
642 : );
643 2104 : return dir;
644 1052 : }
645 :
646 541 : 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 541 : Vector<String> angUnits(5);
653 541 : angUnits[0] = "deg";
654 541 : angUnits[1] = "arcmin";
655 541 : angUnits[2] = "arcsec";
656 541 : angUnits[3] = "marcsec";
657 541 : angUnits[4] = "uarcsec";
658 : // First force position angle to be between 0 and 180 deg
659 541 : if(posangle.getValue() < 0) {
660 32 : posangle += Quantity(180, "deg");
661 : }
662 :
663 541 : String prefUnits;
664 541 : Quantity vmax(max(fabs(major.getValue("arcsec")), fabs(minor.getValue("arcsec"))), "arcsec");
665 :
666 1480 : for (uInt i=0; i<angUnits.size(); i++) {
667 1480 : prefUnits = angUnits[i];
668 1480 : if(vmax.getValue(prefUnits) > 1) {
669 541 : break;
670 : }
671 : }
672 541 : major.convert(prefUnits);
673 541 : minor.convert(prefUnits);
674 541 : majorErr.convert(prefUnits);
675 541 : minorErr.convert(prefUnits);
676 :
677 541 : Double vmaj = major.getValue();
678 541 : Double vmin = minor.getValue();
679 :
680 : // Formatting as "value +/- err" for symmetric errors
681 :
682 541 : Double dmaj = majorErr.getValue();
683 541 : Double dmin = minorErr.getValue();
684 : // position angle is always in degrees cuz users like that
685 541 : Double pa = posangle.getValue("deg");
686 541 : Double dpa = posanErr.getValue("deg");
687 :
688 541 : Vector<Double> majVec(2), minVec(2), paVec(2);
689 541 : majVec[0] = vmaj;
690 541 : majVec[1] = dmaj;
691 541 : minVec[0] = vmin;
692 541 : minVec[1] = dmin;
693 541 : paVec[0] = pa;
694 541 : paVec[1] = dpa;
695 541 : uInt precision1 = precisionForValueErrorPairs(majVec, minVec);
696 541 : uInt precision2 = precisionForValueErrorPairs(paVec, Vector<Double>(0));
697 :
698 541 : ostringstream summary;
699 541 : summary << std::fixed << std::setprecision(precision1);
700 541 : summary << " --- major axis FWHM: " << major.getValue();
701 541 : if (includeUncertainties) {
702 541 : if (majorErr.getValue() == 0) {
703 0 : summary << " " << prefUnits << " (fixed)" << endl;
704 : }
705 : else {
706 541 : summary << " +/- " << majorErr.getValue()
707 541 : << " " << prefUnits << endl;
708 : }
709 : }
710 : else {
711 0 : summary << " " << prefUnits << endl;
712 : }
713 541 : summary << " --- minor axis FWHM: " << minor.getValue();
714 541 : if (includeUncertainties) {
715 541 : if (minorErr.getValue() == 0) {
716 0 : summary << " " << prefUnits << " (fixed)" << endl;
717 : }
718 : else {
719 541 : summary << " +/- " << minorErr.getValue()
720 541 : << " " << prefUnits << endl;
721 : }
722 : }
723 : else {
724 0 : summary << " " << prefUnits << endl;
725 : }
726 541 : summary << std::setprecision(precision2);
727 541 : summary << " --- position angle: " << pa;
728 541 : if (includeUncertainties) {
729 541 : if (dpa == 0) {
730 0 : summary << "deg (fixed)" << endl;
731 : }
732 : else {
733 541 : summary << " +/- " << dpa << " deg" << endl;
734 : }
735 : }
736 : else {
737 0 : summary << " deg" << endl;
738 : }
739 1623 : return summary.str();
740 541 : }
741 :
742 :
743 : // Local Variables:
744 : // compile-command: "gmake OPTLIB=1 TwoSidedShape"
745 : // End:
746 :
747 :
748 : } //# NAMESPACE CASA - END
749 :
|
