1 /* This file is part of the hkl library.
3 * The hkl library is free software: you can redistribute it and/or modify
4 * it under the terms of the GNU General Public License as published by
5 * the Free Software Foundation, either version 3 of the License, or
6 * (at your option) any later version.
8 * The hkl library is distributed in the hope that it will be useful,
9 * but WITHOUT ANY WARRANTY; without even the implied warranty of
10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 * GNU General Public License for more details.
13 * You should have received a copy of the GNU General Public License
14 * along with the hkl library. If not, see <http://www.gnu.org/licenses/>.
16 * Copyright (C) 2003-2015 Synchrotron SOLEIL
17 * L'Orme des Merisiers Saint-Aubin
18 * BP 48 91192 GIF-sur-YVETTE CEDEX
20 * Authors: Picca Frédéric-Emmanuel <picca@synchrotron-soleil.fr>
22 #include <math.h> // for fabs, acos, cos, sin, sqrt, etc
23 #include <stdio.h> // for fprintf, FILE
24 #include <stdlib.h> // for rand, RAND_MAX, free
25 #include <string.h> // for memcpy
26 #include "hkl-macros-private.h" // for HKL_MALLOC
27 #include "hkl-matrix-private.h" // for _HklMatrix
28 #include "hkl-quaternion-private.h" // for _HklQuaternion
29 #include "hkl-vector-private.h" // for HklVector, HklQuaternion
30 #include "hkl.h" // for HklMatrix, HKL_EPSILON, etc
33 * hkl_vector_dup: (skip)
34 * @self: the HklVector to copy
38 * Returns: A copy of self which need to be free using hkl_vector_free
40 HklVector
* hkl_vector_dup (const HklVector
* self
) {
43 dup
= HKL_MALLOC(HklVector
);
44 memcpy(dup
, self
, sizeof (HklVector
));
49 * hkl_vector_free: (skip)
52 * delete an HklVector struct
54 void hkl_vector_free (HklVector
* self
) {
60 * @self: the #HklVector to initialize.
61 * @x: the first coordinate value
62 * @y: the second coordinate value
63 * @z: the third coordinate value
65 * initialize an #HklVector
67 void hkl_vector_init(HklVector
*self
, double x
, double y
, double z
)
75 * hkl_vector_fprintf: (skip)
76 * @file: the stream to print into
77 * @self: the #HklVector to print.
79 * print an #HklVector into a stream
81 void hkl_vector_fprintf(FILE *file
, const HklVector
*self
)
83 fprintf(file
, "|%f, %f, %f|", self
->data
[0], self
->data
[1], self
->data
[2]);
87 * hkl_vector_cmp: (skip)
88 * @self: the first vector
89 * @vector: th vector to compare with
91 * compare two #HklVector. this comparison use HKL_EPSILON
92 * to do the comparison.
94 * Returns: FALSE if both are equals, TRUE otherwise.
96 int hkl_vector_cmp(const HklVector
*self
, const HklVector
*vector
)
101 if ( fabs(self
->data
[i
] - vector
->data
[i
]) > HKL_EPSILON
)
107 * hkl_vector_is_opposite: (skip)
111 * Check if two vectors are oposite.
113 * Returns: TRUE is vector are oposite vectors.
115 int hkl_vector_is_opposite(const HklVector
*self
, const HklVector
*vector
)
120 if ( fabs(self
->data
[i
] + vector
->data
[i
]) > HKL_EPSILON
)
126 * hkl_vector_add_vector: (skip)
127 * @self: the modified #HklVector
128 * @vector: the #hklvector to add
130 * add an #HklVector to another one.
132 void hkl_vector_add_vector(HklVector
*self
, const HklVector
*vector
)
136 self
->data
[i
] += vector
->data
[i
];
140 * hkl_vector_minus_vector: (skip)
141 * @self: the modified #HklVector
142 * @vector: the #hklvector to substract
144 * substract an #HklVector to another one.
146 void hkl_vector_minus_vector(HklVector
*self
, const HklVector
*vector
)
150 self
->data
[i
] -= vector
->data
[i
];
154 * hkl_vector_div_double: (skip)
155 * @self: the #HklVector to divide.
156 * @d: constant use to divide the #HklVector
158 * divide an #HklVector by constant.
160 void hkl_vector_div_double(HklVector
*self
, const double d
)
168 * hkl_vector_times_double: (skip)
169 * @self: the #HklVector to modify
170 * @d: the multiply factor
172 * multiply an #HklVector by a constant value.
174 void hkl_vector_times_double(HklVector
*self
, const double d
)
182 * hkl_vector_times_vector: (skip)
183 * @self: the #HklVector to modify
184 * @vector: the #HklVector use to modify the first one
186 * multiply an #HklVector by another one. This method multiply
187 * coordinate by coordinate.
189 void hkl_vector_times_vector(HklVector
*self
, const HklVector
*vector
)
193 self
->data
[i
] *= vector
->data
[i
];
197 * hkl_vector_times_matrix: (skip)
198 * @self: the #HklVector to multiply
199 * @m: the #HklMatrix use to multiply the #HklVector
201 * multiply an #HklVector by an #HklMatrix.
204 void hkl_vector_times_matrix(HklVector
*self
, const HklMatrix
*m
)
209 self
->data
[0] = tmp
.data
[0] *m
->data
[0][0] + tmp
.data
[1] *m
->data
[1][0] + tmp
.data
[2] *m
->data
[2][0];
210 self
->data
[1] = tmp
.data
[0] *m
->data
[0][1] + tmp
.data
[1] *m
->data
[1][1] + tmp
.data
[2] *m
->data
[2][1];
211 self
->data
[2] = tmp
.data
[0] *m
->data
[0][2] + tmp
.data
[1] *m
->data
[1][2] + tmp
.data
[2] *m
->data
[2][2];
215 * hkl_vector_sum: (skip)
216 * @self: the #HklVector to sum.
218 * compute the #HklVector sum of all its elements.
220 * Returns: the sum of all elements.
222 double hkl_vector_sum(const HklVector
*self
)
224 return self
->data
[0] + self
->data
[1] + self
->data
[2];
228 * hkl_vector_scalar_product: (skip)
229 * @self: the first #HklVector
230 * @vector: the second #HklVector
232 * compute the scalar product of two #HklVector
234 * Returns: the scalar product.
236 double hkl_vector_scalar_product(const HklVector
*self
, const HklVector
*vector
)
242 scalar
+= self
->data
[i
] *vector
->data
[i
];
247 * hkl_vector_vectorial_product: (skip)
248 * @self: the first #HklVector (modify)
249 * @vector: the second #HklVector
251 * compute the vectorial product of two vectors
253 void hkl_vector_vectorial_product(HklVector
*self
, const HklVector
*vector
)
258 self
->data
[0] = tmp
.data
[1] * vector
->data
[2] - tmp
.data
[2] * vector
->data
[1];
259 self
->data
[1] = tmp
.data
[2] * vector
->data
[0] - tmp
.data
[0] * vector
->data
[2];
260 self
->data
[2] = tmp
.data
[0] * vector
->data
[1] - tmp
.data
[1] * vector
->data
[0];
265 * hkl_vector_angle: (skip)
266 * @self: the fist #HklVector
267 * @vector: the second #HklVector
269 * compute the angles beetween two #HklVector
271 * Returns: the return value is in beetween [0, pi]
273 double hkl_vector_angle(const HklVector
*self
, const HklVector
*vector
)
281 norm_self
= hkl_vector_norm2(self
);
282 norm_vector
= hkl_vector_norm2(vector
);
284 if (norm_self
< HKL_EPSILON
|| norm_vector
< HKL_EPSILON
)
287 norm
= norm_self
* norm_vector
;
289 cos_angle
= hkl_vector_scalar_product(self
, vector
) / norm
;
291 /* problem with round */
295 if (cos_angle
<= -1 )
298 angle
= acos(cos_angle
);
303 * hkl_vector_oriented_angle: (skip)
304 * @self: the first #HklVector
305 * @vector: the second #HklVector
306 * @ref: the reference #HklVector
308 * compute the angles beetween two #HklVector and use
309 * a reference #HklVector to orientate the space. That's
310 * way the return value can be in beetween [-pi, pi].
311 * the (self, vector, ref) is a right oriented base.
313 * Returns: the angles [-pi, pi]
315 double hkl_vector_oriented_angle(const HklVector
*self
,
316 const HklVector
*vector
,
317 const HklVector
*ref
)
323 angle
= hkl_vector_angle(self
, vector
);
325 hkl_vector_vectorial_product(&tmp
, vector
);
326 hkl_vector_normalize(&tmp
);
328 hkl_vector_normalize(&ref_u
);
329 if (hkl_vector_is_opposite(&tmp
, &ref_u
))
334 * hkl_vector_oriented_angle_points: (skip)
335 * @self: the first point
336 * @p2: the second point
337 * @p3: the third point
338 * @ref: the reference #HklVector
340 * compute the angles beetween three points (p1, p2, p3) and use
341 * a reference #HklVector to orientate the space. That's
342 * way the return value can be in beetween [-pi, pi].
343 * the (self, vector, ref) is a right oriented base.
345 * Returns: the angles [-pi, pi]
347 double hkl_vector_oriented_angle_points(const HklVector
*self
,
350 const HklVector
*ref
)
357 hkl_vector_minus_vector(&v1
, p2
);
358 hkl_vector_minus_vector(&v2
, p2
);
359 return hkl_vector_oriented_angle(&v1
, &v2
, ref
);
363 * hkl_vector_normalize: (skip)
364 * @self: the #HklVector to normalize
366 * normalize a hkl_vector
368 * Returns: TRUE if the #HklVector can be normalized, FALSE otherwise
370 int hkl_vector_normalize(HklVector
*self
)
372 double norm
= hkl_vector_norm2(self
);
373 if ( norm
<= HKL_EPSILON
)
376 hkl_vector_div_double(self
, norm
);
382 * hkl_vector_is_colinear: (skip)
383 * @self: the first #HklVector
384 * @vector: the second #HklVector
386 * check if two #HklVector are colinears
388 * Returns: TRUE if both are colinear.
390 int hkl_vector_is_colinear(const HklVector
*self
, const HklVector
*vector
)
393 HklVector tmp
= *self
;
395 hkl_vector_vectorial_product(&tmp
, vector
);
396 if (hkl_vector_norm2(&tmp
) < HKL_EPSILON
)
404 * hkl_vector_randomize: (skip)
405 * @self: the #HklVector to randomize
407 * initialize a vector with random values.
408 * coordinates range [-1, 1]
410 void hkl_vector_randomize(HklVector
*self
)
412 self
->data
[0] = -1 + 2 *rand()/(RAND_MAX
+1.0);
413 self
->data
[1] = -1 + 2 *rand()/(RAND_MAX
+1.0);
414 self
->data
[2] = -1 + 2 *rand()/(RAND_MAX
+1.0);
418 * hkl_vector_randomize_vector: (skip)
419 * @self: the #HklVector to randomize
420 * @vector: the #HklVector result to avoid
422 * randomize an #HklVector an be sure that it is not equal
423 * to the #HklVector vector.
425 void hkl_vector_randomize_vector(HklVector
*self
, const HklVector
*vector
)
428 hkl_vector_randomize(self
);
429 while (!hkl_vector_cmp(self
, vector
));
433 * hkl_vector_randomize_vector_vector: (skip)
434 * @self: the #HklVector to randomize
435 * @vector1: the first #HklVector solution to avoid
436 * @vector2: the second #HklVector solution to avoid
438 * randomize an #HklVector an be sure that it is not equal
439 * to the #HklVector vector1 and vector2.
442 void hkl_vector_randomize_vector_vector(HklVector
*self
,
443 const HklVector
*vector1
,
444 const HklVector
*vector2
)
447 hkl_vector_randomize(self
);
448 while (!hkl_vector_cmp(self
, vector1
) || !hkl_vector_cmp(self
, vector2
));
452 * hkl_vector_rotated_around_vector: (skip)
453 * @self: the #HklVector to rotate
454 * @axe: the axe of rotation
455 * @angle: the angle of the rotation
457 * rotate a vector around another one with a given angle.
459 void hkl_vector_rotated_around_vector(HklVector
*self
,
460 const HklVector
*axe
, double angle
)
462 double c
= cos(angle
);
463 double s
= sin(angle
);
468 hkl_vector_normalize(&axe_n
);
472 self
->data
[0] = (c
+ (1 - c
) * axe_n
.data
[0] * axe_n
.data
[0]) * tmp
.data
[0];
473 self
->data
[0] += ((1 - c
) * axe_n
.data
[0] * axe_n
.data
[1] - axe_n
.data
[2] * s
) * tmp
.data
[1];
474 self
->data
[0] += ((1 - c
) * axe_n
.data
[0] * axe_n
.data
[2] + axe_n
.data
[1] * s
) * tmp
.data
[2];
476 self
->data
[1] = ((1 - c
) * axe_n
.data
[0] * axe_n
.data
[1] + axe_n
.data
[2] * s
) * tmp
.data
[0];
477 self
->data
[1] += (c
+ (1 - c
) * axe_n
.data
[1] * axe_n
.data
[1]) * tmp
.data
[1];
478 self
->data
[1] += ((1 - c
) * axe_n
.data
[1] * axe_n
.data
[2] - axe_n
.data
[0] * s
) * tmp
.data
[2];
480 self
->data
[2] = ((1 - c
) * axe_n
.data
[0] * axe_n
.data
[2] - axe_n
.data
[1] * s
) * tmp
.data
[0];
481 self
->data
[2] += ((1 - c
) * axe_n
.data
[1] * axe_n
.data
[2] + axe_n
.data
[0] * s
) * tmp
.data
[1];
482 self
->data
[2] += (c
+ (1 - c
) * axe_n
.data
[2] * axe_n
.data
[2]) * tmp
.data
[2];
486 * hkl_vector_norm2: (skip)
487 * @self: the #hklvector use to compute the norm2
489 * compute the norm2 of an #HklVector
491 * Returns: the sqrt(|v|)
493 double hkl_vector_norm2(const HklVector
*self
)
495 return sqrt(self
->data
[0] * self
->data
[0]
496 + self
->data
[1] * self
->data
[1]
497 + self
->data
[2] * self
->data
[2]);
501 * hkl_vector_rotated_quaternion: (skip)
502 * @self: the #HklVector to rotate
503 * @qr: the #HklQuaternion use to rotate the vector
505 * rotate an #HklVector using an #HklQuaternion.
507 void hkl_vector_rotated_quaternion(HklVector
*self
, const HklQuaternion
*qr
)
509 double v1
= self
->data
[0];
510 double v2
= self
->data
[1];
511 double v3
= self
->data
[2];
512 double a
= qr
->data
[0];
513 double b
= qr
->data
[1];
514 double c
= qr
->data
[2];
515 double d
= qr
->data
[3];
527 self
->data
[0] = 2*( (t8
+ t10
)*v1
+ (t6
- t4
)*v2
+ (t3
+ t7
)*v3
) + v1
;
528 self
->data
[1] = 2*( (t4
+ t6
)*v1
+ (t5
+ t10
)*v2
+ (t9
- t2
)*v3
) + v2
;
529 self
->data
[2] = 2*( (t7
- t3
)*v1
+ (t2
+ t9
)*v2
+ (t5
+ t8
)*v3
) + v3
;
533 * hkl_vector_rotated_around_line: (skip)
534 * @self: the point to rotate around a line
535 * @angle: the angle of the rotation
536 * @c1: the fist point of the line
537 * @c2: the second point of the line
539 * This method rotate a point around a line defined by two points
540 * of a certain amount of angle. The rotation is right handed.
541 * this mean that c2 - c1 gives the direction of the rotation.
543 void hkl_vector_rotated_around_line(HklVector
*self
, double angle
,
544 const HklVector
*c1
, const HklVector
*c2
)
548 if (!self
|| !c1
|| !c2
|| fabs(angle
) < HKL_EPSILON
)
552 hkl_vector_minus_vector(&axis
, c1
);
553 /* the c2 - c1 vector must be non null */
555 hkl_vector_minus_vector(self
, c1
);
556 hkl_vector_rotated_around_vector(self
, &axis
, angle
);
557 hkl_vector_add_vector(self
, c1
);
561 * hkl_vector_is_null: (skip)
562 * @self: the #hklvector to check
564 * check if all the coordinates of an #HklVector are null.
566 * Returns: HKl_TRUE if all |elements| are below HKL_EPSILON, HKl_FALSE otherwise
570 int hkl_vector_is_null(const HklVector
*self
)
574 if ( fabs(self
->data
[i
]) > HKL_EPSILON
)
580 * hkl_vector_project_on_plan: (skip)
581 * @self: the vector to project
582 * @normal: the normal of the plane.
584 * project an #HklVector on a plan of normal which contain
585 * the origin [0, 0, 0]
588 void hkl_vector_project_on_plan(HklVector
*self
,
589 const HklVector
*normal
)
597 hkl_vector_normalize(&tmp
);
598 hkl_vector_times_double(&tmp
, hkl_vector_scalar_product(self
, &tmp
));
599 hkl_vector_minus_vector(self
, &tmp
);
603 * hkl_vector_project_on_plan_with_point: (skip)
604 * @self: the vector to project (modify)
605 * @normal: the normal of the plane.
606 * @point: a point of the plan.
608 * project an #HklVector on a plan of normal #normal which contain #point.
610 void hkl_vector_project_on_plan_with_point(HklVector
*self
,
611 const HklVector
*normal
,
612 const HklVector
*point
)
617 if(!self
|| !normal
|| !point
)
621 hkl_vector_normalize(&tmp
);
622 d1
= hkl_vector_scalar_product(self
, &tmp
);
623 d2
= hkl_vector_scalar_product(point
, &tmp
);
624 hkl_vector_times_double(&tmp
, d1
- d2
);
625 hkl_vector_minus_vector(self
, &tmp
);