Line data Source code
1 : #ifndef CRPROPA_VECTOR3_H
2 : #define CRPROPA_VECTOR3_H
3 : #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
4 :
5 : #include <iostream>
6 : #include <cmath>
7 : #include <vector>
8 : #include <limits>
9 : #include <algorithm>
10 : #ifdef CRPROPA_HAVE_PYTHON
11 : #include <Python.h>
12 : #include <numpy/arrayobject.h>
13 : #endif // CRPROPA_HAVE_PYTHON
14 : #include <unistd.h>
15 :
16 :
17 : namespace crpropa {
18 :
19 : /**
20 : * \addtogroup Core
21 : * @{
22 : */
23 :
24 : /**
25 : @class Vector3
26 : @brief Template class for 3-vectors of type float, double, ...
27 :
28 : Allows accessing and changing the elements x, y, z directly or through the
29 : corresponding get and set methods.
30 :
31 : Angle definitions are
32 : phi [-pi, pi]: azimuthal angle in the x-y plane, 0 pointing in x-direction
33 : theta [0, pi]: zenith angle towards the z axis, 0 pointing in z-direction
34 : */
35 : template<typename T>
36 : class Vector3 {
37 : public:
38 : union {
39 : struct
40 : {
41 : T x;
42 : T y;
43 : T z;
44 : };
45 : T data[3];
46 : };
47 :
48 3187095 : Vector3() : data{0., 0., 0.} {
49 : }
50 :
51 : // avoid creation of default non-conversion constructor
52 20388214 : Vector3(const Vector3 &v) : data{v.data[0], v.data[1], v.data[2]} {
53 0 : }
54 :
55 : // Provides implicit conversion
56 : template<typename U>
57 0 : Vector3(const Vector3<U> &v) {
58 14 : data[0] = v.x;
59 4 : data[1] = v.y;
60 4 : data[2] = v.z;
61 0 : }
62 :
63 : ~Vector3()
64 : {
65 26231545 : }
66 :
67 0 : explicit Vector3(const double *v) {
68 0 : data[0] = v[0];
69 0 : data[1] = v[1];
70 0 : data[2] = v[2];
71 : }
72 :
73 0 : explicit Vector3(const float *v) {
74 0 : data[0] = v[0];
75 0 : data[1] = v[1];
76 0 : data[2] = v[2];
77 : }
78 :
79 12618760 : explicit Vector3(const T &X, const T &Y, const T &Z) : data{X, Y, Z} {
80 : }
81 :
82 18524158 : explicit Vector3(T t) : data{t, t, t} {
83 1 : }
84 :
85 : void setX(const T X) {
86 1 : x = X;
87 : }
88 :
89 : void setY(const T Y) {
90 0 : y = Y;
91 : }
92 :
93 : void setZ(const T Z) {
94 0 : z = Z;
95 : }
96 :
97 : void setXYZ(const T X, const T Y, const T Z) {
98 1351680 : x = X;
99 1351680 : y = Y;
100 1351680 : z = Z;
101 : }
102 :
103 0 : void setR(const T r) {
104 0 : *this *= r / getR();
105 0 : }
106 :
107 3 : void setRThetaPhi(const T r, const T theta, const T phi) {
108 3 : x = r * sin(theta) * cos(phi);
109 3 : y = r * sin(theta) * sin(phi);
110 3 : z = r * cos(theta);
111 3 : }
112 :
113 : T getX() const {
114 370011 : return x;
115 : }
116 :
117 : T getY() const {
118 370010 : return y;
119 : }
120 :
121 : T getZ() const {
122 40010 : return z;
123 : }
124 :
125 : // magnitude (2-norm) of the vector
126 : T getR() const {
127 28451884 : return std::sqrt(x * x + y * y + z * z);
128 : }
129 :
130 : // square of magnitude of the vector
131 : T getR2() const {
132 2883657 : return x * x + y * y + z * z;
133 : }
134 :
135 : // return the azimuth angle
136 19 : T getPhi() const {
137 : T eps = std::numeric_limits < T > ::min();
138 19 : if ((fabs(x) < eps) && (fabs(y) < eps))
139 : return 0.0;
140 : else
141 18 : return std::atan2(y, x);
142 : }
143 :
144 : // return the zenith angle
145 8 : T getTheta() const {
146 : T eps = std::numeric_limits < T > ::min();
147 8 : if ((fabs(x) < eps) && (fabs(y) < eps) && (fabs(z) < eps))
148 : return 0.0;
149 : else
150 8 : return atan2((T) sqrt(x * x + y * y), z);
151 : }
152 :
153 : // return the unit-vector e_r
154 7201366 : Vector3<T> getUnitVector() const {
155 7201366 : return *this / getR();
156 : }
157 :
158 : // return the unit-vector e_theta
159 1 : Vector3<T> getUnitVectorTheta() const {
160 1 : T theta = getTheta();
161 1 : T phi = getPhi();
162 1 : return Vector3<T>(cos(theta) * cos(phi), cos(theta) * sin(phi),
163 1 : -sin(theta));
164 : }
165 :
166 : // return the unit-vector e_phi
167 4 : Vector3<T> getUnitVectorPhi() const {
168 4 : return Vector3<T>(-sin(getPhi()), cos(getPhi()), 0);
169 : }
170 :
171 : // return the angle [0, pi] between the vectors
172 691021 : T getAngleTo(const Vector3<T> &v) const {
173 691021 : T cosdistance = dot(v) / v.getR() / getR();
174 : // In some directions cosdistance is > 1 on some compilers
175 : // This ensures that the correct result is returned
176 691021 : if (cosdistance >= 1.)
177 : return 0;
178 690876 : else if (cosdistance <= -1.)
179 : return M_PI;
180 : else
181 690876 : return acos(cosdistance);
182 : }
183 :
184 : // return true if the angle between the vectors is smaller than a threshold
185 : bool isParallelTo(const Vector3<T> &v, T maxAngle) const {
186 691017 : return getAngleTo(v) < maxAngle;
187 : }
188 :
189 : // linear distance to a given vector
190 110 : T getDistanceTo(const Vector3<T> &point) const {
191 : Vector3<T> d = *this - point;
192 110 : return d.getR();
193 : }
194 :
195 : // return the component parallel to a second vector
196 : // 0 if the second vector has 0 magnitude
197 2 : Vector3<T> getParallelTo(const Vector3<T> &v) const {
198 : T vmag = v.getR();
199 2 : if (vmag == std::numeric_limits < T > ::min())
200 : return Vector3<T>(0.);
201 : return v * dot(v) / vmag;
202 : }
203 :
204 : // return the component perpendicular to a second vector
205 : // 0 if the second vector has 0 magnitude
206 1 : Vector3<T> getPerpendicularTo(const Vector3<T> &v) const {
207 1 : if (v.getR() == std::numeric_limits < T > ::min())
208 : return Vector3<T>(0.);
209 1 : return (*this) - getParallelTo(v);
210 : }
211 :
212 : // rotate the vector around a given axis by a given angle
213 481006 : Vector3<T> getRotated(const Vector3<T> &axis, T angle) const {
214 : Vector3<T> u = axis;
215 481006 : if (u.getR() != 0.)
216 : u = u / u.getR();
217 481006 : T c = cos(angle);
218 481006 : T s = sin(angle);
219 481006 : Vector3<T> Rx(c + u.x * u.x * (1 - c), u.x * u.y * (1 - c) - u.z * s,
220 481006 : u.x * u.z * (1 - c) + u.y * s);
221 481006 : Vector3<T> Ry(u.y * u.x * (1 - c) + u.z * s, c + u.y * u.y * (1 - c),
222 481006 : u.y * u.z * (1 - c) - u.x * s);
223 481006 : Vector3<T> Rz(u.z * u.x * (1 - c) - u.y * s,
224 481006 : u.z * u.y * (1 - c) + u.x * s, c + u.z * u.z * (1 - c));
225 481006 : return Vector3<T>(dot(Rx), dot(Ry), dot(Rz));
226 : }
227 :
228 : // return vector with values limited to the range [lower, upper]
229 0 : Vector3<T> clip(T lower, T upper) const {
230 : Vector3<T> out;
231 0 : out.x = std::max(lower, std::min(x, upper));
232 0 : out.y = std::max(lower, std::min(y, upper));
233 0 : out.z = std::max(lower, std::min(z, upper));
234 0 : return out;
235 : }
236 :
237 : // return vector with absolute values
238 : Vector3<T> abs() const {
239 0 : return Vector3<T>(std::abs(x), std::abs(y), std::abs(z));
240 : }
241 :
242 : // return vector with floored values
243 6 : Vector3<T> floor() const {
244 6 : return Vector3<T>(std::floor(x), std::floor(y), std::floor(z));
245 : }
246 :
247 : // return vector with ceiled values
248 0 : Vector3<T> ceil() const {
249 0 : return Vector3<T>(std::ceil(x), std::ceil(y), std::ceil(z));
250 : }
251 :
252 : // minimum element
253 : T min() const {
254 4 : return std::min(x, std::min(y, z));
255 : }
256 :
257 : // maximum element
258 : T max() const {
259 4 : return std::max(x, std::max(y, z));
260 : }
261 :
262 : // dot product
263 : T dot(const Vector3<T> &v) const {
264 692580 : return x * v.x + y * v.y + z * v.z;
265 : }
266 :
267 : // cross product
268 : Vector3<T> cross(const Vector3<T> &v) const {
269 2132230 : return Vector3<T>(y * v.z - v.y * z, z * v.x - v.z * x,
270 1652230 : x * v.y - v.x * y);
271 : }
272 :
273 : // returns true if all elements of the two vectors are equal
274 : bool operator ==(const Vector3<T> &v) const {
275 33 : if (x != v.x)
276 : return false;
277 33 : if (y != v.y)
278 : return false;
279 33 : if (z != v.z)
280 : return false;
281 : return true;
282 : }
283 :
284 : Vector3<T> operator +(const Vector3<T> &v) const {
285 4337555 : return Vector3(x + v.x, y + v.y, z + v.z);
286 : }
287 :
288 : Vector3<T> operator +(const T &f) const {
289 0 : return Vector3(x + f, y + f, z + f);
290 : }
291 :
292 : Vector3<T> operator -(const Vector3<T> &v) const {
293 10872055 : return Vector3(x - v.x, y - v.y, z - v.z);
294 : }
295 :
296 : Vector3<T> operator -(const T &f) const {
297 0 : return Vector3(x - f, y - f, z - f);
298 : }
299 :
300 : // element-wise multiplication
301 : Vector3<T> operator *(const Vector3<T> &v) const {
302 20 : return Vector3(x * v.x, y * v.y, z * v.z);
303 : }
304 :
305 : Vector3<T> operator *(T v) const {
306 1184226 : return Vector3(data[0] * v, data[1] * v, data[2] * v);
307 : }
308 :
309 : // element-wise division
310 : Vector3<T> operator /(const Vector3<T> &v) const {
311 100097 : return Vector3(x / v.x, y / v.y, z / v.z);
312 : }
313 :
314 : Vector3<T> operator /(const T &f) const {
315 6242748 : return Vector3(x / f, y / f, z / f);
316 : }
317 :
318 : // element-wise modulo operation
319 0 : Vector3<T> operator %(const Vector3<T> &v) const {
320 0 : return Vector3(fmod(x, v.x), fmod(y, v.y), fmod(z, v.z));
321 : }
322 :
323 0 : Vector3<T> operator %(const T &f) const {
324 0 : return Vector3(fmod(x, f), fmod(y, f), fmod(z, f));
325 : }
326 :
327 : Vector3<T> &operator -=(const Vector3<T> &v) {
328 0 : data[0] -= v.x;
329 0 : data[1] -= v.y;
330 0 : data[2] -= v.z;
331 : return *this;
332 : }
333 :
334 : Vector3<T> &operator -=(const T &f) {
335 0 : data[0] -= f;
336 0 : data[1] -= f;
337 0 : data[2] -= f;
338 : return *this;
339 : }
340 :
341 : Vector3<T> &operator +=(const Vector3<T> &v) {
342 20664839 : data[0] += v.x;
343 20664839 : data[1] += v.y;
344 20544720 : data[2] += v.z;
345 : return *this;
346 : }
347 :
348 : Vector3<T> &operator +=(const T &f) {
349 0 : data[0] += f;
350 0 : data[1] += f;
351 0 : data[2] += f;
352 : return *this;
353 : }
354 :
355 : // element-wise multiplication
356 : Vector3<T> &operator *=(const Vector3<T> &v) {
357 0 : data[0] *= v.x;
358 0 : data[1] *= v.y;
359 0 : data[2] *= v.z;
360 : return *this;
361 : }
362 :
363 : Vector3<T> &operator *=(const T &f) {
364 3312562 : data[0] *= f;
365 3312562 : data[1] *= f;
366 3312562 : data[2] *= f;
367 : return *this;
368 : }
369 :
370 : // element-wise division
371 : Vector3<T> &operator /=(const Vector3<T> &v) {
372 0 : data[0] /= v.x;
373 0 : data[1] /= v.y;
374 0 : data[2] /= v.z;
375 : return *this;
376 : }
377 :
378 : Vector3<T> &operator /=(const T &f) {
379 691019 : data[0] /= f;
380 691019 : data[1] /= f;
381 691019 : data[2] /= f;
382 : return *this;
383 : }
384 :
385 : // element-wise modulo operation
386 0 : Vector3<T> &operator %=(const Vector3<T> &v) {
387 0 : data[0] = fmod(x, v.x);
388 0 : data[1] = fmod(y, v.y);
389 0 : data[2] = fmod(z, v.z);
390 0 : return *this;
391 : }
392 :
393 1 : Vector3<T> &operator %=(const T &f) {
394 1 : data[0] = fmod(x, f);
395 1 : data[1] = fmod(y, f);
396 1 : data[2] = fmod(z, f);
397 1 : return *this;
398 : }
399 :
400 : Vector3<T> &operator =(const Vector3<T> &v) {
401 64166297 : data[0] = v.x;
402 64166309 : data[1] = v.y;
403 63680684 : data[2] = v.z;
404 0 : return *this;
405 : }
406 :
407 : //Vector3<T> &operator =(Vector3<T> &&v) noexcept {
408 : // data[0] = v.data[0];
409 : // data[1] = v.data[1];
410 : // data[2] = v.data[2];
411 : // return *this;
412 : //}
413 :
414 : Vector3<T> &operator =(const T &f) {
415 : data[0] = f;
416 : data[1] = f;
417 : data[2] = f;
418 : return *this;
419 : }
420 :
421 2 : const std::string getDescription() {
422 : char buffer[256];
423 2 : snprintf(buffer, 256, "Vector(%.6G, %.6G, %.6G)", data[0], data[1], data[2]);
424 2 : return buffer;
425 : }
426 :
427 : // ----------------------------
428 : // Python numpy interface
429 : // ----------------------------
430 : #ifdef CRPROPA_HAVE_PYTHON
431 0 : PyObject* __array__() {
432 0 : npy_intp dims[1] = {3};
433 : PyObject *array;
434 : // type handling
435 0 : if (typeid(T) == typeid(double))
436 0 : array = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, (void *)data);
437 0 : else if (typeid(T) == typeid(float))
438 0 : array = PyArray_SimpleNewFromData(1, dims, NPY_FLOAT, (void *)data);
439 : else {
440 0 : PyErr_SetString(PyExc_TypeError, "Unsupported type for Vector3");
441 0 : return NULL;
442 : }
443 : return array;
444 : }
445 :
446 : size_t __len__() {
447 : return 3;
448 : }
449 :
450 : // initialize the vector from a numpy array
451 4 : Vector3(PyObject* inputObject) {
452 : // check if input is a numpy array
453 4 : if (!PyArray_Check(inputObject)) {
454 0 : throw std::runtime_error("Expected a numpy array as input");
455 : }
456 :
457 : // convert to PyArrayObject
458 4 : PyArrayObject *array = (PyArrayObject*)PyArray_FROM_O(inputObject);
459 :
460 : // chek dimensions
461 4 : if (PyArray_NDIM(array) != 1 || PyArray_DIM(array, 0) != 3) {
462 0 : throw std::runtime_error("Expected a 1D array of length 3");
463 : }
464 :
465 : // handle different numpy dtypes
466 4 : if (PyArray_TYPE(array) == NPY_DOUBLE) {
467 : // convert pyarray data
468 : double *dataPtr = (double *)PyArray_DATA(array);
469 1 : data[0] = dataPtr[0];
470 1 : data[1] = dataPtr[1];
471 1 : data[2] = dataPtr[2];
472 3 : } else if (PyArray_TYPE(array) == NPY_FLOAT) {
473 : // convert pyarray data
474 : float *dataPtr = (float *)PyArray_DATA(array);
475 1 : data[0] = dataPtr[0];
476 1 : data[1] = dataPtr[1];
477 1 : data[2] = dataPtr[2];
478 2 : } else if (PyArray_TYPE(array) == NPY_INT32) {
479 : // convert pyarray data
480 : int *dataPtr = (int *)PyArray_DATA(array);
481 1 : data[0] = static_cast<T>(dataPtr[0]);
482 1 : data[1] = static_cast<T>(dataPtr[1]);
483 1 : data[2] = static_cast<T>(dataPtr[2]);
484 1 : } else if (PyArray_TYPE(array) == NPY_INT64) {
485 : // convert pyarray data
486 : long long *dataPtr = (long long *)PyArray_DATA(array);
487 1 : data[0] = static_cast<T>(dataPtr[0]);
488 1 : data[1] = static_cast<T>(dataPtr[1]);
489 1 : data[2] = static_cast<T>(dataPtr[2]);
490 : } else {
491 0 : throw std::runtime_error("Unsupported numpy array dtype");
492 : }
493 :
494 : // free the reference to the numpy array
495 : Py_DECREF(array);
496 4 : }
497 : #endif // CRPROPA_HAVE_PYTHON
498 :
499 : };
500 :
501 : #ifndef SWIG
502 : template<typename T>
503 66 : inline std::ostream &operator <<(std::ostream &out, const Vector3<T> &v) {
504 198 : out << v.x << " " << v.y << " " << v.z;
505 66 : return out;
506 : }
507 :
508 : template<typename T>
509 : inline std::istream &operator >>(std::istream &in, Vector3<T> &v) {
510 : in >> v.x >> v.y >> v.z;
511 : return in;
512 : }
513 : #endif
514 :
515 : template<typename T>
516 : inline Vector3<T> operator *(T f, const Vector3<T> &v) {
517 3634379 : return Vector3<T>(v.x * f, v.y * f, v.z * f);
518 : }
519 :
520 : typedef Vector3<double> Vector3d;
521 : typedef Vector3<float> Vector3f;
522 :
523 : /** @}*/
524 : } // namespace crpropa
525 :
526 : #endif // CRPROPA_VECTOR3_H
527 :
528 :
|
