LOOS  v2.3.2
Coord.hpp
1 /*
2  This file is part of LOOS.
3 
4  LOOS (Lightweight Object-Oriented Structure library)
5  Copyright (c) 2008, Tod D. Romo, Alan Grossfield
6  Department of Biochemistry and Biophysics
7  School of Medicine & Dentistry, University of Rochester
8 
9  This package (LOOS) is free software: you can redistribute it and/or modify
10  it under the terms of the GNU General Public License as published by
11  the Free Software Foundation under version 3 of the License.
12 
13  This package is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with this program. If not, see <http://www.gnu.org/licenses/>.
20 */
21 
22 
23 
24 #if !defined(LOOS_COORDS_HPP)
25 #define LOOS_COORDS_HPP
26 
27 #include <iostream>
28 #include <string>
29 #include <stdexcept>
30 #include <utils_random.hpp>
31 
32 namespace loos {
33 
34  // Forward declarations for matrix-vector multiply...
35 
36  template<class T> class Coord;
37  template<class T> class Matrix44;
38  template<class T> Coord<T> operator*(const Matrix44<T>&, const Coord<T>&);
39 
40 
42 
81  template <class T>
82  class Coord {
83  enum { X=0, Y=1, Z=2, MAXCOORD } CoordIndex;
84 
86  static const double epsilon;
87 
88  public:
89 
90  typedef T element_type;
91 
92  // Constructors
93 
94  Coord() { zero(); }
95 
96 
97  Coord(const T ax, const T ay, const T az) { set(ax, ay, az); }
98 
99  Coord(const Coord<T>& o) { copy(o); }
100 
101  Coord(const T x) {
102  uint i;
103  for (i=0; i<MAXCOORD; ++i)
104  v[i] = x;
105  v[i] = 1.0;
106  }
107 
108 
109  // ---------------------------------------
110  // Accessors
111 
112  void x(const T ax) { v[X] = ax; }
113  void y(const T ay) { v[Y] = ay; }
114  void z(const T az) { v[Z] = az; }
115 
116 
117  const T& x(void) const { return(v[X]); }
118  const T& y(void) const { return(v[Y]); }
119  const T& z(void) const { return(v[Z]); }
120 
121 
122 #if !defined(SWIG)
123  T& x(void) { return(v[X]); }
124  T& y(void) { return(v[Y]); }
125  T& z(void) { return(v[Z]); }
126 
127 
129  T& operator[](const unsigned int i) throw(std::out_of_range) {
130  if (i>=MAXCOORD)
131  throw std::out_of_range("Index into Coord<T> is out of range.");
132  return(v[i]);
133  }
134 
136  const T& operator[](const unsigned int i) const throw(std::out_of_range) {
137  if (i>=MAXCOORD)
138  throw std::out_of_range("Index into Coord<T> is out of range.");
139  return(v[i]);
140  }
141 #endif // !defined(SWIG)
142 
144  void set(const T x, const T y, const T z) {
145  v[X] = x;
146  v[Y] = y;
147  v[Z] = z;
148  v[MAXCOORD] = 1;
149  }
150 
151 
152  // ---------------------------------------
153  // I/O
154 
155 #if !defined(SWIG)
156  friend std::ostream& operator<<(std::ostream& os, const Coord<T>&o) {
158  os << "(";
159  for (uint i=0; i<MAXCOORD; ++i)
160  os << o.v[i] << (i < MAXCOORD-1 ? "," : "");
161  os << ")";
162  return(os);
163  }
164 
165  friend std::istream& operator>>(std::istream& is, Coord<T>& i) throw(std::runtime_error){
166  char c;
167  is >> c;
168  if (c != '(')
169  throw(std::runtime_error("Invalid Coord conversion"));
170 
171  is >> i.x();
172  is >> c;
173  if (c != ',')
174  throw(std::runtime_error("Invalid Coord conversion"));
175 
176 
177  is >> i.y();
178  is >> c;
179  if (c != ',')
180  throw(std::runtime_error("Invalid Coord conversion"));
181 
182 
183  is >> i.z();
184  is >> c;
185  if (c != ')')
186  throw(std::runtime_error("Invalid Coord conversion"));
187 
188  return(is);
189  }
190 
191  // ---------------------------------------
192  // Operators
193 
194 
195 
196  const Coord<T>& operator=(const Coord<T>& c) { copy(c); return(*this); }
197 
198 
200  friend Coord<T> operator+(const T lhs, const Coord<T>& rhs) {
201  Coord<T> res(rhs);
202  res += lhs;
203  return(res);
204  }
205 
206 
207 
209  friend Coord<T> operator-(const T lhs, const Coord<T>& rhs) {
210  Coord<T> res;
211 
212  for (uint i=0; i<MAXCOORD; ++i)
213  res.v[i] = lhs - rhs.v[i];
214 
215  return(res);
216  }
217 
218 
220  friend Coord<T> operator*<>(const Matrix44<T>&, const Coord<T>&);
221 
222 
224  friend Coord<T> operator*(const T lhs, const Coord<T>& rhs) {
225  Coord<T> res(rhs);
226  res *= lhs;
227  return(res);
228  }
229 
230 
231 
233  Coord<T>& operator/=(const T rhs) {
234  for (uint i=0; i<MAXCOORD; ++i)
235  v[i] /= rhs;
236 
237  return(*this);
238  }
239 
240  Coord<T> operator/(const T rhs) const {
241  Coord<T> res(*this);
242  res /= rhs;
243  return(res);
244  }
245 
247  friend Coord<T> operator/(const T lhs, const Coord<T>& rhs) {
248  Coord<T> res;
249 
250  for (uint i=0; i<MAXCOORD; ++i)
251  res.v[i] = lhs / rhs.v[i];
252 
253  return(res);
254  }
255 
256 
257 
258 #endif // !defined(SWIG)
259 
260 
261  Coord<T>& operator+=(const T rhs) {
262  for (uint i=0; i<MAXCOORD; ++i)
263  v[i] += rhs;
264 
265  return(*this);
266  }
267 
268  Coord<T> operator+(const T rhs) const {
269  Coord<T> res(*this);
270  res += rhs;
271  return(res);
272  }
273 
275  Coord<T>& operator+=(const Coord<T>& rhs) {
276  for (uint i=0; i<MAXCOORD; i++)
277  v[i] += rhs.v[i];
278 
279  return(*this);
280  }
281 
282  Coord<T> operator+(const Coord<T>& rhs) const {
283  Coord<T> res(*this);
284  res += rhs;
285  return(res);
286  }
287 
289  Coord<T>& operator-=(const Coord<T>& rhs) {
290  for (uint i=0; i<MAXCOORD; ++i)
291  v[i] -= rhs.v[i];
292 
293  return(*this);
294  }
295 
296  Coord<T> operator-(const Coord<T>& rhs) const {
297  Coord<T> res(*this);
298  res -= rhs;
299  return(res);
300  }
301 
302  Coord<T>& operator-=(const T rhs) {
303  for (uint i=0; i<MAXCOORD; ++i)
304  v[i] -= rhs;
305 
306  return(*this);
307  }
308 
309  Coord<T> operator-(const T rhs) const {
310  Coord<T> res(*this);
311 
312  res -= rhs;
313  return(res);
314  }
315 
318  Coord<T> res(*this);
319 
320  for (uint i=0; i<MAXCOORD; ++i)
321  res.v[i] = -res.v[i];
322  return(res);
323  }
324 
325 
326 
328  Coord<T>& operator*=(const T rhs) {
329 
330  for (uint i=0; i<MAXCOORD; ++i)
331  v[i] *= rhs;
332 
333  return(*this);
334  }
335 
336  Coord<T> operator*(const T rhs) const {
337  Coord<T> res(*this);
338  res *= rhs;
339  return(res);
340  }
341 
342 
344  T dot(const Coord<T>& rhs) const {
345  T s = v[0] * rhs.v[0];
346  for (uint i=1; i<MAXCOORD; ++i)
347  s += v[i] * rhs.v[i];
348 
349  return(s);
350  }
351 
352 
353 
354  T operator*(const Coord<T>rhs) const {
355  return(dot(rhs));
356  }
357 
358 
360  Coord<T> cross(const Coord<T>& rhs) const {
361  Coord<T> res;
362  res.v[X] = v[Y] * rhs.v[Z] - v[Z] * rhs.v[Y];
363  res.v[Y] = v[Z] * rhs.v[X] - v[X] * rhs.v[Z];
364  res.v[Z] = v[X] * rhs.v[Y] - v[Y] * rhs.v[X];
365 
366  return(res);
367  }
368 
370  Coord<T>& operator^=(const Coord<T>& rhs) {
371  Coord<T> tmp = cross(rhs);
372  *this = tmp;
373  return(*this);
374  }
375 
377  Coord<T> operator^(const Coord<T>& rhs) const {
378  return(cross(rhs));
379  }
380 
381 
383  Coord<T>& operator%=(const Coord<T>& rhs) {
384 
385  for (uint i=0; i<MAXCOORD; ++i)
386  v[i] = fmod(v[i], rhs.v[i]);
387 
388  return(*this);
389  }
390 
391  Coord<T> operator%(const Coord<T>& rhs) const {
392  Coord<T> res(*this);
393 
394  res %= rhs;
395  return(res);
396  }
397 
398 
399 
400  //-----------------------------
401  // Misc
402 
404  void reimage(const Coord<T>& box) {
405  for (uint i=0; i<MAXCOORD; ++i) {
406  int n = (int)(fabs(v[i]) / box.v[i] + 0.5);
407  v[i] = (v[i] >= 0) ? v[i] - n*(box.v[i]) : v[i] + n*(box.v[i]);
408  }
409  }
410 
411 
413  double length2(void) const {
414  double d = 0.0;
415 
416  for (uint i=0; i<MAXCOORD; ++i)
417  d += v[i]*v[i];
418 
419  return(d);
420  }
421 
423  double length(void) const {
424  return(sqrt(length2()));
425  }
426 
427 
429  double distance2(const Coord<T>& o) const {
430  Coord<T> d = o - *this;
431  return(d.length2());
432  }
433 
436  double distance2(const Coord<T>& o, const Coord<T>& box) const {
437  Coord<T> d = o - *this;
438  d.reimage(box);
439  return(d.length2());
440  }
441 
443  double distance(const Coord<T>& o) const {
444  return(sqrt(distance2(o)));
445  }
446 
449  double distance(const Coord<T>& o, const Coord<T>& box) const {
450  return(sqrt(distance2(o, box)));
451  }
452 
454 
461  void random(void) {
462  base_generator_type& rng = rng_singleton();
463  boost::uniform_on_sphere<> uni(3);
464  boost::variate_generator<base_generator_type&,
465  boost::uniform_on_sphere<> > func(rng, uni);
466  std::vector<double> vec = func();
467  for (uint i=0; i<MAXCOORD; i++) {
468  v[i] = static_cast<T>(vec[i]);
469  }
470 
471  }
472 
473 
474 
476  void zero(void) {
477  uint i;
478  for (i=0; i<MAXCOORD; ++i)
479  v[i] = 0;
480  v[i] = 1;
481  }
482 
484  bool operator==(const Coord<T>& rhs) const {
485  double d = distance2(rhs);
486  if (d < epsilon * epsilon)
487  return(true);
488  return(false);
489  }
490 
492  bool operator!=(const Coord<T>& rhs) const {
493  return(!(operator==(rhs)));
494  }
495 
496 
497 
498 
499  private:
500  void copy(const Coord<T>& c) {
501  if (&c == this)
502  return;
503 
504  uint i;
505  for (i=0; i<MAXCOORD; ++i)
506  v[i] = c.v[i];
507 
508  v[i] = 1;
509  };
510 
511 
512  T v[MAXCOORD+1];
513  };
514 
515  template<typename T>
516  const double Coord<T>::epsilon = 1e-15;
517 
518 }
519 
520 #endif
T & operator[](const unsigned int i)
Retrieve an element from the Coord with range-checking.
Definition: Coord.hpp:129
friend Coord< T > operator*(const T lhs, const Coord< T > &rhs)
Handle T * Coord
Definition: Coord.hpp:224
void set(const T x, const T y, const T z)
Short-cut to set the cartesian coordinates...
Definition: Coord.hpp:144
const T & operator[](const unsigned int i) const
Retrieve an element from a const Coord with range-checking.
Definition: Coord.hpp:136
Coord< T > operator^(const Coord< T > &rhs) const
Cross-product (note precedence issues)
Definition: Coord.hpp:377
Basic 3-D coordinates class.
Definition: Coord.hpp:36
void random(void)
Generate a random vector on a unit sphere.
Definition: Coord.hpp:461
T dot(const Coord< T > &rhs) const
Dot product.
Definition: Coord.hpp:344
double length(void) const
Length of the coordinate (as a vector)
Definition: Coord.hpp:423
Coord< T > & operator^=(const Coord< T > &rhs)
Mutating cross-product (note precedence issues)
Definition: Coord.hpp:370
double distance2(const Coord< T > &o) const
Distance squared between two coordinates.
Definition: Coord.hpp:429
Specialized 4x4 Matrix class for handling coordinate transforms.
Definition: Coord.hpp:37
bool operator!=(const Coord< T > &rhs) const
Compute inequality based on ! ==.
Definition: Coord.hpp:492
Coord< T > & operator*=(const T rhs)
Multiplication by a constant.
Definition: Coord.hpp:328
Coord< T > & operator/=(const T rhs)
Division by a constant.
Definition: Coord.hpp:233
double distance(const Coord< T > &o) const
Distance between two coordinates.
Definition: Coord.hpp:443
Coord< T > operator*(const Matrix44< T > &, const Coord< T > &)
Definition: Matrix44.hpp:244
void zero(void)
Zero out the coordinates (while keeping it homogenous)
Definition: Coord.hpp:476
Coord< T > & operator-=(const Coord< T > &rhs)
Subtraction.
Definition: Coord.hpp:289
base_generator_type & rng_singleton(void)
Suite-wide random number generator singleton.
Coord< T > operator-()
Unary negation.
Definition: Coord.hpp:317
Coord< T > & operator+=(const Coord< T > &rhs)
Handle addition.
Definition: Coord.hpp:275
bool operator==(const Coord< T > &rhs) const
Compute equality based on norm(u-v) < epsilon.
Definition: Coord.hpp:484
Coord< T > & operator%=(const Coord< T > &rhs)
Modulo of each component of the Coord with a constant.
Definition: Coord.hpp:383
double distance2(const Coord< T > &o, const Coord< T > &box) const
Definition: Coord.hpp:436
Namespace for most things not already encapsulated within a class.
double length2(void) const
Length of the Coord (as a vector) squared.
Definition: Coord.hpp:413
Coord< T > cross(const Coord< T > &rhs) const
Cross-product. Returns a new Coord
Definition: Coord.hpp:360
friend Coord< T > operator+(const T lhs, const Coord< T > &rhs)
Handle the case of T + Coord
Definition: Coord.hpp:200
friend Coord< T > operator/(const T lhs, const Coord< T > &rhs)
T / Coord case... This may not actually be a good idea?
Definition: Coord.hpp:247
friend Coord< T > operator-(const T lhs, const Coord< T > &rhs)
Handle the case of T - Coord
Definition: Coord.hpp:209
double distance(const Coord< T > &o, const Coord< T > &box) const
Definition: Coord.hpp:449
void reimage(const Coord< T > &box)
Handle coordinates with periodic boundary conditions.
Definition: Coord.hpp:404
friend Coord< T > operator*(const Matrix44< T > &, const Coord< T > &)
For matrix-vector multiply.
Definition: Matrix44.hpp:244