1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
|
//============================================================================
//
// 3-dim real vector class
// $Id$
// Copyright (C) 2004 Georg Drenkhahn
//
// This file is free software; you can redistribute it and/or modify it under
// the terms of the GNU General Public License version 2 as published by the
// Free Software Foundation.
//
//============================================================================
#include <math.h>
#if !defined(NAN)
static inline double nan__()
{
static const unsigned int one = 1;
static const bool BigEndian = (*((unsigned char *) &one) == 0);
static const unsigned char be_nan_bytes[] = { 0x7f, 0xf8, 0, 0, 0, 0, 0, 0 };
static const unsigned char le_nan_bytes[] = { 0, 0, 0, 0, 0, 0, 0xf8, 0x7f };
return *( ( const double * )( BigEndian ? be_nan_bytes : le_nan_bytes ) );
}
# define NAN (::nan__())
#endif
#include <config.h>
#ifdef HAVE_NUMERIC_LIMITS
#include <limits>
#endif
#include <cerrno>
#include "vec3.h"
template<typename T>
vec3<T>& vec3<T>::normalize()
{
T n = norm(*this);
if (n != 0)
{
(*this) /= n;
}
else
{
errno = EDOM; // indicate domain error
// TODO: throw an exception?
}
return *this;
}
template<typename T>
vec3<T>& vec3<T>::rotate(const vec3<T>& r)
{
T phi = norm(r);
if (phi != 0)
{
// part of vector which is parallel to r
vec3<T> par(r*(*this)/(r*r) * r);
// part of vector which is perpendicular to r
vec3<T> perp(*this);
std::valarray<T> &perp_ = perp;
perp -= std::valarray<T>(par);
// rotation direction, size of perp
vec3<T> rotdir(norm(perp) * normalized(crossprod(r,perp)));
*this = par + cos(phi)*perp + sin(phi)*rotdir;
}
return *this;
}
/*--- static member functions ---*/
template<typename T>
T vec3<T>::cos_angle(const vec3<T>& a, const vec3<T>& b)
{
T den = norm(a) * norm(b);
T ret = 0;
// if |a|=0 or |b|=0 then angle is not defined. We return NAN in this case.
if (den != 0.0)
{
ret = a*b/den;
}
else
{
errno = EDOM; // indicate domain error
#ifdef HAVE_NUMERIC_LIMITS
// TODO test
ret = std::numeric_limits<T>::quiet_NaN();
#else
ret = NAN; // return NAN from ISO C99
#endif
}
return ret;
}
template<typename T>
T vec3<T>::angle(const vec3<T>& a, const vec3<T>& b, const vec3<T>& c)
{
// if |a|=0 or |b|=0 then angle is not defined. We return NAN in this case.
T ang = vec3<T>::angle(a,b);
return (crossprod(a,b)*c<0) ?
T(2.*M_PI)-ang : ang;
}
// explicite instantiation
template class vec3<double>;
|
