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-2010 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>
25 #include <gsl/gsl_math.h>
27 #include <hkl/hkl-vector.h>
28 #include <hkl/hkl-matrix.h>
29 #include <hkl/hkl-quaternion.h>
33 * @self: the #HklVector to initialize.
34 * @x: the first coordinate value
35 * @y: the second coordinate value
36 * @z: the third coordinate value
38 * initialize an #HklVector
40 void hkl_vector_init(HklVector
*self
, double x
, double y
, double z
)
49 * @file: the stream to print into
50 * @self: the #HklVector to print.
52 * print an #HklVector into a stream
54 void hkl_vector_fprintf(FILE *file
, const HklVector
*self
)
56 fprintf(file
, "|%f, %f, %f|", self
->data
[0], self
->data
[1], self
->data
[2]);
61 * @self: the first vector
62 * @vector: th vector to compare with
64 * compare two #HklVector. this comparison use HKL_EPSILON
65 * to do the comparison.
67 * Returns: HKL_FALSE if both are equals, HKL_TRUE otherwise.
69 int hkl_vector_cmp(const HklVector
*self
, const HklVector
*vector
)
74 if ( fabs(self
->data
[i
] - vector
->data
[i
]) > HKL_EPSILON
)
80 int hkl_vector_is_opposite(const HklVector
*self
, const HklVector
*vector
)
85 if ( fabs(self
->data
[i
] + vector
->data
[i
]) > HKL_EPSILON
)
91 * hkl_vector_add_vector:
92 * @self: the modified #HklVector
93 * @vector: the #hklvector to add
95 * add an #HklVector to another one.
97 void hkl_vector_add_vector(HklVector
*self
, const HklVector
*vector
)
101 self
->data
[i
] += vector
->data
[i
];
105 * hkl_vector_minus_vector:
106 * @self: the modified #HklVector
107 * @vector: the #hklvector to substract
109 * substract an #HklVector to another one.
111 void hkl_vector_minus_vector(HklVector
*self
, const HklVector
*vector
)
115 self
->data
[i
] -= vector
->data
[i
];
119 * hkl_vector_div_double:
120 * @self: the #HklVector to divide.
121 * @d: constant use to divide the #HklVector
123 * divide an #HklVector by constant.
125 void hkl_vector_div_double(HklVector
*self
, const double d
)
133 * hkl_vector_times_double:
134 * @self: the #HklVector to modify
135 * @d: the multiply factor
137 * multiply an #HklVector by a constant value.
139 void hkl_vector_times_double(HklVector
*self
, const double d
)
147 * hkl_vector_times_vector:
148 * @self: the #HklVector to modify
149 * @vector: the #HklVector use to modify the first one
151 * multiply an #HklVector by another one. This method multiply
152 * coordinate by coordinate.
154 void hkl_vector_times_vector(HklVector
*self
, const HklVector
*vector
)
158 self
->data
[i
] *= vector
->data
[i
];
162 * hkl_vector_times_smatrix:
163 * @self: the #HklVector to multiply
164 * @m: the #HklMatrix use to multiply the #HklVector
166 * multiply an #HklVector by an #HklMatrix.
169 void hkl_vector_times_matrix(HklVector
*self
, const HklMatrix
*m
)
174 self
->data
[0] = tmp
.data
[0] *m
->data
[0][0] + tmp
.data
[1] *m
->data
[1][0] + tmp
.data
[2] *m
->data
[2][0];
175 self
->data
[1] = tmp
.data
[0] *m
->data
[0][1] + tmp
.data
[1] *m
->data
[1][1] + tmp
.data
[2] *m
->data
[2][1];
176 self
->data
[2] = tmp
.data
[0] *m
->data
[0][2] + tmp
.data
[1] *m
->data
[1][2] + tmp
.data
[2] *m
->data
[2][2];
181 * @self: the #HklVector to sum.
183 * compute the #HklVector sum of all its elements.
185 * Returns: the sum of all elements.
187 double hkl_vector_sum(const HklVector
*self
)
189 return self
->data
[0] + self
->data
[1] + self
->data
[2];
193 * hkl_vector_scalar_product:
194 * @self: the first #HklVector
195 * @vector: the second #HklVector
197 * compute the scalar product of two #HklVector
199 * Returns: the scalar product.
201 double hkl_vector_scalar_product(const HklVector
*self
, const HklVector
*vector
)
207 scalar
+= self
->data
[i
] *vector
->data
[i
];
212 * hkl_vector_vectorial_product:
213 * @self: the first #HklVector (modify)
214 * @vector: the second #HklVector
216 * compute the vectorial product of two vectors
218 void hkl_vector_vectorial_product(HklVector
*self
, const HklVector
*vector
)
223 self
->data
[0] = tmp
.data
[1] * vector
->data
[2] - tmp
.data
[2] * vector
->data
[1];
224 self
->data
[1] = tmp
.data
[2] * vector
->data
[0] - tmp
.data
[0] * vector
->data
[2];
225 self
->data
[2] = tmp
.data
[0] * vector
->data
[1] - tmp
.data
[1] * vector
->data
[0];
231 * @self: the fist #HklVector
232 * @vector: the second #HklVector
234 * compute the angles beetween two #HklVector
236 * Returns: the return value is in beetween [0, pi]
238 double hkl_vector_angle(const HklVector
*self
, const HklVector
*vector
)
246 norm_self
= hkl_vector_norm2(self
);
247 norm_vector
= hkl_vector_norm2(vector
);
249 if (norm_self
< HKL_EPSILON
|| norm_vector
< HKL_EPSILON
)
252 norm
= norm_self
* norm_vector
;
254 cos_angle
= hkl_vector_scalar_product(self
, vector
) / norm
;
256 /* problem with round */
260 if (cos_angle
<= -1 )
263 angle
= acos(cos_angle
);
268 * hkl_vector_oriented_angle:
269 * @self: the first #HklVector
270 * @vector: the second #HklVector
271 * @ref: the reference #HklVector
273 * compute the angles beetween two #HklVector and use
274 * a reference #HklVector to orientate the space. That's
275 * way the return value can be in beetween [-pi, pi].
276 * the (self, vector, ref) is a right oriented base.
278 * Returns: the angles [-pi, pi]
280 double hkl_vector_oriented_angle(const HklVector
*self
,
281 const HklVector
*vector
,
282 const HklVector
*ref
)
288 angle
= hkl_vector_angle(self
, vector
);
290 hkl_vector_vectorial_product(&tmp
, vector
);
291 hkl_vector_normalize(&tmp
);
293 hkl_vector_normalize(&ref_u
);
294 if (hkl_vector_is_opposite(&tmp
, &ref_u
))
300 * hkl_vector_normalize:
301 * @self: the #HklVector to normalize
303 * normalize a hkl_vector
305 * Returns: HKL_SUCCESS if the #HklVector can be normalized, HKL_FAIL otherwise
307 int hkl_vector_normalize(HklVector
*self
)
309 int status
= HKL_FAIL
;
311 double norm
= hkl_vector_norm2(self
);
312 if ( norm
> HKL_EPSILON
)
314 hkl_vector_div_double(self
, norm
);
315 status
= HKL_SUCCESS
;
322 * hkl_vector_is_colinear:
323 * @self: the first #HklVector
324 * @vector: the second #HklVector
326 * check if two #HklVector are colinears
328 * Returns: HKL_TRUE if both are colinear.
330 int hkl_vector_is_colinear(const HklVector
*self
, const HklVector
*vector
)
333 HklVector tmp
= *self
;
335 hkl_vector_vectorial_product(&tmp
, vector
);
336 if (hkl_vector_norm2(&tmp
) < HKL_EPSILON
)
344 * hkl_vector_randomize:
345 * @self: the #HklVector to randomize
347 * initialize a vector with random values.
348 * coordinates range [-1, 1]
350 void hkl_vector_randomize(HklVector
*self
)
352 self
->data
[0] = -1 + 2 *rand()/(RAND_MAX
+1.0);
353 self
->data
[1] = -1 + 2 *rand()/(RAND_MAX
+1.0);
354 self
->data
[2] = -1 + 2 *rand()/(RAND_MAX
+1.0);
358 * hkl_vector_randomize_vector:
359 * @self: the #HklVector to randomize
360 * @vector: the #HklVector result to avoid
362 * randomize an #HklVector an be sure that it is not equal
363 * to the #HklVector vector.
365 void hkl_vector_randomize_vector(HklVector
*self
, const HklVector
*vector
)
368 hkl_vector_randomize(self
);
369 while (!hkl_vector_cmp(self
, vector
));
373 * hkl_vector_randomize_vector_vector:
374 * @self: the #HklVector to randomize
375 * @vector1: the first #HklVector solution to avoid
376 * @vector2: the second #HklVector solution to avoid
378 * randomize an #HklVector an be sure that it is not equal
379 * to the #HklVector vector1 and vector2.
382 void hkl_vector_randomize_vector_vector(HklVector
*self
,
383 const HklVector
*vector1
,
384 const HklVector
*vector2
)
387 hkl_vector_randomize(self
);
388 while (!hkl_vector_cmp(self
, vector1
) || !hkl_vector_cmp(self
, vector2
));
392 * hkl_vector_rotated_around_vector:
393 * @self: the #HklVector to rotate
394 * @axe: the axe of rotation
395 * @angle: the angle of the rotation
397 * rotate a vector around another one with a given angle.
399 void hkl_vector_rotated_around_vector(HklVector
*self
,
400 const HklVector
*axe
, double angle
)
402 double c
= cos(angle
);
403 double s
= sin(angle
);
408 hkl_vector_normalize(&axe_n
);
412 self
->data
[0] = (c
+ (1 - c
) * axe_n
.data
[0] * axe_n
.data
[0]) * tmp
.data
[0];
413 self
->data
[0] += ((1 - c
) * axe_n
.data
[0] * axe_n
.data
[1] - axe_n
.data
[2] * s
) * tmp
.data
[1];
414 self
->data
[0] += ((1 - c
) * axe_n
.data
[0] * axe_n
.data
[2] + axe_n
.data
[1] * s
) * tmp
.data
[2];
416 self
->data
[1] = ((1 - c
) * axe_n
.data
[0] * axe_n
.data
[1] + axe_n
.data
[2] * s
) * tmp
.data
[0];
417 self
->data
[1] += (c
+ (1 - c
) * axe_n
.data
[1] * axe_n
.data
[1]) * tmp
.data
[1];
418 self
->data
[1] += ((1 - c
) * axe_n
.data
[1] * axe_n
.data
[2] - axe_n
.data
[0] * s
) * tmp
.data
[2];
420 self
->data
[2] = ((1 - c
) * axe_n
.data
[0] * axe_n
.data
[2] - axe_n
.data
[1] * s
) * tmp
.data
[0];
421 self
->data
[2] += ((1 - c
) * axe_n
.data
[1] * axe_n
.data
[2] + axe_n
.data
[0] * s
) * tmp
.data
[1];
422 self
->data
[2] += (c
+ (1 - c
) * axe_n
.data
[2] * axe_n
.data
[2]) * tmp
.data
[2];
427 * @self: the #hklvector use to compute the norm2
429 * compute the norm2 of an #HklVector
431 * Returns: the sqrt(|v|)
433 double hkl_vector_norm2(const HklVector
*self
)
435 return sqrt(self
->data
[0] * self
->data
[0]
436 + self
->data
[1] * self
->data
[1]
437 + self
->data
[2] * self
->data
[2]);
441 * hkl_vector_rotated_quaternion:
442 * @self: the #HklVector to rotate
443 * @qr: the #HklQuaternion use to rotate the vector
445 * rotate an #HklVector using an #HklQuaternion.
447 void hkl_vector_rotated_quaternion(HklVector
*self
, const HklQuaternion
*qr
)
449 double v1
= self
->data
[0];
450 double v2
= self
->data
[1];
451 double v3
= self
->data
[2];
452 double a
= qr
->data
[0];
453 double b
= qr
->data
[1];
454 double c
= qr
->data
[2];
455 double d
= qr
->data
[3];
467 self
->data
[0] = 2*( (t8
+ t10
)*v1
+ (t6
- t4
)*v2
+ (t3
+ t7
)*v3
) + v1
;
468 self
->data
[1] = 2*( (t4
+ t6
)*v1
+ (t5
+ t10
)*v2
+ (t9
- t2
)*v3
) + v2
;
469 self
->data
[2] = 2*( (t7
- t3
)*v1
+ (t2
+ t9
)*v2
+ (t5
+ t8
)*v3
) + v3
;
473 * hkl_vector_is_null:
474 * @self: the #hklvector to check
476 * check if all the coordinates of an #HklVector are null.
478 * Returns: HKl_TRUE if all |elements| are below HKL_EPSILON, HKl_FALSE otherwise
482 int hkl_vector_is_null(const HklVector
*self
)
486 if ( fabs(self
->data
[i
]) > HKL_EPSILON
)
492 * hkl_vector_project_on_plan:
493 * @self: the vector to project (modify)
494 * @plan: the normal of the plane.
496 * project an #HklVector on a plan.
500 void hkl_vector_project_on_plan(HklVector
*self
,
501 const HklVector
*plan
)
506 hkl_vector_normalize(&tmp
);
507 hkl_vector_times_double(&tmp
, hkl_vector_scalar_product(self
, &tmp
));
508 hkl_vector_minus_vector(self
, &tmp
);