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-2017 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;
 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(i=0; i<3; ++i) {
 233                                 double tmp = UBh.data[i] - reflection->_hkl.data[i];
 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         return TRUE;
 561 }
 562
 563 /**
 564  * hkl_sample_uz_set:
 565  * @self: the this ptr
 566  * @uz: the uz parameter to set
 567  * @error: return location for a GError, or NULL
 568  *
 569  * set the uz part of the U matrix.
 570  *
 571  * Returns: TRUE on success, FALSE if an error occurred
 572  **/
 573 int hkl_sample_uz_set(HklSample *self, const HklParameter *uz,
 574                       GError **error)
 575 {
 576         hkl_error (error == NULL || *error == NULL);
 577
 578         if(!hkl_parameter_init_copy(self->uz, uz, error)){
 579                 g_assert (error == NULL || *error != NULL);
 580                 return FALSE;
 581         }
 582         g_assert (error == NULL || *error == NULL);
 583
 584         hkl_matrix_init_from_euler(&self->U,
 585                                    hkl_parameter_value_get(self->ux, HKL_UNIT_DEFAULT),
 586                                    hkl_parameter_value_get(self->uy, HKL_UNIT_DEFAULT),
 587                                    hkl_parameter_value_get(self->uz, HKL_UNIT_DEFAULT));
 588         hkl_sample_compute_UB(self);
 589
 590         return TRUE;
 591 }
 592
 593 /**
 594  * hkl_sample_U_get:
 595  * @self: the this ptr
 596  *
 597  * Return value: the U matrix of the sample
 598  **/
 599 const HklMatrix *hkl_sample_U_get(const HklSample *self)
 600 {
 601         return &self->U;
 602 }
 603
 604 /*
 605  * TODO implemente the error
 606  */
 607 void hkl_sample_U_set(HklSample *self, const HklMatrix *U, GError **error)
 608 {
 609         double x, y, z;
 610
 611         hkl_matrix_matrix_set(&self->U, U);
 612         hkl_sample_compute_UB(self);
 613         hkl_matrix_to_euler(&self->U, &x, &y, &z);
 614         hkl_parameter_value_set(self->ux, x, HKL_UNIT_DEFAULT, NULL);
 615         hkl_parameter_value_set(self->uy, y, HKL_UNIT_DEFAULT, NULL);
 616         hkl_parameter_value_set(self->uz, z, HKL_UNIT_DEFAULT, NULL);
 617 }
 618
 619 /**
 620  * hkl_sample_UB_get:
 621  * @self: the this ptr.
 622  *
 623  * Return value: get the UB matrix of the sample
 624  **/
 625 const HklMatrix *hkl_sample_UB_get(const HklSample *self)
 626 {
 627         return &self->UB;
 628 }
 629
 630 /**
 631  * hkl_sample_UB_set:
 632  * @self: the sample to modify
 633  * @UB: the UB matrix to set
 634  * @error: error set in case of impossibility
 635  *
 636  * Set the UB matrix using an external UB matrix. In fact you give
 637  * the UB matrix but only the U matrix of the sample is affected by
 638  * this operation. We keep the B matrix constant.
 639  * U * B = UB -> U = UB * B^-1
 640  * TODO implemente the error
 641  **/
 642 int hkl_sample_UB_set(HklSample *self, const HklMatrix *UB,
 643                       GError **error)
 644 {
 645         struct set_UB_t params = {
 646                 .sample = self,
 647                 .UB = *UB
 648         };
 649
 650         return minimize(self, set_UB_fitness, &params, error);
 651 }
 652
 653 /**
 654  * hkl_sample_n_reflections_get: (skip)
 655  * @self: the this ptr
 656  *
 657  * return the number of reflections of the sample
 658  *
 659  * Returns:
 660  **/
 661 size_t hkl_sample_n_reflections_get(const HklSample *self)
 662 {
 663         return self->n_reflections;
 664 }
 665
 666 /**
 667  * hkl_sample_reflections_first_get: (skip)
 668  * @self: the this ptr
 669  *
 670  * Return value: the first HklSampleReflection of the sample.
 671  **/
 672 HklSampleReflection *hkl_sample_reflections_first_get(HklSample *self)
 673 {
 674         return list_top(&self->reflections, HklSampleReflection, list);
 675 }
 676
 677 /**
 678  * hkl_sample_reflections_next_get: (skip)
 679  * @self: the this ptr
 680  * @reflection: the current reflection
 681  *
 682  * Return value: (allow-none): the next reflection or NULL if no more reflection
 683  **/
 684 HklSampleReflection *hkl_sample_reflections_next_get(HklSample *self,
 685                                                      HklSampleReflection *reflection)
 686 {
 687         return list_next(&self->reflections, reflection, list);
 688 }
 689
 690 /**
 691  * hkl_sample_add_reflection:
 692  * @self: the this ptr
 693  * @reflection: The reflection to add
 694  *
 695  * add a reflection to the sample, if the reflection is already part
 696  * of the sample reflection list, this method does nothing.
 697  **/
 698 void hkl_sample_add_reflection(HklSample *self,
 699                                HklSampleReflection *reflection)
 700 {
 701         HklSampleReflection *ref;
 702
 703         list_for_each(&self->reflections, ref, list){
 704                 if (ref == reflection)
 705                         return;
 706         }
 707
 708         list_add_tail(&self->reflections, &reflection->list);
 709         self->n_reflections++;
 710 }
 711
 712 /**
 713  * hkl_sample_del_reflection:
 714  * @self: the this ptr
 715  * @reflection: the reflection to remove.
 716  *
 717  * remove an HklSampleRefelction from the reflections list.
 718  **/
 719 void hkl_sample_del_reflection(HklSample *self,
 720                                HklSampleReflection *reflection)
 721 {
 722         list_del(&reflection->list);
 723         hkl_sample_reflection_free(reflection);
 724         self->n_reflections--;
 725 }
 726
 727 /**
 728  * hkl_sample_compute_UB_busing_levy:
 729  * @self: the this ptr
 730  * @r1: the first #HklsampleReflection
 731  * @r2: the second #HklSampleReflection
 732  *
 733  * compute the UB matrix using the Busing and Levy method
 734  * #todo: add ref
 735  *
 736  * Returns: 0 or 1 if it succeed
 737  **/
 738 int hkl_sample_compute_UB_busing_levy(HklSample *self,
 739                                       const HklSampleReflection *r1,
 740                                       const HklSampleReflection *r2,
 741                                       GError **error)
 742
 743 {
 744         hkl_error (error == NULL || *error == NULL);
 745
 746         if (!hkl_vector_is_colinear(&r1->hkl, &r2->hkl)) {
 747                 HklVector h1c;
 748                 HklVector h2c;
 749                 HklMatrix B;
 750                 HklMatrix Tc;
 751
 752                 /* Compute matrix Tc from r1 and r2. */
 753                 h1c = r1->hkl;
 754                 h2c = r2->hkl;
 755                 hkl_lattice_get_B(self->lattice, &B);
 756                 hkl_matrix_times_vector(&B, &h1c);
 757                 hkl_matrix_times_vector(&B, &h2c);
 758                 hkl_matrix_init_from_two_vector(&Tc, &h1c, &h2c);
 759                 hkl_matrix_transpose(&Tc);
 760
 761                 /* compute U */
 762                 hkl_matrix_init_from_two_vector(&self->U,
 763                                                 &r1->_hkl, &r2->_hkl);
 764                 hkl_matrix_times_matrix(&self->U, &Tc);
 765                 hkl_sample_compute_UxUyUz(self);
 766                 hkl_sample_compute_UB(self);
 767         }else{
 768                 g_set_error(error,
 769                             HKL_SAMPLE_ERROR,
 770                             HKL_SAMPLE_ERROR_COMPUTE_UB_BUSING_LEVY,
 771                             "It is not possible to compute the UB matrix when the given reflections are colinear");
 772
 773                 return FALSE;
 774         }
 775         g_assert (error == NULL || *error == NULL);
 776
 777         return TRUE;
 778 }
 779
 780 /**
 781  * hkl_sample_affine:
 782  * @self: the this ptr
 783  *
 784  * affine the sample
 785  *
 786  * Returns: the fitness of the affined #HklSample
 787  **/
 788 int hkl_sample_affine(HklSample *self, GError **error)
 789 {
 790         return minimize(self, mono_crystal_fitness, self, error);
 791 }
 792
 793 /**
 794  * hkl_sample_get_reflection_measured_angle:
 795  * @self: the this ptr
 796  * @r1: the first #HklSampleReflection
 797  * @r2: the second #HklSampleReflection
 798  *
 799  * get the measured angles between two #HklSampleReflection
 800  *
 801  * Returns: the measured angle beetween them
 802  **/
 803 double hkl_sample_get_reflection_measured_angle(const HklSample *self,
 804                                                 const HklSampleReflection *r1,
 805                                                 const HklSampleReflection *r2)
 806 {
 807         return hkl_vector_angle(&r1->_hkl,
 808                                 &r2->_hkl);
 809 }
 810
 811 /**
 812  * hkl_sample_get_reflection_theoretical_angle:
 813  * @self: the this ptr
 814  * @r1: the first #HklSampleReflection
 815  * @r2: the second #HklSampleReflection
 816  *
 817  * get the theoretical angles between two #HklSampleReflection
 818  *
 819  * Returns: the theoretical angle beetween them
 820  **/
 821 double hkl_sample_get_reflection_theoretical_angle(const HklSample *self,
 822                                                    const HklSampleReflection *r1,
 823                                                    const HklSampleReflection *r2)
 824 {
 825         HklVector hkl1;
 826         HklVector hkl2;
 827
 828         hkl1 = r1->hkl;
 829         hkl2 = r2->hkl;
 830         hkl_matrix_times_vector(&self->UB, &hkl1);
 831         hkl_matrix_times_vector(&self->UB, &hkl2);
 832
 833         return hkl_vector_angle(&hkl1, &hkl2);
 834 }
 835
 836 /**
 837  * hkl_sample_fprintf: (skip)
 838  * @f:
 839  * @self:
 840  *
 841  * print to a file a sample
 842  **/
 843 void hkl_sample_fprintf(FILE *f, const HklSample *self)
 844 {
 845         if(!self)
 846                 return;
 847
 848         fprintf(f, "\nSample name: \"%s\"", self->name);
 849
 850         fprintf(f, "\nLattice parameters:");
 851         fprintf(f, "\n ");
 852         hkl_parameter_fprintf(f, self->lattice->a);
 853         fprintf(f, "\n ");
 854         hkl_parameter_fprintf(f, self->lattice->b);
 855         fprintf(f, "\n ");
 856         hkl_parameter_fprintf(f, self->lattice->c);
 857         fprintf(f, "\n ");
 858         hkl_parameter_fprintf(f, self->lattice->alpha);
 859         fprintf(f, "\n ");
 860         hkl_parameter_fprintf(f, self->lattice->beta);
 861         fprintf(f, "\n ");
 862         hkl_parameter_fprintf(f, self->lattice->gamma);
 863         fprintf(f, "\n ");
 864         hkl_parameter_fprintf(f, self->ux);
 865         fprintf(f, "\n ");
 866         hkl_parameter_fprintf(f, self->uy);
 867         fprintf(f, "\n ");
 868         hkl_parameter_fprintf(f, self->uz);
 869         fprintf(f, "\nUB:\n");
 870         hkl_matrix_fprintf(f, &self->UB);
 871
 872         if (!list_empty(&self->reflections)){
 873                 HklSampleReflection *reflection;
 874                 HklParameter **axis;
 875
 876                 reflection  = list_top(&self->reflections, HklSampleReflection, list);
 877
 878                 fprintf(f, "Reflections:");
 879                 fprintf(f, "\n");
 880                 fprintf(f, "%-10.6s %-10.6s %-10.6s", "h", "k", "l");
 881                 darray_foreach(axis, reflection->geometry->axes){
 882                         fprintf(f, " %-10.6s", (*axis)->name);
 883                 }
 884
 885                 list_for_each(&self->reflections, reflection, list){
 886                         fprintf(f, "\n%-10.6f %-10.6f %-10.6f",
 887                                 reflection->hkl.data[0],
 888                                 reflection->hkl.data[1],
 889                                 reflection->hkl.data[2]);
 890                         darray_foreach(axis, reflection->geometry->axes){
 891                                 fprintf(f, " %-10.6f", hkl_parameter_value_get(*axis, HKL_UNIT_USER));
 892                         }
 893                 }
 894         }
 895 }
 896
 897 /***********************/
 898 /* HklSampleReflection */
 899 /***********************/
 900
 901 /**
 902  * hkl_sample_reflection_new:
 903  * @geometry:
 904  * @detector:
 905  * @h:
 906  * @k:
 907  * @l:
 908  *
 909  * constructeur
 910  *
 911  * Returns:
 912  **/
 913 HklSampleReflection *hkl_sample_reflection_new(const HklGeometry *geometry,
 914                                                const HklDetector *detector,
 915                                                double h, double k, double l,
 916                                                GError **error)
 917 {
 918         HklSampleReflection *self;
 919
 920         if (!geometry || !detector)
 921                 return NULL;
 922
 923         self = HKL_MALLOC(HklSampleReflection);
 924
 925         self->geometry = hkl_geometry_new_copy(geometry);
 926         self->detector = hkl_detector_new_copy(detector);
 927         self->hkl.data[0] = h;
 928         self->hkl.data[1] = k;
 929         self->hkl.data[2] = l;
 930         self->flag = TRUE;
 931
 932         hkl_sample_reflection_update(self);
 933
 934         return self;
 935 }
 936
 937 /**
 938  * hkl_sample_reflection_new_copy: (skip)
 939  * @self:
 940  *
 941  * copy constructor
 942  *
 943  * Returns:
 944  **/
 945 HklSampleReflection *hkl_sample_reflection_new_copy(const HklSampleReflection *self)
 946 {
 947         HklSampleReflection *dup = NULL;
 948
 949         dup = HKL_MALLOC(HklSampleReflection);
 950
 951         dup->geometry = hkl_geometry_new_copy(self->geometry);
 952         dup->detector = hkl_detector_new_copy(self->detector);
 953         dup->hkl = self->hkl;
 954         dup->_hkl = self->_hkl;
 955         dup->flag = self->flag;
 956
 957         return dup;
 958 }
 959
 960 /**
 961  * hkl_sample_reflection_free: (skip)
 962  * @self:
 963  *
 964  * destructor
 965  **/
 966 void hkl_sample_reflection_free(HklSampleReflection *self)
 967 {
 968         hkl_geometry_free(self->geometry);
 969         hkl_detector_free(self->detector);
 970         free(self);
 971 }
 972
 973 /**
 974  * hkl_sample_reflection_hkl_get:
 975  * @self: the this ptr
 976  * @h: (out caller-allocates): the h-coordinate of the #HklSampleReflection
 977  * @k: (out caller-allocates): the k-coordinate of the #HklSampleReflection
 978  * @l: (out caller-allocates): the l-coordinate of the #HklSampleReflection
 979  *
 980  * get the hkl coordinates of the #HklSampleReflection
 981  **/
 982 void hkl_sample_reflection_hkl_get(const HklSampleReflection *self,
 983                                    double *h, double *k, double *l)
 984 {
 985         *h = self->hkl.data[0];
 986         *k = self->hkl.data[1];
 987         *l = self->hkl.data[2];
 988 }
 989
 990 /**
 991  * hkl_sample_reflection_hkl_set:
 992  * @self: the this ptr
 993  * @h: the h-coordinate of the #HklSampleReflection
 994  * @k: the k-coordinate of the #HklSampleReflection
 995  * @l: the l-coordinate of the #HklSampleReflection
 996  * @error: return location for a GError, or NULL
 997  *
 998  * set the hkl coordinates of the #HklSampleReflection
 999  *
1000  * Returns: TRUE on success, FALSE if an error occurred
1001  **/
1002 int hkl_sample_reflection_hkl_set(HklSampleReflection *self,
1003                                   double h, double k, double l,
1004                                   GError **error)
1005 {
1006         hkl_error (error == NULL || *error == NULL);
1007
1008         if((fabs(h) + fabs(k) + fabs(l) < HKL_EPSILON)){
1009                 g_set_error(error,
1010                             HKL_SAMPLE_REFLECTION_ERROR,
1011                             HKL_SAMPLE_REFLECTION_ERROR_HKL_SET,
1012                             "it is not allow to set a null hkl reflection\n");
1013                 return FALSE;
1014         }
1015
1016         self->hkl.data[0] = h;
1017         self->hkl.data[1] = k;
1018         self->hkl.data[2] = l;
1019
1020         return TRUE;
1021 }
1022
1023 /**
1024  * hkl_sample_reflection_flag_get:
1025  * @self: the this ptr
1026  *
1027  * get the flag of the reflection
1028  *
1029  * Returns: the flag value
1030  **/
1031 int hkl_sample_reflection_flag_get(const HklSampleReflection *self)
1032 {
1033         return self->flag;
1034 }
1035
1036 /**
1037  * hkl_sample_reflection_flag_set:
1038  * @self: the this ptr
1039  * @flag: the value of the flag to set
1040  *
1041  * set the flag of the reflection.
1042  **/
1043 void hkl_sample_reflection_flag_set(HklSampleReflection *self, int flag)
1044 {
1045         self->flag = flag;
1046 }
1047
1048 /**
1049  * hkl_sample_reflection_geometry_get:
1050  * @self: the this ptr
1051  *
1052  * get the #HklGeometry stored in the #HklSampleReflection
1053  *
1054  * Returns: the geometry saved into the reflection.
1055  **/
1056 const HklGeometry *hkl_sample_reflection_geometry_get(HklSampleReflection *self)
1057 {
1058         return self->geometry;
1059 }
1060
1061 /**
1062  * hkl_sample_reflection_geometry_set:
1063  * @self: the this ptr
1064  * @geometry: the geometry to set in the  #HklSampleReflection
1065  *
1066  * set the geometry of the reflection
1067  **/
1068 void hkl_sample_reflection_geometry_set(HklSampleReflection *self,
1069                                         const HklGeometry *geometry)
1070 {
1071         if(self->geometry){
1072                 if(self->geometry != geometry){
1073                         hkl_geometry_free(self->geometry);
1074                         self->geometry = hkl_geometry_new_copy(geometry);
1075                 }
1076         }else
1077                 self->geometry = hkl_geometry_new_copy(geometry);
1078
1079         hkl_sample_reflection_update(self);
1080 }