math.hh 6.32 KB
Newer Older
Dennis Gläser's avatar
Dennis Gläser committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
/*!
 * \file
21
22
 * \ingroup Common
 * \brief Defines some mathematical functions.
Dennis Gläser's avatar
Dennis Gläser committed
23
24
25
26
 */
#ifndef FRACKIT_MATH_HH
#define FRACKIT_MATH_HH

27
#include <cmath>
Dennis Gläser's avatar
Dennis Gläser committed
28
#include <vector>
29

Dennis Gläser's avatar
Dennis Gläser committed
30
#include <frackit/geometry/vector.hh>
Dennis Gläser's avatar
cleanup    
Dennis Gläser committed
31
#include <frackit/precision/precision.hh>
Dennis Gläser's avatar
Dennis Gläser committed
32
33
34
35

namespace Frackit {

/*!
36
 * \ingroup Common
Dennis Gläser's avatar
Dennis Gläser committed
37
38
 * \brief Returns the vector resulting from the cross product
 *        of two 3-dimensional vectors.
Dennis Gläser's avatar
Dennis Gläser committed
39
40
41
42
43
44
45
46
47
48
 */
template<class ctype>
Vector<ctype, 3> crossProduct(const Vector<ctype, 3>& v1,
                              const Vector<ctype, 3>& v2)
{
    return { v1.y()*v2.z() - v1.z()*v2.y(),
             v1.z()*v2.x() - v1.x()*v2.z(),
             v1.x()*v2.y() - v1.y()*v2.x() };
}

Dennis Gläser's avatar
Dennis Gläser committed
49
/*!
50
 * \ingroup Common
Dennis Gläser's avatar
Dennis Gläser committed
51
52
 * \brief Returns the scalar resulting from the cross product
 *        of two 2-dimensional vectors.
Dennis Gläser's avatar
Dennis Gläser committed
53
54
55
56
57
58
 */
template <class ctype>
ctype crossProduct(const Vector<ctype, 2>& v1,
                   const Vector<ctype, 2>& v2)
{ return v1.x()*v2.y() - v1.y()*v2.x(); }

Dennis Gläser's avatar
Dennis Gläser committed
59
/*!
60
 * \ingroup Common
Dennis Gläser's avatar
Dennis Gläser committed
61
62
 * \brief Returns the scalar resulting from the box product
 *        of three 3-dimensional vectors.
Dennis Gläser's avatar
Dennis Gläser committed
63
64
65
66
67
68
69
70
 */
template<class ctype>
ctype boxProduct(const Vector<ctype, 3>& v1,
                 const Vector<ctype, 3>& v2,
                 const Vector<ctype, 3>& v3)
{ return v1*crossProduct(v2, v3); }

/*!
71
 * \ingroup Common
Dennis Gläser's avatar
Dennis Gläser committed
72
73
 * \brief Returns true if the given 3-dimensional
 *        vectors form a right-hand system.
Dennis Gläser's avatar
Dennis Gläser committed
74
75
76
77
78
79
 */
template<class ctype>
bool isRightHandSystem(const Vector<ctype, 3>& v1,
                       const Vector<ctype, 3>& v2,
                       const Vector<ctype, 3>& v3)
{
Dennis Gläser's avatar
cleanup    
Dennis Gläser committed
80
81
82
    assert(std::abs(boxProduct(v1, v2, v3)) > v1.length()*Precision<ctype>::confusion()
                                              *v2.length()*Precision<ctype>::confusion()
                                              *v3.length()*Precision<ctype>::confusion());
Dennis Gläser's avatar
Dennis Gläser committed
83
84
85
    return boxProduct(v1, v2, v3) > 0.0;
}

86
/*!
87
 * \ingroup Common
Dennis Gläser's avatar
Dennis Gläser committed
88
 * \brief Converts a given angle in radians to degrees.
89
90
91
92
93
94
 */
template<class ctype>
ctype toDegrees(const ctype radians)
{ return radians*180.0/M_PI; }

/*!
95
 * \ingroup Common
Dennis Gläser's avatar
Dennis Gläser committed
96
 * \brief Converts a given angle in degrees to radians.
97
98
99
100
101
 */
template<class ctype>
ctype toRadians(const ctype degrees)
{ return degrees*M_PI/180.0; }

102
/*!
103
 * \ingroup Common
Dennis Gläser's avatar
Dennis Gläser committed
104
105
106
107
108
109
 * \brief Rotate a vector around an axis by a given angle.
 * \param v The vector to be rotated
 * \param axis The rotation axis
 * \param angle The rotation angle
 * \note The rotation matrix is taken from:
 *       https://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
 */
template<class ctype>
void rotate(Vector<ctype, 3>& v,
            const Direction<ctype, 3>& axis,
            const ctype angle)
{
    using std::sin;
    using std::cos;

    const auto cosPhi = cos(angle);
    const auto sinPhi = sin(angle);

    Vector<ctype, 3> Rx(cosPhi + axis.x()*axis.x()*(1.0 - cosPhi),
                        axis.x()*axis.y()*(1.0 - cosPhi) - axis.z()*sinPhi,
                        axis.x()*axis.z()*(1.0 - cosPhi) - axis.y()*sinPhi);
    Vector<ctype, 3> Ry(axis.y()*axis.x()*(1.0 - cosPhi) + axis.z()*sinPhi,
                        cosPhi + axis.y()*axis.y()*(1.0 - cosPhi),
                        axis.y()*axis.z()*(1.0 - cosPhi) - axis.x()*sinPhi);
    Vector<ctype, 3> Rz(axis.z()*axis.x()*(1.0 - cosPhi) - axis.y()*sinPhi,
                        axis.z()*axis.y()*(1.0 - cosPhi) + axis.x()*sinPhi,
                        cosPhi + axis.z()*axis.z()*(1.0 - cosPhi));
    v = Vector<ctype, 3>(Rx*v, Ry*v, Rz*v);
}

/*!
135
 * \ingroup Common
Dennis Gläser's avatar
Dennis Gläser committed
136
137
138
139
140
141
142
143
144
 * \brief Rotate several vectors around an axis by a given angle.
 * \param vectors The vectors to be rotated
 * \param axis The rotation axis
 * \param angle The rotation angle
 * \note This overload is more efficient than calling the overload
 *       for a single vector several times, as the rotation matrix
 *       is only constructed once!
 * \note The rotation matrix is taken from:
 *       https://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
 */
template<class ctype>
void rotate(std::vector<Vector<ctype, 3>>& vectors,
            const Direction<ctype, 3>& axis,
            const ctype angle)
{
    using std::sin;
    using std::cos;

    const auto cosPhi = cos(angle);
    const auto sinPhi = sin(angle);

    Vector<ctype, 3> Rx(cosPhi + axis.x()*axis.x()*(1.0 - cosPhi),
                        axis.x()*axis.y()*(1.0 - cosPhi) - axis.z()*sinPhi,
                        axis.x()*axis.z()*(1.0 - cosPhi) - axis.y()*sinPhi);
    Vector<ctype, 3> Ry(axis.y()*axis.x()*(1.0 - cosPhi) + axis.z()*sinPhi,
                        cosPhi + axis.y()*axis.y()*(1.0 - cosPhi),
                        axis.y()*axis.z()*(1.0 - cosPhi) - axis.x()*sinPhi);
    Vector<ctype, 3> Rz(axis.z()*axis.x()*(1.0 - cosPhi) - axis.y()*sinPhi,
                        axis.z()*axis.y()*(1.0 - cosPhi) + axis.x()*sinPhi,
                        cosPhi + axis.z()*axis.z()*(1.0 - cosPhi));

    for (auto& v : vectors)
        v = Vector<ctype, 3>(Rx*v, Ry*v, Rz*v);
}

Dennis Gläser's avatar
Dennis Gläser committed
171
172
173
} // end namespace Frackit

#endif // FRACKIT_MATH_HH