Line data Source code
1 : #ifndef CRPROPA_VECTOR3_H
2 : #define CRPROPA_VECTOR3_H
3 :
4 : #include <iostream>
5 : #include <cmath>
6 : #include <vector>
7 : #include <limits>
8 : #include <algorithm>
9 :
10 : namespace crpropa {
11 :
12 : /**
13 : * \addtogroup Core
14 : * @{
15 : */
16 :
17 : /**
18 : @class Vector3
19 : @brief Template class for 3-vectors of type float, double, ...
20 :
21 : Allows accessing and changing the elements x, y, z directly or through the
22 : corresponding get and set methods.
23 :
24 : Angle definitions are
25 : phi [-pi, pi]: azimuthal angle in the x-y plane, 0 pointing in x-direction
26 : theta [0, pi]: zenith angle towards the z axis, 0 pointing in z-direction
27 : */
28 : template<typename T>
29 : class Vector3 {
30 : public:
31 : union {
32 : struct
33 : {
34 : T x;
35 : T y;
36 : T z;
37 : };
38 : T data[3];
39 : };
40 :
41 3187045 : Vector3() : data{0., 0., 0.} {
42 : }
43 :
44 : // avoid creation of default non-conversion constructor
45 20418471 : Vector3(const Vector3 &v) : data{v.data[0], v.data[1], v.data[2]} {
46 0 : }
47 :
48 : // Provides implicit conversion
49 : template<typename U>
50 0 : Vector3(const Vector3<U> &v) {
51 14 : data[0] = v.x;
52 4 : data[1] = v.y;
53 4 : data[2] = v.z;
54 0 : }
55 :
56 : ~Vector3()
57 : {
58 26256343 : }
59 :
60 0 : explicit Vector3(const double *v) {
61 0 : data[0] = v[0];
62 0 : data[1] = v[1];
63 0 : data[2] = v[2];
64 : }
65 :
66 0 : explicit Vector3(const float *v) {
67 0 : data[0] = v[0];
68 0 : data[1] = v[1];
69 0 : data[2] = v[2];
70 : }
71 :
72 12648908 : explicit Vector3(const T &X, const T &Y, const T &Z) : data{X, Y, Z} {
73 : }
74 :
75 18519663 : explicit Vector3(T t) : data{t, t, t} {
76 1 : }
77 :
78 : void setX(const T X) {
79 1 : x = X;
80 : }
81 :
82 : void setY(const T Y) {
83 0 : y = Y;
84 : }
85 :
86 : void setZ(const T Z) {
87 0 : z = Z;
88 : }
89 :
90 : void setXYZ(const T X, const T Y, const T Z) {
91 1351680 : x = X;
92 1351680 : y = Y;
93 1351680 : z = Z;
94 : }
95 :
96 0 : void setR(const T r) {
97 0 : *this *= r / getR();
98 0 : }
99 :
100 3 : void setRThetaPhi(const T r, const T theta, const T phi) {
101 3 : x = r * sin(theta) * cos(phi);
102 3 : y = r * sin(theta) * sin(phi);
103 3 : z = r * cos(theta);
104 3 : }
105 :
106 : T getX() const {
107 370007 : return x;
108 : }
109 :
110 : T getY() const {
111 370006 : return y;
112 : }
113 :
114 : T getZ() const {
115 40006 : return z;
116 : }
117 :
118 : // magnitude (2-norm) of the vector
119 : T getR() const {
120 25117663 : return std::sqrt(x * x + y * y + z * z);
121 : }
122 :
123 : // square of magnitude of the vector
124 : T getR2() const {
125 2883584 : return x * x + y * y + z * z;
126 : }
127 :
128 : // return the azimuth angle
129 19 : T getPhi() const {
130 : T eps = std::numeric_limits < T > ::min();
131 19 : if ((fabs(x) < eps) && (fabs(y) < eps))
132 : return 0.0;
133 : else
134 18 : return std::atan2(y, x);
135 : }
136 :
137 : // return the zenith angle
138 8 : T getTheta() const {
139 : T eps = std::numeric_limits < T > ::min();
140 8 : if ((fabs(x) < eps) && (fabs(y) < eps) && (fabs(z) < eps))
141 : return 0.0;
142 : else
143 8 : return atan2((T) sqrt(x * x + y * y), z);
144 : }
145 :
146 : // return the unit-vector e_r
147 7201361 : Vector3<T> getUnitVector() const {
148 7201361 : return *this / getR();
149 : }
150 :
151 : // return the unit-vector e_theta
152 1 : Vector3<T> getUnitVectorTheta() const {
153 1 : T theta = getTheta();
154 1 : T phi = getPhi();
155 1 : return Vector3<T>(cos(theta) * cos(phi), cos(theta) * sin(phi),
156 1 : -sin(theta));
157 : }
158 :
159 : // return the unit-vector e_phi
160 4 : Vector3<T> getUnitVectorPhi() const {
161 4 : return Vector3<T>(-sin(getPhi()), cos(getPhi()), 0);
162 : }
163 :
164 : // return the angle [0, pi] between the vectors
165 691025 : T getAngleTo(const Vector3<T> &v) const {
166 691025 : T cosdistance = dot(v) / v.getR() / getR();
167 : // In some directions cosdistance is > 1 on some compilers
168 : // This ensures that the correct result is returned
169 691025 : if (cosdistance >= 1.)
170 : return 0;
171 690880 : else if (cosdistance <= -1.)
172 : return M_PI;
173 : else
174 690880 : return acos(cosdistance);
175 : }
176 :
177 : // return true if the angle between the vectors is smaller than a threshold
178 : bool isParallelTo(const Vector3<T> &v, T maxAngle) const {
179 691017 : return getAngleTo(v) < maxAngle;
180 : }
181 :
182 : // linear distance to a given vector
183 110 : T getDistanceTo(const Vector3<T> &point) const {
184 : Vector3<T> d = *this - point;
185 110 : return d.getR();
186 : }
187 :
188 : // return the component parallel to a second vector
189 : // 0 if the second vector has 0 magnitude
190 2 : Vector3<T> getParallelTo(const Vector3<T> &v) const {
191 : T vmag = v.getR();
192 2 : if (vmag == std::numeric_limits < T > ::min())
193 : return Vector3<T>(0.);
194 : return v * dot(v) / vmag;
195 : }
196 :
197 : // return the component perpendicular to a second vector
198 : // 0 if the second vector has 0 magnitude
199 1 : Vector3<T> getPerpendicularTo(const Vector3<T> &v) const {
200 1 : if (v.getR() == std::numeric_limits < T > ::min())
201 : return Vector3<T>(0.);
202 1 : return (*this) - getParallelTo(v);
203 : }
204 :
205 : // rotate the vector around a given axis by a given angle
206 481006 : Vector3<T> getRotated(const Vector3<T> &axis, T angle) const {
207 : Vector3<T> u = axis;
208 481006 : if (u.getR() != 0.)
209 : u = u / u.getR();
210 481006 : T c = cos(angle);
211 481006 : T s = sin(angle);
212 481006 : Vector3<T> Rx(c + u.x * u.x * (1 - c), u.x * u.y * (1 - c) - u.z * s,
213 481006 : u.x * u.z * (1 - c) + u.y * s);
214 481006 : Vector3<T> Ry(u.y * u.x * (1 - c) + u.z * s, c + u.y * u.y * (1 - c),
215 481006 : u.y * u.z * (1 - c) - u.x * s);
216 481006 : Vector3<T> Rz(u.z * u.x * (1 - c) - u.y * s,
217 481006 : u.z * u.y * (1 - c) + u.x * s, c + u.z * u.z * (1 - c));
218 481006 : return Vector3<T>(dot(Rx), dot(Ry), dot(Rz));
219 : }
220 :
221 : // return vector with values limited to the range [lower, upper]
222 0 : Vector3<T> clip(T lower, T upper) const {
223 : Vector3<T> out;
224 0 : out.x = std::max(lower, std::min(x, upper));
225 0 : out.y = std::max(lower, std::min(y, upper));
226 0 : out.z = std::max(lower, std::min(z, upper));
227 0 : return out;
228 : }
229 :
230 : // return vector with absolute values
231 : Vector3<T> abs() const {
232 0 : return Vector3<T>(std::abs(x), std::abs(y), std::abs(z));
233 : }
234 :
235 : // return vector with floored values
236 6 : Vector3<T> floor() const {
237 6 : return Vector3<T>(std::floor(x), std::floor(y), std::floor(z));
238 : }
239 :
240 : // return vector with ceiled values
241 0 : Vector3<T> ceil() const {
242 0 : return Vector3<T>(std::ceil(x), std::ceil(y), std::ceil(z));
243 : }
244 :
245 : // minimum element
246 : T min() const {
247 4 : return std::min(x, std::min(y, z));
248 : }
249 :
250 : // maximum element
251 : T max() const {
252 4 : return std::max(x, std::max(y, z));
253 : }
254 :
255 : // dot product
256 : T dot(const Vector3<T> &v) const {
257 692583 : return x * v.x + y * v.y + z * v.z;
258 : }
259 :
260 : // cross product
261 : Vector3<T> cross(const Vector3<T> &v) const {
262 2132230 : return Vector3<T>(y * v.z - v.y * z, z * v.x - v.z * x,
263 1652230 : x * v.y - v.x * y);
264 : }
265 :
266 : // returns true if all elements of the two vectors are equal
267 : bool operator ==(const Vector3<T> &v) const {
268 33 : if (x != v.x)
269 : return false;
270 33 : if (y != v.y)
271 : return false;
272 33 : if (z != v.z)
273 0 : return false;
274 : return true;
275 : }
276 :
277 : Vector3<T> operator +(const Vector3<T> &v) const {
278 4367246 : return Vector3(x + v.x, y + v.y, z + v.z);
279 : }
280 :
281 : Vector3<T> operator +(const T &f) const {
282 0 : return Vector3(x + f, y + f, z + f);
283 : }
284 :
285 : Vector3<T> operator -(const Vector3<T> &v) const {
286 7540036 : return Vector3(x - v.x, y - v.y, z - v.z);
287 : }
288 :
289 : Vector3<T> operator -(const T &f) const {
290 0 : return Vector3(x - f, y - f, z - f);
291 : }
292 :
293 : // element-wise multiplication
294 : Vector3<T> operator *(const Vector3<T> &v) const {
295 20 : return Vector3(x * v.x, y * v.y, z * v.z);
296 : }
297 :
298 : Vector3<T> operator *(T v) const {
299 1216220 : return Vector3(data[0] * v, data[1] * v, data[2] * v);
300 : }
301 :
302 : // element-wise division
303 : Vector3<T> operator /(const Vector3<T> &v) const {
304 100097 : return Vector3(x / v.x, y / v.y, z / v.z);
305 : }
306 :
307 : Vector3<T> operator /(const T &f) const {
308 6242639 : return Vector3(x / f, y / f, z / f);
309 : }
310 :
311 : // element-wise modulo operation
312 0 : Vector3<T> operator %(const Vector3<T> &v) const {
313 0 : return Vector3(fmod(x, v.x), fmod(y, v.y), fmod(z, v.z));
314 : }
315 :
316 0 : Vector3<T> operator %(const T &f) const {
317 0 : return Vector3(fmod(x, f), fmod(y, f), fmod(z, f));
318 : }
319 :
320 : Vector3<T> &operator -=(const Vector3<T> &v) {
321 0 : data[0] -= v.x;
322 0 : data[1] -= v.y;
323 0 : data[2] -= v.z;
324 : return *this;
325 : }
326 :
327 : Vector3<T> &operator -=(const T &f) {
328 0 : data[0] -= f;
329 0 : data[1] -= f;
330 0 : data[2] -= f;
331 : return *this;
332 : }
333 :
334 : Vector3<T> &operator +=(const Vector3<T> &v) {
335 20664839 : data[0] += v.x;
336 20664839 : data[1] += v.y;
337 20544720 : data[2] += v.z;
338 : return *this;
339 : }
340 :
341 : Vector3<T> &operator +=(const T &f) {
342 0 : data[0] += f;
343 0 : data[1] += f;
344 0 : data[2] += f;
345 : return *this;
346 : }
347 :
348 : // element-wise multiplication
349 : Vector3<T> &operator *=(const Vector3<T> &v) {
350 0 : data[0] *= v.x;
351 0 : data[1] *= v.y;
352 0 : data[2] *= v.z;
353 : return *this;
354 : }
355 :
356 : Vector3<T> &operator *=(const T &f) {
357 3312538 : data[0] *= f;
358 3312538 : data[1] *= f;
359 3312538 : data[2] *= f;
360 : return *this;
361 : }
362 :
363 : // element-wise division
364 : Vector3<T> &operator /=(const Vector3<T> &v) {
365 0 : data[0] /= v.x;
366 0 : data[1] /= v.y;
367 0 : data[2] /= v.z;
368 : return *this;
369 : }
370 :
371 : Vector3<T> &operator /=(const T &f) {
372 691019 : data[0] /= f;
373 691019 : data[1] /= f;
374 691019 : data[2] /= f;
375 : return *this;
376 : }
377 :
378 : // element-wise modulo operation
379 0 : Vector3<T> &operator %=(const Vector3<T> &v) {
380 0 : data[0] = fmod(x, v.x);
381 0 : data[1] = fmod(y, v.y);
382 0 : data[2] = fmod(z, v.z);
383 0 : return *this;
384 : }
385 :
386 1 : Vector3<T> &operator %=(const T &f) {
387 1 : data[0] = fmod(x, f);
388 1 : data[1] = fmod(y, f);
389 1 : data[2] = fmod(z, f);
390 1 : return *this;
391 : }
392 :
393 : Vector3<T> &operator =(const Vector3<T> &v) {
394 64218078 : data[0] = v.x;
395 64218092 : data[1] = v.y;
396 63732471 : data[2] = v.z;
397 0 : return *this;
398 : }
399 :
400 : //Vector3<T> &operator =(Vector3<T> &&v) noexcept {
401 : // data[0] = v.data[0];
402 : // data[1] = v.data[1];
403 : // data[2] = v.data[2];
404 : // return *this;
405 : //}
406 :
407 : Vector3<T> &operator =(const T &f) {
408 : data[0] = f;
409 : data[1] = f;
410 : data[2] = f;
411 : return *this;
412 : }
413 : };
414 :
415 : #ifndef SWIG
416 : template<typename T>
417 38 : inline std::ostream &operator <<(std::ostream &out, const Vector3<T> &v) {
418 114 : out << v.x << " " << v.y << " " << v.z;
419 38 : return out;
420 : }
421 :
422 : template<typename T>
423 : inline std::istream &operator >>(std::istream &in, Vector3<T> &v) {
424 : in >> v.x >> v.y >> v.z;
425 : return in;
426 : }
427 : #endif
428 :
429 : template<typename T>
430 : inline Vector3<T> operator *(T f, const Vector3<T> &v) {
431 3632080 : return Vector3<T>(v.x * f, v.y * f, v.z * f);
432 : }
433 :
434 : typedef Vector3<double> Vector3d;
435 : typedef Vector3<float> Vector3f;
436 :
437 : /** @}*/
438 : } // namespace crpropa
439 :
440 : #endif // CRPROPA_VECTOR3_H
|
