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