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>
24 #define _XOPEN_SOURCE 500
27 #include <gsl/gsl_multimin.h>
29 #include <hkl/hkl-sample.h>
30 #include <hkl/hkl-matrix.h>
34 static void hkl_sample_reflection_update(HklSampleReflection
*self
)
42 /* compute the _hkl using only the axes of the geometry */
43 /* first Q from angles */
44 hkl_source_compute_ki(&self
->geometry
->source
, &ki
);
46 hkl_vector_rotated_quaternion(&self
->_hkl
, &self
->geometry
->holders
[self
->detector
.idx
].q
);
47 hkl_vector_minus_vector(&self
->_hkl
, &ki
);
49 q
= self
->geometry
->holders
[0].q
;
50 hkl_quaternion_conjugate(&q
);
51 hkl_vector_rotated_quaternion(&self
->_hkl
, &q
);
54 static HklSampleReflection
*hkl_sample_reflection_new(HklGeometry
*geometry
,
55 HklDetector
const *detector
,
56 double h
, double k
, double l
)
58 HklSampleReflection
*self
;
60 if (!geometry
|| !detector
)
63 self
= HKL_MALLOC(HklSampleReflection
);
65 hkl_geometry_update(geometry
);
67 self
->geometry
= hkl_geometry_new_copy(geometry
);
68 self
->detector
= *detector
;
69 self
->hkl
.data
[0] = h
;
70 self
->hkl
.data
[1] = k
;
71 self
->hkl
.data
[2] = l
;
72 self
->flag
= HKL_TRUE
;
74 hkl_sample_reflection_update(self
);
79 static HklSampleReflection
*hkl_sample_reflection_new_copy(HklSampleReflection
const *src
)
81 HklSampleReflection
*self
= NULL
;
83 self
= HKL_MALLOC(HklSampleReflection
);
85 self
->geometry
= hkl_geometry_new_copy(src
->geometry
);
86 self
->detector
= src
->detector
;
88 self
->_hkl
= src
->_hkl
;
89 self
->flag
= src
->flag
;
94 static void hkl_sample_reflection_free(HklSampleReflection
*self
)
96 hkl_geometry_free(self
->geometry
);
100 static void hkl_sample_compute_UxUyUz(HklSample
*self
)
106 hkl_matrix_to_euler(&self
->U
, &ux
, &uy
, &uz
);
107 hkl_parameter_set_value(self
->ux
, ux
);
108 hkl_parameter_set_value(self
->uy
, uy
);
109 hkl_parameter_set_value(self
->uz
, uz
);
112 static int hkl_sample_compute_UB(HklSample
*self
)
116 if (hkl_lattice_get_B(self
->lattice
, &B
))
120 hkl_matrix_times_matrix(&self
->UB
, &B
);
126 * this structure is used by the minimization gsl algorithm.
127 * in the set_UB method
135 static double set_UB_fitness(gsl_vector
const *x
, void *params
)
142 struct set_UB_t
*parameters
= params
;
143 HklSample
*sample
= parameters
->sample
;
144 const HklMatrix
*UB
= parameters
->UB
;
146 sample
->ux
->value
= euler_x
= gsl_vector_get(x
, 0);
147 sample
->uy
->value
= euler_y
= gsl_vector_get(x
, 1);
148 sample
->uz
->value
= euler_z
= gsl_vector_get(x
, 2);
149 sample
->lattice
->a
->value
= gsl_vector_get(x
, 3);
150 sample
->lattice
->b
->value
= gsl_vector_get(x
, 4);
151 sample
->lattice
->c
->value
= gsl_vector_get(x
, 5);
152 sample
->lattice
->alpha
->value
= gsl_vector_get(x
, 6);
153 sample
->lattice
->beta
->value
= gsl_vector_get(x
, 7);
154 sample
->lattice
->gamma
->value
= gsl_vector_get(x
, 8);
155 hkl_matrix_init_from_euler(&sample
->U
, euler_x
, euler_y
, euler_z
);
156 if (hkl_sample_compute_UB(sample
))
162 double tmp
= UB
->data
[i
][j
] - sample
->UB
.data
[i
][j
];
163 fitness
+= tmp
* tmp
;
168 static double mono_crystal_fitness(gsl_vector
const *x
, void *params
)
175 HklSample
*sample
= params
;
177 sample
->ux
->value
= euler_x
= gsl_vector_get(x
, 0);
178 sample
->uy
->value
= euler_y
= gsl_vector_get(x
, 1);
179 sample
->uz
->value
= euler_z
= gsl_vector_get(x
, 2);
180 sample
->lattice
->a
->value
= gsl_vector_get(x
, 3);
181 sample
->lattice
->b
->value
= gsl_vector_get(x
, 4);
182 sample
->lattice
->c
->value
= gsl_vector_get(x
, 5);
183 sample
->lattice
->alpha
->value
= gsl_vector_get(x
, 6);
184 sample
->lattice
->beta
->value
= gsl_vector_get(x
, 7);
185 sample
->lattice
->gamma
->value
= gsl_vector_get(x
, 8);
186 hkl_matrix_init_from_euler(&sample
->U
, euler_x
, euler_y
, euler_z
);
187 if (hkl_sample_compute_UB(sample
))
191 for(i
=0; i
<HKL_LIST_LEN(sample
->reflections
); ++i
) {
192 HklSampleReflection
*reflection
;
194 reflection
= sample
->reflections
[i
];
195 if(reflection
->flag
== HKL_TRUE
){
198 UBh
= reflection
->hkl
;
199 hkl_matrix_times_vector(&sample
->UB
, &UBh
);
202 double tmp
= UBh
.data
[j
] - reflection
->_hkl
.data
[j
];
203 fitness
+= tmp
* tmp
;
210 static double minimize(HklSample
*sample
, double (* f
) (const gsl_vector
* x
, void * params
), void *params
)
212 gsl_multimin_fminimizer_type
const *T
= gsl_multimin_fminimizer_nmsimplex
;
213 gsl_multimin_fminimizer
*s
= NULL
;
215 gsl_multimin_function minex_func
;
224 x
= gsl_vector_alloc (9);
225 gsl_vector_set (x
, 0, sample
->ux
->value
);
226 gsl_vector_set (x
, 1, sample
->uy
->value
);
227 gsl_vector_set (x
, 2, sample
->uz
->value
);
228 gsl_vector_set (x
, 3, sample
->lattice
->a
->value
);
229 gsl_vector_set (x
, 4, sample
->lattice
->b
->value
);
230 gsl_vector_set (x
, 5, sample
->lattice
->c
->value
);
231 gsl_vector_set (x
, 6, sample
->lattice
->alpha
->value
);
232 gsl_vector_set (x
, 7, sample
->lattice
->beta
->value
);
233 gsl_vector_set (x
, 8, sample
->lattice
->gamma
->value
);
235 /* Set initial step sizes to 1 */
236 ss
= gsl_vector_alloc (9);
237 gsl_vector_set (ss
, 0, sample
->ux
->fit
);
238 gsl_vector_set (ss
, 1, sample
->uy
->fit
);
239 gsl_vector_set (ss
, 2, sample
->uz
->fit
);
240 gsl_vector_set (ss
, 3, sample
->lattice
->a
->fit
);
241 gsl_vector_set (ss
, 4, sample
->lattice
->b
->fit
);
242 gsl_vector_set (ss
, 5, sample
->lattice
->c
->fit
);
243 gsl_vector_set (ss
, 6, sample
->lattice
->alpha
->fit
);
244 gsl_vector_set (ss
, 7, sample
->lattice
->beta
->fit
);
245 gsl_vector_set (ss
, 8, sample
->lattice
->gamma
->fit
);
247 /* Initialize method and iterate */
250 minex_func
.params
= params
;
251 s
= gsl_multimin_fminimizer_alloc (T
, 9);
252 gsl_set_error_handler_off();
253 gsl_multimin_fminimizer_set (s
, &minex_func
, x
, ss
);
256 status
= gsl_multimin_fminimizer_iterate(s
);
259 size
= gsl_multimin_fminimizer_size (s
);
260 status
= gsl_multimin_test_size (size
, HKL_EPSILON
/ 2.);
261 } while (status
== GSL_CONTINUE
&& iter
< 10000);
264 gsl_multimin_fminimizer_free(s
);
265 gsl_set_error_handler (NULL
);
274 HklSample
* hkl_sample_new(char const *name
, HklSampleType type
)
276 HklSample
*self
= NULL
;
278 /* check parameters */
282 self
= HKL_MALLOC(HklSample
);
284 self
->name
= strdup(name
);
286 self
->lattice
= hkl_lattice_new_default();
287 hkl_matrix_init(&self
->U
,1, 0, 0, 0, 1, 0, 0, 0, 1);
288 hkl_matrix_init(&self
->UB
,1, 0, 0, 0, 1, 0, 0, 0, 1);
290 self
->ux
= hkl_parameter_new("ux", -M_PI
, 0., M_PI
,
293 &hkl_unit_angle_deg
);
294 self
->uy
= hkl_parameter_new("uy", -M_PI
, 0., M_PI
,
297 &hkl_unit_angle_deg
);
298 self
->uz
= hkl_parameter_new("uz", -M_PI
, 0., M_PI
,
301 &hkl_unit_angle_deg
);
303 hkl_sample_compute_UB(self
);
304 HKL_LIST_INIT(self
->reflections
);
309 HklSample
*hkl_sample_new_copy(HklSample
const *src
)
311 HklSample
*self
= NULL
;
315 /* check parameters */
319 self
= HKL_MALLOC(HklSample
);
321 self
->name
= strdup(src
->name
);
322 self
->type
= src
->type
;
323 self
->lattice
= hkl_lattice_new_copy(src
->lattice
);
326 self
->ux
= hkl_parameter_new_copy(src
->ux
);
327 self
->uy
= hkl_parameter_new_copy(src
->uy
);
328 self
->uz
= hkl_parameter_new_copy(src
->uz
);
330 /* copy the reflections */
331 len
= HKL_LIST_LEN(src
->reflections
);
332 HKL_LIST_ALLOC(self
->reflections
, len
);
334 self
->reflections
[i
] = hkl_sample_reflection_new_copy(src
->reflections
[i
]);
339 void hkl_sample_free(HklSample
*self
)
345 hkl_lattice_free(self
->lattice
);
346 hkl_parameter_free(self
->ux
);
347 hkl_parameter_free(self
->uy
);
348 hkl_parameter_free(self
->uz
);
349 HKL_LIST_FREE_DESTRUCTOR(self
->reflections
, hkl_sample_reflection_free
);
353 void hkl_sample_set_name(HklSample
*self
, char const *name
)
360 self
->name
= strdup(name
);
363 int hkl_sample_set_lattice(HklSample
*self
,
364 double a
, double b
, double c
,
365 double alpha
, double beta
, double gamma
)
373 status
= hkl_lattice_set(self
->lattice
, a
, b
, c
, alpha
, beta
, gamma
);
374 if (status
== HKL_SUCCESS
)
375 hkl_sample_compute_UB(self
);
380 int hkl_sample_set_U_from_euler(HklSample
*self
,
381 double x
, double y
, double z
)
386 hkl_matrix_init_from_euler(&self
->U
, x
, y
, z
);
387 hkl_sample_compute_UB(self
);
388 hkl_parameter_set_value(self
->ux
, x
);
389 hkl_parameter_set_value(self
->uy
, y
);
390 hkl_parameter_set_value(self
->uz
, z
);
395 void hkl_sample_get_UB(HklSample
*self
, HklMatrix
*UB
)
400 hkl_sample_compute_UB(self
);
405 * hkl_sample_set_UB: set the UB matrix of the sample
406 * @self: the sample to modify
407 * @UB: the UB matrix to set
409 * Set the UB matrix using an external UB matrix. In fact you give
410 * the UB matrix but only the U matrix of the sample is affected by
411 * this operation. We keep the B matrix constant.
412 * U * B = UB -> U = UB * B^-1
414 double hkl_sample_set_UB(HklSample
*self
, const HklMatrix
*UB
)
416 struct set_UB_t params
;
421 params
.sample
= self
;
424 return minimize(self
, set_UB_fitness
, ¶ms
);
427 HklSampleReflection
*hkl_sample_add_reflection(HklSample
*self
,
428 HklGeometry
*geometry
,
429 HklDetector
const *detector
,
430 double h
, double k
, double l
)
432 HklSampleReflection
*ref
;
434 ref
= hkl_sample_reflection_new(geometry
, detector
, h
, k
, l
);
437 HKL_LIST_ADD_VALUE(self
->reflections
, ref
);
442 HklSampleReflection
* hkl_sample_get_ith_reflection(HklSample
const *self
, size_t idx
)
447 return self
->reflections
[idx
];
450 int hkl_sample_del_reflection(HklSample
*self
, size_t idx
)
455 hkl_sample_reflection_free(self
->reflections
[idx
]);
456 HKL_LIST_DEL(self
->reflections
, idx
);
461 int hkl_sample_compute_UB_busing_levy(HklSample
*self
, size_t idx1
, size_t idx2
)
463 HklSampleReflection
*r1
;
464 HklSampleReflection
*r2
;
467 || idx1
>= HKL_LIST_LEN(self
->reflections
)
468 || idx2
>= HKL_LIST_LEN(self
->reflections
))
471 r1
= self
->reflections
[idx1
];
472 r2
= self
->reflections
[idx2
];
474 if (!hkl_vector_is_colinear(&r1
->hkl
, &r2
->hkl
)) {
480 /* Compute matrix Tc from r1 and r2. */
483 hkl_lattice_get_B(self
->lattice
, &B
);
484 hkl_matrix_times_vector(&B
, &h1c
);
485 hkl_matrix_times_vector(&B
, &h2c
);
486 hkl_matrix_init_from_two_vector(&Tc
, &h1c
, &h2c
);
487 hkl_matrix_transpose(&Tc
);
490 hkl_matrix_init_from_two_vector(&self
->U
,
491 &r1
->_hkl
, &r2
->_hkl
);
492 hkl_matrix_times_matrix(&self
->U
, &Tc
);
493 hkl_sample_compute_UxUyUz(self
);
494 hkl_sample_compute_UB(self
);
501 double hkl_sample_affine(HklSample
*self
)
506 return minimize(self
, mono_crystal_fitness
, self
);
509 double hkl_sample_get_reflection_mesured_angle(HklSample
const *self
,
510 size_t idx1
, size_t idx2
)
513 || idx1
>= HKL_LIST_LEN(self
->reflections
)
514 || idx2
>= HKL_LIST_LEN(self
->reflections
))
517 return hkl_vector_angle(&self
->reflections
[idx1
]->_hkl
,
518 &self
->reflections
[idx2
]->_hkl
);
521 double hkl_sample_get_reflection_theoretical_angle(HklSample
const *self
,
522 size_t idx1
, size_t idx2
)
528 || idx1
>= HKL_LIST_LEN(self
->reflections
)
529 || idx2
>= HKL_LIST_LEN(self
->reflections
))
532 hkl1
= self
->reflections
[idx1
]->hkl
;
533 hkl2
= self
->reflections
[idx2
]->hkl
;
534 hkl_matrix_times_vector(&self
->UB
, &hkl1
);
535 hkl_matrix_times_vector(&self
->UB
, &hkl2
);
537 return hkl_vector_angle(&hkl1
, &hkl2
);
540 void hkl_sample_fprintf(FILE *f
, HklSample
const *self
)
544 fprintf(f
, "\nSample name: \"%s\"", self
->name
);
546 fprintf(f
, "\nLattice parameters:");
548 hkl_parameter_fprintf(f
, self
->lattice
->a
);
550 hkl_parameter_fprintf(f
, self
->lattice
->b
);
552 hkl_parameter_fprintf(f
, self
->lattice
->c
);
554 hkl_parameter_fprintf(f
, self
->lattice
->alpha
);
556 hkl_parameter_fprintf(f
, self
->lattice
->beta
);
558 hkl_parameter_fprintf(f
, self
->lattice
->gamma
);
560 hkl_parameter_fprintf(f
, self
->ux
);
562 hkl_parameter_fprintf(f
, self
->uy
);
564 hkl_parameter_fprintf(f
, self
->uz
);
565 fprintf(f
, "\nUB:\n");
566 hkl_matrix_fprintf(f
, &self
->UB
);
568 len
= HKL_LIST_LEN(self
->reflections
);
570 HklSampleReflection
*reflection
;
574 reflection
= hkl_sample_get_ith_reflection(self
, 0);
576 fprintf(f
, "Reflections:");
578 fprintf(f
, "i %-10.6s %-10.6s %-10.6s", "h", "k", "l");
579 axes
= reflection
->geometry
->axes
;
580 axes_len
= HKL_LIST_LEN(reflection
->geometry
->axes
);
581 for(i
=0; i
<axes_len
; ++i
)
582 fprintf(f
, " %-10.6s", hkl_axis_get_name(&axes
[i
]));
584 for(i
=0; i
<len
; ++i
){
587 reflection
= hkl_sample_get_ith_reflection(self
, i
);
588 axes
= reflection
->geometry
->axes
;
589 axes_len
= HKL_LIST_LEN(reflection
->geometry
->axes
);
590 fprintf(f
, "\n%d %-10.6f %-10.6f %-10.6f", i
,
591 reflection
->hkl
.data
[0], reflection
->hkl
.data
[1], reflection
->hkl
.data
[2]);
592 for(j
=0; j
<axes_len
; ++j
)
593 fprintf(f
, " %-10.6f", hkl_axis_get_value_unit(&axes
[j
]));
598 /***********************/
599 /* HklSampleReflection */
600 /***********************/
602 void hkl_sample_reflection_set_hkl(HklSampleReflection
*self
, double h
, double k
, double l
)
605 || (fabs(h
) + fabs(k
) + fabs(l
) < HKL_EPSILON
))
608 self
->hkl
.data
[0] = h
;
609 self
->hkl
.data
[1] = k
;
610 self
->hkl
.data
[2] = l
;
613 void hkl_sample_reflection_set_flag(HklSampleReflection
*self
, int flag
)
620 void hkl_sample_reflection_set_geometry(HklSampleReflection
*self
, HklGeometry
*geometry
)
622 if(!self
|| !geometry
)
626 if(self
->geometry
!= geometry
){
627 hkl_geometry_free(self
->geometry
);
628 self
->geometry
= hkl_geometry_new_copy(geometry
);
631 self
->geometry
= hkl_geometry_new_copy(geometry
);
633 hkl_sample_reflection_update(self
);
640 HklSampleList
*hkl_sample_list_new(void)
642 HklSampleList
*self
= NULL
;
643 self
= HKL_MALLOC(HklSampleList
);
645 HKL_LIST_INIT(self
->samples
);
646 self
->current
= NULL
;
651 void hkl_sample_list_free(HklSampleList
*self
)
654 HKL_LIST_FREE_DESTRUCTOR(self
->samples
, hkl_sample_free
);
655 self
->current
= NULL
;
660 void hkl_sample_list_clear(HklSampleList
*self
)
663 HKL_LIST_FREE_DESTRUCTOR(self
->samples
, hkl_sample_free
);
664 self
->current
= NULL
;
665 HKL_LIST_INIT(self
->samples
);
669 HklSample
*hkl_sample_list_append(HklSampleList
*self
, HklSample
*sample
)
672 || hkl_sample_list_get_idx_from_name(self
, sample
->name
) != HKL_FAIL
)
675 HKL_LIST_ADD_VALUE(self
->samples
, sample
);
680 void hkl_sample_list_del(HklSampleList
*self
, HklSample
*sample
)
683 if (self
->current
== sample
)
684 self
->current
= NULL
;
685 HKL_LIST_DEL_ITEM_DESTRUCTOR(self
->samples
, sample
, hkl_sample_free
);
690 size_t hkl_sample_list_len(HklSampleList
const *self
)
692 return HKL_LIST_LEN(self
->samples
);
696 HklSample
*hkl_sample_list_get_ith(HklSampleList
*self
, size_t idx
)
698 HklSample
*sample
= NULL
;
699 if (self
&& idx
< HKL_LIST_LEN(self
->samples
))
700 sample
= self
->samples
[idx
];
706 * hkl_sample_list_get_by_name:
707 * @self: the #HklSampleList
708 * @name: the name of the #HklSample you are looking for.
710 * get the @name named #HklSample from the #HklSampleList.
712 * Returns: an #HklSample or NULL if not present in the #HklSampleList
716 HklSample
*hkl_sample_list_get_by_name(HklSampleList
*self
, char const *name
)
718 HklSample
*sample
= NULL
;
724 idx
= hkl_sample_list_get_idx_from_name(self
, name
);
726 sample
= self
->samples
[idx
];
732 * hkl_sample_list_get_idx_from_name:
733 * @self: the #HklSampleList
734 * @name: the name of the #HklSample.
736 * find the named @name #HklSample in the #HklSampleList and return its index.
738 * Returns: the index or -1 if the #HklSample is not present.
740 * @todo: this method should return an int and not a size_t, as -1 is returned if the
741 * sample doesn not exist in the HklSampleList BUG
744 size_t hkl_sample_list_get_idx_from_name(HklSampleList
*self
, char const *name
)
748 if (!self
|| !name
|| !self
->samples
)
751 for(idx
=0; idx
<HKL_LIST_LEN(self
->samples
); ++idx
)
752 if (!strcmp(self
->samples
[idx
]->name
, name
))
758 int hkl_sample_list_select_current(HklSampleList
*self
, char const *name
)
763 if(!self
|| !name
|| !self
->samples
)
766 idx
= hkl_sample_list_get_idx_from_name(self
, name
);
767 if (idx
!= HKL_FAIL
){
768 self
->current
= self
->samples
[idx
];
775 void hkl_sample_list_fprintf(FILE *f
, HklSampleList
const *self
)
778 for(i
=0; i
<HKL_LIST_LEN(self
->samples
); ++i
)
779 hkl_sample_fprintf(f
, self
->samples
[i
]);