LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
Angle.h
Go to the documentation of this file.
1 #if !defined(LSST_AFW_GEOM_ANGLE_H)
2 #define LSST_AFW_GEOM_ANGLE_H
3 
4 #include <limits>
5 #include <iostream>
6 #include <boost/math/constants/constants.hpp>
7 #include <cmath>
8 #include <boost/static_assert.hpp>
9 
10 namespace lsst { namespace afw { namespace geom {
11 
12 /************************************************************************************************************/
13 /*
14  * None of C99, C++98, and C++0x define M_PI, so we'll do it ourselves
15  */
16 #pragma clang diagnostic push
17 #pragma clang diagnostic ignored "-Wunused-variable"
18 double const PI = boost::math::constants::pi<double>();
19 double const TWOPI = boost::math::constants::pi<double>() * 2.0;
20 double const HALFPI = boost::math::constants::pi<double>() * 0.5;
21 double const ONE_OVER_PI = 1.0 / boost::math::constants::pi<double>();
22 double const SQRTPI = sqrt(boost::math::constants::pi<double>());
23 double const INVSQRTPI = 1.0/sqrt(boost::math::constants::pi<double>());
24 double const ROOT2 = boost::math::constants::root_two<double>(); // sqrt(2)
25 #pragma clang diagnostic pop
26 
27 // These shouldn't be necessary if the Angle class is used, but sometimes you just need
28 // them. Better to define them once here than have *180/PI throughout the code...
29 inline double degToRad(double x) {
30  return x * PI / 180.;
31 }
32 inline double radToDeg(double x) {
33  return x * 180. / PI;
34 }
35 inline double radToArcsec(double x) {
36  return x * 3600. * 180. / PI;
37 }
38 inline double radToMas(double x) {
39  return x * 1000. * 3600. * 180. / PI;
40 }
41 inline double arcsecToRad(double x) {
42  return (x / 3600.) * PI / 180.;
43 }
44 inline double masToRad(double x) {
45  return (x / (1000. * 3600.)) * PI / 180.;
46 }
47 
48 // NOTE, if you add things here, you must also add them to
49 // python/lsst/afw/geom/__init__.py
50 // (if you want them to accessible from python)
51 
52 #if 0 && !defined(M_PI) // a good idea, but with ramifications
53 # define M_PI ::lsst::afw::geom::PI
54 #endif
55 
56 /************************************************************************************************************/
57 
58 class Angle;
71 class AngleUnit {
72  friend class Angle;
73  template<typename T> friend const Angle operator *(T lhs, AngleUnit const rhs);
74 public:
75  explicit AngleUnit(double val) : _val(val) {}
76 
77  bool operator==(AngleUnit const &rhs) const;
78 private:
79  double _val;
80 };
81 
83  return (_val == rhs._val);
84 }
85 
86 // NOTE, if you add things here, remember to also add them to
87 // python/lsst/afw/geom/__init__.py
88 
89 // swig likes this way of initialising the constant, so don't mess with it;
90 // N.b. swig 1.3 doesn't like PI/(60*180)
91 AngleUnit const radians = AngleUnit(1.0);
92 AngleUnit const degrees = AngleUnit(PI/180.0); // constant with units of degrees
93 AngleUnit const hours = AngleUnit(PI*15.0/180.0); // constant with units of hours
94 AngleUnit const arcminutes = AngleUnit(PI/60/180.0); // constant with units of arcminutes
95 AngleUnit const arcseconds = AngleUnit(PI/180.0/3600.0); // constant with units of arcseconds
96 
97 /************************************************************************************************************/
104 class Angle {
105  friend class AngleUnit;
106 public:
108  explicit Angle(double val, AngleUnit units=radians) : _val(val*units._val) {}
109  Angle() : _val(0) {}
111  Angle(Angle const& other) : _val(other._val) {}
113  operator double() const { return _val; }
115  //operator float() const { return _val; }
116 
118  double asAngularUnits(AngleUnit const& units) const {
119  return _val/units._val;
120  }
122  double asRadians() const { return asAngularUnits(radians); }
124  double asDegrees() const { return asAngularUnits(degrees); }
126  double asHours() const { return asAngularUnits(hours); }
128  double asArcminutes() const { return asAngularUnits(arcminutes); }
130  double asArcseconds() const { return asAngularUnits(arcseconds); }
131 
132  double toUnitSphereDistanceSquared() const { return 2. * (1. - std::cos(asRadians())); }
133  // == 4.0 * pow(std::sin(0.5 * asRadians()), 2.0)
135  return (std::acos(1. - d2/2.)) * radians;
136  // == 2.0 * asin(0.5 * sqrt(d2))
137  }
138 
146  void wrap() {
147  _val = std::fmod(_val, TWOPI);
148  // _val is now in the range (-TWOPI, TWOPI)
149  if (_val < 0.0)
150  _val += TWOPI;
151  // if _val is small enough, adding 2 pi gives 2 pi
152  if (_val >= TWOPI)
153  _val = 0.0;
154  }
155 
162  void wrapCtr() {
163  _val = std::fmod(_val, TWOPI);
164  // _val is now in the range [-TWOPI, TWOPI]
165  if (_val < -PI) {
166  _val += TWOPI;
167  if (_val >= PI) {
168  // handle roundoff error, however unlikely
169  _val = -PI;
170  }
171  } else if (_val >= PI) {
172  _val -= TWOPI;
173  if (_val < -PI) {
174  // handle roundoff error, however unlikely
175  _val = -PI;
176  }
177  }
178  }
179 
186  void wrapNear(
187  Angle const & refAng
188  ) {
189  // compute this = (this - refAng).wrapCtr() + refAng
190  // which is correct except for roundoff error at the edges
191  double refAngRad = refAng.asRadians();
192  *this -= refAng;
193  wrapCtr();
194  _val += refAngRad;
195 
196  // roundoff can cause slightly out-of-range values; fix those
197  if (_val - refAngRad >= PI) {
198  _val -= TWOPI;
199  }
200  // maximum relative roundoff error for subtraction is 2 epsilon
201  if (_val - refAngRad < -PI) {
202  _val -= _val * 2.0 * std::numeric_limits<double>::epsilon();
203  }
204  }
205 
206 
207 #define ANGLE_OPUP_TYPE(OP, TYPE) \
208  Angle& operator OP(TYPE const& d) { \
209  _val OP d; \
210  return *this; \
211  }
212 
213 ANGLE_OPUP_TYPE(*=, double)
215 ANGLE_OPUP_TYPE(+=, double)
217 ANGLE_OPUP_TYPE(-=, double)
219 
220 #undef ANGLE_OPUP_TYPE
221 
222 #define ANGLE_COMP(OP) \
223  bool operator OP ( const Angle& rhs ) { \
224  return _val OP rhs._val; \
225  }
226 
233 
234 #undef ANGLE_COMP
235 
236 private:
237  double _val;
238 };
239 
240 Angle const NullAngle = Angle(-1000000., degrees);
241 
242 
243 /************************************************************************************************************/
244 /*
245  * Operators for Angles.
246  *
247  * N.b. We need both int and double versions to avoid ambiguous overloading due to implicit conversion of
248  * Angle to double
249  */
250 #define ANGLE_OP(OP) \
251  inline const Angle operator OP(Angle const a, Angle const d) { \
252  return Angle(static_cast<double>(a) OP static_cast<double>(d)); \
253  }
254 
255 #define ANGLE_OP_TYPE(OP, TYPE) \
256  inline const Angle operator OP(Angle const a, TYPE d) { \
257  return Angle(static_cast<double>(a) OP d); \
258  } \
259  \
260  inline const Angle operator OP(TYPE d, Angle const a) { \
261  return Angle(d OP static_cast<double>(a)); \
262  }
263 
267 ANGLE_OP_TYPE(*, double)
269 
270 #undef ANGLE_OP
271 #undef ANGLE_OP_TYPE
272 
273 // Division is different. Don't allow division by an Angle
274 inline const Angle operator /(Angle const a, int d) {
275  return Angle(static_cast<double>(a)/d);
276 }
277 
278 inline const Angle operator /(Angle const a, double d) {
279  return Angle(static_cast<double>(a)/d);
280 }
281 
282 template<typename T>
283  double operator /(T const lhs, Angle const rhs);
284 
285 /************************************************************************************************************/
289 template<typename T>
290 inline bool isAngle(T) {
291  return false;
292 };
293 
294 inline bool isAngle(Angle const&) {
295  return true;
296 };
297 
298 /************************************************************************************************************/
302 template<typename T>
303 inline
304 const Angle operator *(T lhs,
305  AngleUnit const rhs
306  ) {
307  BOOST_STATIC_ASSERT_MSG(std::numeric_limits<T>::is_specialized,
308  "Only numeric types may be converted to Angles using degrees/radians!");
309  return Angle(lhs*rhs._val);
310 }
314 std::ostream& operator<<(std::ostream &s,
315  Angle const a
316  );
317 
318 }}}
319 #endif // if !defined(LSST_AFW_GEOM_ANGLE_H)
double toUnitSphereDistanceSquared() const
Definition: Angle.h:132
const Angle operator/(Angle const a, int d)
Definition: Angle.h:274
double degToRad(double x)
Definition: Angle.h:29
double asArcminutes() const
Definition: Angle.h:128
double asAngularUnits(AngleUnit const &units) const
Definition: Angle.h:118
AngleUnit const radians
constant with units of radians
Definition: Angle.h:91
std::ostream & operator<<(std::ostream &os, lsst::afw::geom::AffineTransform const &transform)
double asArcseconds() const
Definition: Angle.h:130
AngleUnit const hours
Definition: Angle.h:93
double const TWOPI
Definition: Angle.h:19
double const SQRTPI
Definition: Angle.h:22
A class used to convert scalar POD types such as double to Angle.
Definition: Angle.h:71
Angle(Angle const &other)
Definition: Angle.h:111
double asDegrees() const
Definition: Angle.h:124
double radToDeg(double x)
Definition: Angle.h:32
double asRadians() const
Definition: Angle.h:122
#define ANGLE_COMP(OP)
Definition: Angle.h:222
AngleUnit const degrees
Definition: Angle.h:92
void wrapNear(Angle const &refAng)
Definition: Angle.h:186
const Angle operator*(Angle const a, Angle const d)
Definition: Angle.h:266
double radToArcsec(double x)
Definition: Angle.h:35
static Angle fromUnitSphereDistanceSquared(double d2)
Definition: Angle.h:134
bool operator==(AngleUnit const &rhs) const
Definition: Angle.h:82
bool isAngle(T)
Allow a user to check if they have an angle (yes; they could do this themselves via trivial TMP) ...
Definition: Angle.h:290
double const ONE_OVER_PI
Definition: Angle.h:21
AngleUnit(double val)
Definition: Angle.h:75
Angle(double val, AngleUnit units=radians)
Definition: Angle.h:108
int x
#define ANGLE_OP(OP)
Definition: Angle.h:250
lsst::afw::geom::Angle Angle
Definition: misc.h:39
AngleUnit const arcminutes
Definition: Angle.h:94
double const HALFPI
Definition: Angle.h:20
double asHours() const
Definition: Angle.h:126
AngleUnit const arcseconds
Definition: Angle.h:95
double const ROOT2
Definition: Angle.h:24
double arcsecToRad(double x)
Definition: Angle.h:41
double radToMas(double x)
Definition: Angle.h:38
double masToRad(double x)
Definition: Angle.h:44
double const PI
The ratio of a circle&#39;s circumference to diameter.
Definition: Angle.h:18
friend const Angle operator*(T lhs, AngleUnit const rhs)
Use AngleUnit to convert a POD (e.g. int, double) to an Angle; e.g. 180*afwGeomdegrees.
Definition: Angle.h:304
#define ANGLE_OPUP_TYPE(OP, TYPE)
Definition: Angle.h:207
ImageT val
Definition: CR.cc:154
double const INVSQRTPI
Definition: Angle.h:23
#define ANGLE_OP_TYPE(OP, TYPE)
Definition: Angle.h:255
Angle const NullAngle
Definition: Angle.h:240