hkl/hkl-sample.c

   1 /* This file is part of the hkl library.
   2  *
   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.
   7  *
   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.
  12  *
  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/>.
  15  *
  16  * Copyright (C) 2003-2014 Synchrotron SOLEIL
  17  *                         L'Orme des Merisiers Saint-Aubin
  18  *                         BP 48 91192 GIF-sur-YVETTE CEDEX
  19  *
  20  * Authors: Picca Frédéric-Emmanuel <picca@synchrotron-soleil.fr>
  21  */
  22
  23 /* for strdup */
  24 #define _XOPEN_SOURCE 500
  25 #include <gsl/gsl_errno.h>              // for gsl_set_error_handler, etc
  26 #include <gsl/gsl_multimin.h>           // for gsl_multimin_function, etc
  27 #include <gsl/gsl_nan.h>                // for GSL_NAN
  28 #include <gsl/gsl_vector_double.h>      // for gsl_vector_get, etc
  29 #include <math.h>                       // for M_PI, fabs
  30 #include <stddef.h>                     // for size_t
  31 #include <stdio.h>                      // for fprintf, FILE
  32 #include <stdlib.h>                     // for free
  33 #include <string.h>                     // for NULL, strdup
  34 #include "hkl-detector-private.h"       // for hkl_detector_new_copy, etc
  35 #include "hkl-geometry-private.h"       // for _HklGeometry, etc
  36 #include "hkl-lattice-private.h"        // for _HklLattice, etc
  37 #include "hkl-macros-private.h"         // for HKL_MALLOC
  38 #include "hkl-matrix-private.h"         // for hkl_matrix_init_from_euler, etc
  39 #include "hkl-parameter-private.h"      // for hkl_parameter_fprintf, etc
  40 #include "hkl-quaternion-private.h"     // for hkl_quaternion_conjugate, etc
  41 #include "hkl-sample-private.h"         // for _HklSample, etc
  42 #include "hkl-source-private.h"         // for hkl_source_compute_ki
  43 #include "hkl-unit-private.h"           // for hkl_unit_angle_deg, etc
  44 #include "hkl-vector-private.h"         // for HklVector, hkl_vector_angle, etc
  45 #include "hkl.h"                        // for HklSample, etc
  46 #include "hkl/ccan/darray/darray.h"     // for darray_foreach, darray_item
  47 #include "hkl/ccan/list/list.h"         // for list_head, list_add_tail, etc
  48
  49 /* private */
  50
  51 static void hkl_sample_reflection_update(HklSampleReflection *self)
  52 {
  53         HklVector ki;
  54         HklQuaternion q;
  55
  56         if(!self)
  57                 return;
  58
  59         /* compute the _hkl using only the axes of the geometry */
  60         /* first Q from angles */
  61         hkl_source_compute_ki(&self->geometry->source, &ki);
  62         self->_hkl = ki;
  63         hkl_vector_rotated_quaternion(&self->_hkl,
  64                                       &darray_item(self->geometry->holders, self->detector->idx)->q);
  65         hkl_vector_minus_vector(&self->_hkl, &ki);
  66
  67         q = darray_item(self->geometry->holders, 0)->q;
  68         hkl_quaternion_conjugate(&q);
  69         hkl_vector_rotated_quaternion(&self->_hkl, &q);
  70 }
  71
  72 static void hkl_sample_compute_UxUyUz(HklSample *self)
  73 {
  74         double ux;
  75         double uy;
  76         double uz;
  77
  78         hkl_matrix_to_euler(&self->U, &ux, &uy, &uz);
  79         hkl_parameter_value_set(self->ux, ux, HKL_UNIT_DEFAULT, NULL);
  80         hkl_parameter_value_set(self->uy, uy, HKL_UNIT_DEFAULT, NULL);
  81         hkl_parameter_value_set(self->uz, uz, HKL_UNIT_DEFAULT, NULL);
  82 }
  83
  84 static int hkl_sample_compute_UB(HklSample *self)
  85 {
  86         HklMatrix B;
  87
  88         if (!hkl_lattice_get_B(self->lattice, &B))
  89                 return FALSE;
  90
  91         self->UB = self->U;
  92         hkl_matrix_times_matrix(&self->UB, &B);
  93
  94         return TRUE;
  95 }
  96
  97 /*
  98  * this structure is used by the minimization gsl algorithm.
  99  * in the set_UB method
 100  */
 101 struct set_UB_t
 102 {
 103         HklSample *sample;
 104         const HklMatrix *UB;
 105 };
 106
 107 static double set_UB_fitness(const gsl_vector *x, void *params)
 108 {
 109         size_t i, j;
 110         double fitness;
 111         double euler_x;
 112         double euler_y;
 113         double euler_z;
 114         struct set_UB_t *parameters = params;
 115         HklSample *sample = parameters->sample;
 116         const HklMatrix *UB = parameters->UB;
 117
 118         euler_x = gsl_vector_get(x, 0);
 119         euler_y = gsl_vector_get(x, 1);
 120         euler_z = gsl_vector_get(x, 2);
 121
 122         hkl_parameter_value_set(sample->ux, euler_x, HKL_UNIT_DEFAULT, NULL);
 123         hkl_parameter_value_set(sample->uy, euler_y, HKL_UNIT_DEFAULT, NULL);
 124         hkl_parameter_value_set(sample->uz, euler_z, HKL_UNIT_DEFAULT, NULL);
 125         hkl_parameter_value_set(sample->lattice->a, gsl_vector_get(x, 3), HKL_UNIT_DEFAULT, NULL);
 126         hkl_parameter_value_set(sample->lattice->b, gsl_vector_get(x, 4), HKL_UNIT_DEFAULT, NULL);
 127         hkl_parameter_value_set(sample->lattice->c, gsl_vector_get(x, 5), HKL_UNIT_DEFAULT, NULL);
 128         hkl_parameter_value_set(sample->lattice->alpha, gsl_vector_get(x, 6), HKL_UNIT_DEFAULT, NULL);
 129         hkl_parameter_value_set(sample->lattice->beta, gsl_vector_get(x, 7), HKL_UNIT_DEFAULT, NULL);
 130         hkl_parameter_value_set(sample->lattice->gamma, gsl_vector_get(x, 8), HKL_UNIT_DEFAULT, NULL);
 131
 132         hkl_matrix_init_from_euler(&sample->U, euler_x, euler_y, euler_z);
 133         if (!hkl_sample_compute_UB(sample))
 134                 return GSL_NAN;
 135
 136         fitness = 0.;
 137         for(i=0; i<3; ++i)
 138                 for(j=0; j<3; ++j){
 139                         double tmp = UB->data[i][j] - sample->UB.data[i][j];
 140                         fitness += tmp * tmp;
 141                 }
 142         return fitness;
 143 }
 144
 145 static double mono_crystal_fitness(const gsl_vector *x, void *params)
 146 {
 147         size_t i, j;
 148         double fitness;
 149         double euler_x;
 150         double euler_y;
 151         double euler_z;
 152         HklSample *sample = params;
 153         HklSampleReflection *reflection;
 154
 155         euler_x = gsl_vector_get(x, 0);
 156         euler_y = gsl_vector_get(x, 1);
 157         euler_z = gsl_vector_get(x, 2);
 158
 159         hkl_parameter_value_set(sample->ux, euler_x, HKL_UNIT_DEFAULT, NULL);
 160         hkl_parameter_value_set(sample->uy, euler_y, HKL_UNIT_DEFAULT, NULL);
 161         hkl_parameter_value_set(sample->uz, euler_z, HKL_UNIT_DEFAULT, NULL);
 162         hkl_parameter_value_set(sample->lattice->a, gsl_vector_get(x, 3), HKL_UNIT_DEFAULT, NULL);
 163         hkl_parameter_value_set(sample->lattice->b, gsl_vector_get(x, 4), HKL_UNIT_DEFAULT, NULL);
 164         hkl_parameter_value_set(sample->lattice->c, gsl_vector_get(x, 5), HKL_UNIT_DEFAULT, NULL);
 165         hkl_parameter_value_set(sample->lattice->alpha, gsl_vector_get(x, 6), HKL_UNIT_DEFAULT, NULL);
 166         hkl_parameter_value_set(sample->lattice->beta, gsl_vector_get(x, 7), HKL_UNIT_DEFAULT, NULL);
 167         hkl_parameter_value_set(sample->lattice->gamma, gsl_vector_get(x, 8), HKL_UNIT_DEFAULT, NULL);
 168         hkl_matrix_init_from_euler(&sample->U, euler_x, euler_y, euler_z);
 169         if (!hkl_sample_compute_UB(sample))
 170                 return GSL_NAN;
 171
 172         fitness = 0.;
 173         list_for_each(&sample->reflections, reflection, list){
 174                 if(reflection->flag){
 175                         HklVector UBh;
 176
 177                         UBh = reflection->hkl;
 178                         hkl_matrix_times_vector(&sample->UB, &UBh);
 179
 180                         for(j=0; j<3; ++j) {
 181                                 double tmp = UBh.data[j] - reflection->_hkl.data[j];
 182                                 fitness += tmp * tmp;
 183                         }
 184                 }
 185         }
 186         return fitness;
 187 }
 188
 189 static int minimize(HklSample *sample,
 190                     double (* f) (const gsl_vector * x, void * params),
 191                     void *params, GError **error)
 192 {
 193         gsl_multimin_fminimizer_type const *T = gsl_multimin_fminimizer_nmsimplex;
 194         gsl_multimin_fminimizer *s = NULL;
 195         gsl_vector *ss, *x;
 196         gsl_multimin_function minex_func;
 197         size_t iter = 0;
 198         int status;
 199         double size = 0;
 200
 201         hkl_error (error == NULL || *error == NULL);
 202
 203         /* Starting point */
 204         x = gsl_vector_alloc (9);
 205         gsl_vector_set (x, 0, hkl_parameter_value_get(sample->ux, HKL_UNIT_DEFAULT));
 206         gsl_vector_set (x, 1, hkl_parameter_value_get(sample->uy, HKL_UNIT_DEFAULT));
 207         gsl_vector_set (x, 2, hkl_parameter_value_get(sample->uz, HKL_UNIT_DEFAULT));
 208         gsl_vector_set (x, 3, hkl_parameter_value_get(sample->lattice->a, HKL_UNIT_DEFAULT));
 209         gsl_vector_set (x, 4, hkl_parameter_value_get(sample->lattice->b, HKL_UNIT_DEFAULT));
 210         gsl_vector_set (x, 5, hkl_parameter_value_get(sample->lattice->c, HKL_UNIT_DEFAULT));
 211         gsl_vector_set (x, 6, hkl_parameter_value_get(sample->lattice->alpha, HKL_UNIT_DEFAULT));
 212         gsl_vector_set (x, 7, hkl_parameter_value_get(sample->lattice->beta, HKL_UNIT_DEFAULT));
 213         gsl_vector_set (x, 8, hkl_parameter_value_get(sample->lattice->gamma, HKL_UNIT_DEFAULT));
 214
 215         /* Set initial step sizes to 1 */
 216         ss = gsl_vector_alloc (9);
 217         gsl_vector_set (ss, 0, sample->ux->fit);
 218         gsl_vector_set (ss, 1, sample->uy->fit);
 219         gsl_vector_set (ss, 2, sample->uz->fit);
 220         gsl_vector_set (ss, 3, sample->lattice->a->fit);
 221         gsl_vector_set (ss, 4, sample->lattice->b->fit);
 222         gsl_vector_set (ss, 5, sample->lattice->c->fit);
 223         gsl_vector_set (ss, 6, sample->lattice->alpha->fit);
 224         gsl_vector_set (ss, 7, sample->lattice->beta->fit);
 225         gsl_vector_set (ss, 8, sample->lattice->gamma->fit);
 226
 227         /* Initialize method and iterate */
 228         minex_func.n = 9;
 229         minex_func.f = f;
 230         minex_func.params = params;
 231         s = gsl_multimin_fminimizer_alloc (T, 9);
 232         gsl_set_error_handler_off();
 233         gsl_multimin_fminimizer_set (s, &minex_func, x, ss);
 234         do {
 235                 ++iter;
 236                 status = gsl_multimin_fminimizer_iterate(s);
 237                 if (status)
 238                         break;
 239                 size = gsl_multimin_fminimizer_size (s);
 240                 status = gsl_multimin_test_size (size, HKL_EPSILON / 2.);
 241         } while (status == GSL_CONTINUE && iter < 10000);
 242         gsl_vector_free(x);
 243         gsl_vector_free(ss);
 244         gsl_multimin_fminimizer_free(s);
 245         gsl_set_error_handler (NULL);
 246
 247         if (status == GSL_CONTINUE)
 248                 return FALSE;
 249
 250         return TRUE;
 251 }
 252
 253 /*************/
 254 /* HklSample */
 255 /*************/
 256
 257 /**
 258  * hkl_sample_new:
 259  * @name:
 260  *
 261  * constructor
 262  *
 263  * Returns:
 264  **/
 265 HklSample* hkl_sample_new(const char *name)
 266 {
 267         HklSample *self = NULL;
 268
 269         /* check parameters */
 270         if(!name)
 271                 return self;
 272
 273         self = HKL_MALLOC(HklSample);
 274
 275         self->name = strdup(name);
 276         self->lattice = hkl_lattice_new_default();
 277         hkl_matrix_init(&self->U,1, 0, 0, 0, 1, 0, 0, 0, 1);
 278         hkl_matrix_init(&self->UB,1, 0, 0, 0, 1, 0, 0, 0, 1);
 279
 280         self->ux = hkl_parameter_new("ux", -M_PI, 0., M_PI,
 281                                      TRUE, TRUE,
 282                                      &hkl_unit_angle_rad,
 283                                      &hkl_unit_angle_deg);
 284         self->uy = hkl_parameter_new("uy", -M_PI, 0., M_PI,
 285                                      TRUE, TRUE,
 286                                      &hkl_unit_angle_rad,
 287                                      &hkl_unit_angle_deg);
 288         self->uz = hkl_parameter_new("uz", -M_PI, 0., M_PI,
 289                                      TRUE, TRUE,
 290                                      &hkl_unit_angle_rad,
 291                                      &hkl_unit_angle_deg);
 292
 293         hkl_sample_compute_UB(self);
 294         list_head_init(&self->reflections);
 295         self->n_reflections = 0;
 296
 297         return self;
 298 }
 299
 300 /**
 301  * hkl_sample_new_copy: (skip)
 302  * @self:
 303  *
 304  * copy constructor
 305  *
 306  * Returns:
 307  **/
 308 HklSample *hkl_sample_new_copy(const HklSample *self)
 309 {
 310         HklSample *dup = NULL;
 311         HklSampleReflection *reflection;
 312
 313         /* check parameters */
 314         if(!self)
 315                 return dup;
 316
 317         dup = HKL_MALLOC(HklSample);
 318
 319         dup->name = strdup(self->name);
 320         dup->lattice = hkl_lattice_new_copy(self->lattice);
 321         dup->U = self->U;
 322         dup->UB = self->UB;
 323         dup->ux = hkl_parameter_new_copy(self->ux);
 324         dup->uy = hkl_parameter_new_copy(self->uy);
 325         dup->uz = hkl_parameter_new_copy(self->uz);
 326
 327         /* copy the reflections */
 328         list_head_init(&dup->reflections);
 329         list_for_each(&self->reflections, reflection, list){
 330                 list_add_tail(&dup->reflections,
 331                               &hkl_sample_reflection_new_copy(reflection)->list);
 332         }
 333         dup->n_reflections = self->n_reflections;
 334
 335         return dup;
 336 }
 337
 338 /**
 339  * hkl_sample_free: (skip)
 340  * @self:
 341  *
 342  * destructor
 343  **/
 344 void hkl_sample_free(HklSample *self)
 345 {
 346         HklSampleReflection *reflection;
 347         HklSampleReflection *next;
 348
 349         if (!self)
 350                 return;
 351
 352         free(self->name);
 353         hkl_lattice_free(self->lattice);
 354         hkl_parameter_free(self->ux);
 355         hkl_parameter_free(self->uy);
 356         hkl_parameter_free(self->uz);
 357         list_for_each_safe(&self->reflections, reflection, next, list){
 358                 list_del(&reflection->list);
 359                 hkl_sample_reflection_free(reflection);
 360         }
 361         free(self);
 362 }
 363
 364 /**
 365  * hkl_sample_name_get:
 366  * @self: the this ptr
 367  *
 368  * Return value: the name of the sample
 369  **/
 370 const char *hkl_sample_name_get(const HklSample *self)
 371 {
 372         return self->name;
 373 }
 374
 375 /**
 376  * hkl_sample_name_set:
 377  * @self: the this ptr
 378  * @name: the new name to set
 379  *
 380  * set the name of the sample
 381  **/
 382 void hkl_sample_name_set(HklSample *self, const char *name)
 383 {
 384         if(self->name)
 385                 free(self->name);
 386         self->name = strdup(name);
 387 }
 388
 389 /**
 390  * hkl_sample_lattice_get:
 391  * @self: the this ptr
 392  *
 393  * Return value: the lattice parameters of the sample.
 394  **/
 395 const HklLattice *hkl_sample_lattice_get(HklSample *self)
 396 {
 397         return self->lattice;
 398 }
 399
 400 /**
 401  * hkl_sample_lattice_set:
 402  * @self: the this ptr
 403  * @lattice: the lattice to set
 404  **/
 405 void hkl_sample_lattice_set(HklSample *self, const HklLattice *lattice)
 406 {
 407         hkl_lattice_lattice_set(self->lattice, lattice);
 408         hkl_sample_compute_UB(self);
 409 }
 410
 411 /**
 412  * hkl_sample_ux_get:
 413  * @self: the this ptr
 414  *
 415  * Return value: the ux part of the U matrix.
 416  **/
 417 const HklParameter *hkl_sample_ux_get(const HklSample *self)
 418 {
 419         return self->ux;
 420 }
 421
 422 /**
 423  * hkl_sample_uy_get:
 424  * @self: the this ptr
 425  *
 426  * Return value: the uy part of the U matrix.
 427  **/
 428 const HklParameter *hkl_sample_uy_get(const HklSample *self)
 429 {
 430         return self->uy;
 431 }
 432
 433 /**
 434  * hkl_sample_uz_get:
 435  * @self: the this ptr
 436  *
 437  * Return value: the uz part of the U matrix.
 438  **/
 439 const HklParameter *hkl_sample_uz_get(const HklSample *self)
 440 {
 441         return self->uz;
 442 }
 443
 444 /**
 445  * hkl_sample_ux_set:
 446  * @self: the this ptr
 447  * @ux: the ux parameter to set
 448  * @error: return location for a GError, or NULL
 449  *
 450  * set the ux part of the U matrix.
 451  *
 452  * Returns: TRUE on success, FALSE if an error occurred
 453  **/
 454 int hkl_sample_ux_set(HklSample *self,
 455                       const HklParameter *ux,
 456                       GError **error)
 457 {
 458         hkl_error (error == NULL || *error == NULL);
 459
 460         if(!hkl_parameter_init_copy(self->ux, ux, error)){
 461                 g_assert (error == NULL || *error != NULL);
 462                 return FALSE;
 463         }
 464         g_assert (error == NULL || *error == NULL);
 465
 466         hkl_matrix_init_from_euler(&self->U,
 467                                    hkl_parameter_value_get(self->ux, HKL_UNIT_DEFAULT),
 468                                    hkl_parameter_value_get(self->uy, HKL_UNIT_DEFAULT),
 469                                    hkl_parameter_value_get(self->uz, HKL_UNIT_DEFAULT));
 470         hkl_sample_compute_UB(self);
 471
 472         return TRUE;
 473 }
 474
 475 /**
 476  * hkl_sample_uy_set:
 477  * @self: the this ptr
 478  * @uy: the uy parameter to set
 479  * @error: return location for a GError, or NULL
 480  *
 481  * set the uy part of the U matrix.
 482  *
 483  * Returns: TRUE on success, FALSE if an error occurred
 484  **/
 485 int hkl_sample_uy_set(HklSample *self, const HklParameter *uy,
 486                       GError **error)
 487 {
 488         hkl_error (error == NULL || *error == NULL);
 489
 490         if(!hkl_parameter_init_copy(self->uy, uy, error)){
 491                 g_assert (error == NULL || *error != NULL);
 492                 return FALSE;
 493         }
 494         g_assert (error == NULL || *error == NULL);
 495
 496         hkl_matrix_init_from_euler(&self->U,
 497                                    hkl_parameter_value_get(self->ux, HKL_UNIT_DEFAULT),
 498                                    hkl_parameter_value_get(self->uy, HKL_UNIT_DEFAULT),
 499                                    hkl_parameter_value_get(self->uz, HKL_UNIT_DEFAULT));
 500         hkl_sample_compute_UB(self);
 501 }
 502
 503 /**
 504  * hkl_sample_uz_set:
 505  * @self: the this ptr
 506  * @uz: the uz parameter to set
 507  * @error: return location for a GError, or NULL
 508  *
 509  * set the uz part of the U matrix.
 510  *
 511  * Returns: TRUE on success, FALSE if an error occurred
 512  **/
 513 int hkl_sample_uz_set(HklSample *self, const HklParameter *uz,
 514                       GError **error)
 515 {
 516         hkl_error (error == NULL || *error == NULL);
 517
 518         if(!hkl_parameter_init_copy(self->uz, uz, error)){
 519                 g_assert (error == NULL || *error != NULL);
 520                 return FALSE;
 521         }
 522         g_assert (error == NULL || *error == NULL);
 523
 524         hkl_matrix_init_from_euler(&self->U,
 525                                    hkl_parameter_value_get(self->ux, HKL_UNIT_DEFAULT),
 526                                    hkl_parameter_value_get(self->uy, HKL_UNIT_DEFAULT),
 527                                    hkl_parameter_value_get(self->uz, HKL_UNIT_DEFAULT));
 528         hkl_sample_compute_UB(self);
 529 }
 530
 531 /**
 532  * hkl_sample_U_get:
 533  * @self: the this ptr
 534  *
 535  * Return value: the U matrix of the sample
 536  **/
 537 const HklMatrix *hkl_sample_U_get(const HklSample *self)
 538 {
 539         return &self->U;
 540 }
 541
 542 void hkl_sample_U_set(HklSample *self, const HklMatrix *U)
 543 {
 544         double x, y, z;
 545
 546         hkl_matrix_matrix_set(&self->U, U);
 547         hkl_sample_compute_UB(self);
 548         hkl_matrix_to_euler(&self->U, &x, &y, &z);
 549         hkl_parameter_value_set(self->ux, x, HKL_UNIT_DEFAULT, NULL);
 550         hkl_parameter_value_set(self->uy, y, HKL_UNIT_DEFAULT, NULL);
 551         hkl_parameter_value_set(self->uz, z, HKL_UNIT_DEFAULT, NULL);
 552 }
 553
 554 /**
 555  * hkl_sample_UB_get:
 556  * @self: the this ptr.
 557  *
 558  * Return value: get the UB matrix of the sample
 559  **/
 560 const HklMatrix *hkl_sample_UB_get(const HklSample *self)
 561 {
 562         return &self->UB;
 563 }
 564
 565 /**
 566  * hkl_sample_UB_set:
 567  * @self: the sample to modify
 568  * @UB: the UB matrix to set
 569  *
 570  * Set the UB matrix using an external UB matrix. In fact you give
 571  * the UB matrix but only the U matrix of the sample is affected by
 572  * this operation. We keep the B matrix constant.
 573  * U * B = UB -> U = UB * B^-1
 574  **/
 575 int hkl_sample_UB_set(HklSample *self, const HklMatrix *UB,
 576                       GError **error)
 577 {
 578         struct set_UB_t params = {
 579                 .sample = self,
 580                 .UB = UB
 581         };
 582
 583         return minimize(self, set_UB_fitness, &params, error);
 584 }
 585
 586 /**
 587  * hkl_sample_n_reflections_get: (skip)
 588  * @self: the this ptr
 589  *
 590  * return the number of reflections of the sample
 591  *
 592  * Returns:
 593  **/
 594 size_t hkl_sample_n_reflections_get(const HklSample *self)
 595 {
 596         return self->n_reflections;
 597 }
 598
 599 /**
 600  * hkl_sample_reflections_first_get: (skip)
 601  * @self: the this ptr
 602  *
 603  * Return value: the first HklSampleReflection of the sample.
 604  **/
 605 HklSampleReflection *hkl_sample_reflections_first_get(HklSample *self)
 606 {
 607         return list_top(&self->reflections, HklSampleReflection, list);
 608 }
 609
 610 /**
 611  * hkl_sample_reflections_next_get: (skip)
 612  * @self: the this ptr
 613  * @reflection: the current reflection
 614  *
 615  * Return value: (allow-none): the next reflection or NULL if no more reflection
 616  **/
 617 HklSampleReflection *hkl_sample_reflections_next_get(HklSample *self,
 618                                                      HklSampleReflection *reflection)
 619 {
 620         return list_next(&self->reflections, reflection, list);
 621 }
 622
 623 /**
 624  * hkl_sample_add_reflection: (skip)
 625  * @self: the this ptr
 626  * @reflection: The reflection to add
 627  *
 628  * add a reflection to the sample, if the reflection is already part
 629  * of the sample reflection list, this method does nothing.
 630  **/
 631 void hkl_sample_add_reflection(HklSample *self,
 632                                HklSampleReflection *reflection)
 633 {
 634         HklSampleReflection *ref;
 635
 636         list_for_each(&self->reflections, ref, list){
 637                 if (ref == reflection)
 638                         return;
 639         }
 640
 641         list_add_tail(&self->reflections, &reflection->list);
 642         self->n_reflections++;
 643 }
 644
 645 /**
 646  * hkl_sample_del_reflection:
 647  * @self: the this ptr
 648  * @reflection: the reflection to remove.
 649  *
 650  * remove an HklSampleRefelction from the reflections list.
 651  **/
 652 void hkl_sample_del_reflection(HklSample *self,
 653                                HklSampleReflection *reflection)
 654 {
 655         list_del(&reflection->list);
 656         hkl_sample_reflection_free(reflection);
 657         self->n_reflections--;
 658 }
 659
 660 /**
 661  * hkl_sample_compute_UB_busing_levy:
 662  * @self: the this ptr
 663  * @r1: the first #HklsampleReflection
 664  * @r2: the second #HklSampleReflection
 665  *
 666  * compute the UB matrix using the Busing and Levy method
 667  * #todo: add ref
 668  *
 669  * Returns: 0 or 1 if it succeed
 670  **/
 671 int hkl_sample_compute_UB_busing_levy(HklSample *self,
 672                                       const HklSampleReflection *r1,
 673                                       const HklSampleReflection *r2,
 674                                       GError **error)
 675
 676 {
 677         if (!hkl_vector_is_colinear(&r1->hkl, &r2->hkl)) {
 678                 HklVector h1c;
 679                 HklVector h2c;
 680                 HklMatrix B;
 681                 HklMatrix Tc;
 682
 683                 /* Compute matrix Tc from r1 and r2. */
 684                 h1c = r1->hkl;
 685                 h2c = r2->hkl;
 686                 hkl_lattice_get_B(self->lattice, &B);
 687                 hkl_matrix_times_vector(&B, &h1c);
 688                 hkl_matrix_times_vector(&B, &h2c);
 689                 hkl_matrix_init_from_two_vector(&Tc, &h1c, &h2c);
 690                 hkl_matrix_transpose(&Tc);
 691
 692                 /* compute U */
 693                 hkl_matrix_init_from_two_vector(&self->U,
 694                                                 &r1->_hkl, &r2->_hkl);
 695                 hkl_matrix_times_matrix(&self->U, &Tc);
 696                 hkl_sample_compute_UxUyUz(self);
 697                 hkl_sample_compute_UB(self);
 698         } else
 699                 return FALSE;
 700
 701         return TRUE;
 702 }
 703
 704 /**
 705  * hkl_sample_affine:
 706  * @self: the this ptr
 707  *
 708  * affine the sample
 709  *
 710  * Returns: the fitness of the affined #HklSample
 711  **/
 712 int hkl_sample_affine(HklSample *self, GError **error)
 713 {
 714         return minimize(self, mono_crystal_fitness, self, error);
 715 }
 716
 717 /**
 718  * hkl_sample_get_reflection_mesured_angle:
 719  * @self: the this ptr
 720  * @r1: the first #HklSampleReflection
 721  * @r2: the second #HklSampleReflection
 722  *
 723  * get the mesured angles between two #HklSampleReflection
 724  *
 725  * Returns: the mesured angle beetween them
 726  **/
 727 double hkl_sample_get_reflection_mesured_angle(const HklSample *self,
 728                                                const HklSampleReflection *r1,
 729                                                const HklSampleReflection *r2)
 730 {
 731         return hkl_vector_angle(&r1->_hkl,
 732                                 &r2->_hkl);
 733 }
 734
 735 /**
 736  * hkl_sample_get_reflection_theoretical_angle:
 737  * @self: the this ptr
 738  * @r1: the first #HklSampleReflection
 739  * @r2: the second #HklSampleReflection
 740  *
 741  * get the theoretical angles between two #HklSampleReflection
 742  *
 743  * Returns: the theoretical angle beetween them
 744  **/
 745 double hkl_sample_get_reflection_theoretical_angle(const HklSample *self,
 746                                                    const HklSampleReflection *r1,
 747                                                    const HklSampleReflection *r2)
 748 {
 749         HklVector hkl1;
 750         HklVector hkl2;
 751
 752         hkl1 = r1->hkl;
 753         hkl2 = r2->hkl;
 754         hkl_matrix_times_vector(&self->UB, &hkl1);
 755         hkl_matrix_times_vector(&self->UB, &hkl2);
 756
 757         return hkl_vector_angle(&hkl1, &hkl2);
 758 }
 759
 760 /**
 761  * hkl_sample_fprintf: (skip)
 762  * @f:
 763  * @self:
 764  *
 765  * print to a file a sample
 766  **/
 767 void hkl_sample_fprintf(FILE *f, const HklSample *self)
 768 {
 769         if(!self)
 770                 return;
 771
 772         fprintf(f, "\nSample name: \"%s\"", self->name);
 773
 774         fprintf(f, "\nLattice parameters:");
 775         fprintf(f, "\n ");
 776         hkl_parameter_fprintf(f, self->lattice->a);
 777         fprintf(f, "\n ");
 778         hkl_parameter_fprintf(f, self->lattice->b);
 779         fprintf(f, "\n ");
 780         hkl_parameter_fprintf(f, self->lattice->c);
 781         fprintf(f, "\n ");
 782         hkl_parameter_fprintf(f, self->lattice->alpha);
 783         fprintf(f, "\n ");
 784         hkl_parameter_fprintf(f, self->lattice->beta);
 785         fprintf(f, "\n ");
 786         hkl_parameter_fprintf(f, self->lattice->gamma);
 787         fprintf(f, "\n ");
 788         hkl_parameter_fprintf(f, self->ux);
 789         fprintf(f, "\n ");
 790         hkl_parameter_fprintf(f, self->uy);
 791         fprintf(f, "\n ");
 792         hkl_parameter_fprintf(f, self->uz);
 793         fprintf(f, "\nUB:\n");
 794         hkl_matrix_fprintf(f, &self->UB);
 795
 796         if (!list_empty(&self->reflections)){
 797                 HklSampleReflection *reflection;
 798                 HklParameter **axis;
 799
 800                 reflection  = list_top(&self->reflections, HklSampleReflection, list);
 801
 802                 fprintf(f, "Reflections:");
 803                 fprintf(f, "\n");
 804                 fprintf(f, "%-10.6s %-10.6s %-10.6s", "h", "k", "l");
 805                 darray_foreach(axis, reflection->geometry->axes){
 806                         fprintf(f, " %-10.6s", (*axis)->name);
 807                 }
 808
 809                 list_for_each(&self->reflections, reflection, list){
 810                         fprintf(f, "\n%-10.6f %-10.6f %-10.6f",
 811                                 reflection->hkl.data[0],
 812                                 reflection->hkl.data[1],
 813                                 reflection->hkl.data[2]);
 814                         darray_foreach(axis, reflection->geometry->axes){
 815                                 fprintf(f, " %-10.6f", hkl_parameter_value_get(*axis, HKL_UNIT_USER));
 816                         }
 817                 }
 818         }
 819 }
 820
 821 /***********************/
 822 /* HklSampleReflection */
 823 /***********************/
 824
 825 /**
 826  * hkl_sample_reflection_new:
 827  * @geometry:
 828  * @detector:
 829  * @h:
 830  * @k:
 831  * @l:
 832  *
 833  * constructeur
 834  *
 835  * Returns:
 836  **/
 837 HklSampleReflection *hkl_sample_reflection_new(const HklGeometry *geometry,
 838                                                const HklDetector *detector,
 839                                                double h, double k, double l,
 840                                                GError **error)
 841 {
 842         HklSampleReflection *self;
 843
 844         if (!geometry || !detector)
 845                 return NULL;
 846
 847         self = HKL_MALLOC(HklSampleReflection);
 848
 849         self->geometry = hkl_geometry_new_copy(geometry);
 850         self->detector = hkl_detector_new_copy(detector);
 851         self->hkl.data[0] = h;
 852         self->hkl.data[1] = k;
 853         self->hkl.data[2] = l;
 854         self->flag = TRUE;
 855
 856         hkl_sample_reflection_update(self);
 857
 858         return self;
 859 }
 860
 861 /**
 862  * hkl_sample_reflection_new_copy: (skip)
 863  * @self:
 864  *
 865  * copy constructor
 866  *
 867  * Returns:
 868  **/
 869 HklSampleReflection *hkl_sample_reflection_new_copy(const HklSampleReflection *self)
 870 {
 871         HklSampleReflection *dup = NULL;
 872
 873         dup = HKL_MALLOC(HklSampleReflection);
 874
 875         dup->geometry = hkl_geometry_new_copy(self->geometry);
 876         dup->detector = hkl_detector_new_copy(self->detector);
 877         dup->hkl = self->hkl;
 878         dup->_hkl = self->_hkl;
 879         dup->flag = self->flag;
 880
 881         return dup;
 882 }
 883
 884 /**
 885  * hkl_sample_reflection_free: (skip)
 886  * @self:
 887  *
 888  * destructor
 889  **/
 890 void hkl_sample_reflection_free(HklSampleReflection *self)
 891 {
 892         hkl_geometry_free(self->geometry);
 893         hkl_detector_free(self->detector);
 894         free(self);
 895 }
 896
 897 /**
 898  * hkl_sample_reflection_hkl_get:
 899  * @self: the this ptr
 900  * @h: (out caller-allocates): the h-coordinate of the #HklSampleReflection
 901  * @k: (out caller-allocates): the k-coordinate of the #HklSampleReflection
 902  * @l: (out caller-allocates): the l-coordinate of the #HklSampleReflection
 903  *
 904  * get the hkl coordinates of the #HklSampleReflection
 905  **/
 906 void hkl_sample_reflection_hkl_get(const HklSampleReflection *self,
 907                                    double *h, double *k, double *l)
 908 {
 909         *h = self->hkl.data[0];
 910         *k = self->hkl.data[1];
 911         *l = self->hkl.data[2];
 912 }
 913
 914 /**
 915  * hkl_sample_reflection_hkl_set:
 916  * @self: the this ptr
 917  * @h: the h-coordinate of the #HklSampleReflection
 918  * @k: the k-coordinate of the #HklSampleReflection
 919  * @l: the l-coordinate of the #HklSampleReflection
 920  * @error: return location for a GError, or NULL
 921  *
 922  * set the hkl coordinates of the #HklSampleReflection
 923  *
 924  * Returns: TRUE on success, FALSE if an error occurred
 925  **/
 926 int hkl_sample_reflection_hkl_set(HklSampleReflection *self,
 927                                   double h, double k, double l,
 928                                   GError **error)
 929 {
 930         hkl_error (error == NULL || *error == NULL);
 931
 932         if((fabs(h) + fabs(k) + fabs(l) < HKL_EPSILON)){
 933                 g_set_error(error,
 934                             HKL_SAMPLE_REFLECTION_ERROR,
 935                             HKL_SAMPLE_REFLECTION_ERROR_HKL_SET,
 936                             "it is not allow to set a null hkl reflection\n");
 937                 return FALSE;
 938         }
 939
 940         self->hkl.data[0] = h;
 941         self->hkl.data[1] = k;
 942         self->hkl.data[2] = l;
 943
 944         return TRUE;
 945 }
 946
 947 /**
 948  * hkl_sample_reflection_flag_get: (skip)
 949  * @self:
 950  *
 951  * get the flag of the reflection
 952  **/
 953 int hkl_sample_reflection_flag_get(const HklSampleReflection *self)
 954 {
 955         return self->flag;
 956 }
 957
 958 /**
 959  * hkl_sample_reflection_flag_set: (skip)
 960  * @self:
 961  * @flag:
 962  *
 963  * set the flag of the reglection
 964  **/
 965 void hkl_sample_reflection_flag_set(HklSampleReflection *self, int flag)
 966 {
 967         self->flag = flag;
 968 }
 969
 970 /**
 971  * hkl_sample_reflection_geometry_get: (skip)
 972  * @self:
 973  *
 974  * set the geometry of the reflection
 975  **/
 976 const HklGeometry *hkl_sample_reflection_geometry_get(HklSampleReflection *self)
 977 {
 978         return self->geometry;
 979 }
 980
 981 /**
 982  * hkl_sample_reflection_geometry_set: (skip)
 983  * @self:
 984  * @geometry:
 985  *
 986  * set the geometry of the reflection
 987  **/
 988 void hkl_sample_reflection_geometry_set(HklSampleReflection *self,
 989                                         const HklGeometry *geometry)
 990 {
 991         if(self->geometry){
 992                 if(self->geometry != geometry){
 993                         hkl_geometry_free(self->geometry);
 994                         self->geometry = hkl_geometry_new_copy(geometry);
 995                 }
 996         }else
 997                 self->geometry = hkl_geometry_new_copy(geometry);
 998
 999         hkl_sample_reflection_update(self);
1000 }