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", -M_PI, 0., M_PI,
 349                                      TRUE, TRUE,
 350                                      &hkl_unit_angle_rad,
 351                                      &hkl_unit_angle_deg);
 352         self->uy = hkl_parameter_new("uy", -M_PI, 0., M_PI,
 353                                      TRUE, TRUE,
 354                                      &hkl_unit_angle_rad,
 355                                      &hkl_unit_angle_deg);
 356         self->uz = hkl_parameter_new("uz", -M_PI, 0., M_PI,
 357                                      TRUE, TRUE,
 358                                      &hkl_unit_angle_rad,
 359                                      &hkl_unit_angle_deg);
 360
 361         hkl_sample_compute_UB(self);
 362         list_head_init(&self->reflections);
 363         self->n_reflections = 0;
 364
 365         return self;
 366 }
 367
 368 /**
 369  * hkl_sample_new_copy: (skip)
 370  * @self:
 371  *
 372  * copy constructor
 373  *
 374  * Returns:
 375  **/
 376 HklSample *hkl_sample_new_copy(const HklSample *self)
 377 {
 378         HklSample *dup = NULL;
 379
 380         /* check parameters */
 381         if(!self)
 382                 return dup;
 383
 384         dup = HKL_MALLOC(HklSample);
 385
 386         dup->name = strdup(self->name);
 387         dup->lattice = hkl_lattice_new_copy(self->lattice);
 388         dup->U = self->U;
 389         dup->UB = self->UB;
 390         dup->ux = hkl_parameter_new_copy(self->ux);
 391         dup->uy = hkl_parameter_new_copy(self->uy);
 392         dup->uz = hkl_parameter_new_copy(self->uz);
 393
 394         hkl_sample_copy_all_reflections(dup, self);
 395
 396         return dup;
 397 }
 398
 399 /**
 400  * hkl_sample_free: (skip)
 401  * @self:
 402  *
 403  * destructor
 404  **/
 405 void hkl_sample_free(HklSample *self)
 406 {
 407         if (!self)
 408                 return;
 409
 410         free(self->name);
 411         hkl_lattice_free(self->lattice);
 412         hkl_parameter_free(self->ux);
 413         hkl_parameter_free(self->uy);
 414         hkl_parameter_free(self->uz);
 415         hkl_sample_clear_all_reflections(self);
 416         free(self);
 417 }
 418
 419 /**
 420  * hkl_sample_name_get:
 421  * @self: the this ptr
 422  *
 423  * Return value: the name of the sample
 424  **/
 425 const char *hkl_sample_name_get(const HklSample *self)
 426 {
 427         return self->name;
 428 }
 429
 430 /**
 431  * hkl_sample_name_set:
 432  * @self: the this ptr
 433  * @name: the new name to set
 434  *
 435  * set the name of the sample
 436  **/
 437 void hkl_sample_name_set(HklSample *self, const char *name)
 438 {
 439         if(self->name)
 440                 free(self->name);
 441         self->name = strdup(name);
 442 }
 443
 444 /**
 445  * hkl_sample_lattice_get:
 446  * @self: the this ptr
 447  *
 448  * Return value: the lattice parameters of the sample.
 449  **/
 450 const HklLattice *hkl_sample_lattice_get(HklSample *self)
 451 {
 452         return self->lattice;
 453 }
 454
 455 /**
 456  * hkl_sample_lattice_set:
 457  * @self: the this ptr
 458  * @lattice: the lattice to set
 459  **/
 460 void hkl_sample_lattice_set(HklSample *self, const HklLattice *lattice)
 461 {
 462         hkl_lattice_lattice_set(self->lattice, lattice);
 463         hkl_sample_compute_UB(self);
 464 }
 465
 466 /**
 467  * hkl_sample_ux_get:
 468  * @self: the this ptr
 469  *
 470  * Return value: the ux part of the U matrix.
 471  **/
 472 const HklParameter *hkl_sample_ux_get(const HklSample *self)
 473 {
 474         return self->ux;
 475 }
 476
 477 /**
 478  * hkl_sample_uy_get:
 479  * @self: the this ptr
 480  *
 481  * Return value: the uy part of the U matrix.
 482  **/
 483 const HklParameter *hkl_sample_uy_get(const HklSample *self)
 484 {
 485         return self->uy;
 486 }
 487
 488 /**
 489  * hkl_sample_uz_get:
 490  * @self: the this ptr
 491  *
 492  * Return value: the uz part of the U matrix.
 493  **/
 494 const HklParameter *hkl_sample_uz_get(const HklSample *self)
 495 {
 496         return self->uz;
 497 }
 498
 499 /**
 500  * hkl_sample_ux_set:
 501  * @self: the this ptr
 502  * @ux: the ux parameter to set
 503  * @error: return location for a GError, or NULL
 504  *
 505  * set the ux part of the U matrix.
 506  *
 507  * Returns: TRUE on success, FALSE if an error occurred
 508  **/
 509 int hkl_sample_ux_set(HklSample *self,
 510                       const HklParameter *ux,
 511                       GError **error)
 512 {
 513         hkl_error (error == NULL || *error == NULL);
 514
 515         if(!hkl_parameter_init_copy(self->ux, ux, error)){
 516                 g_assert (error == NULL || *error != NULL);
 517                 return FALSE;
 518         }
 519         g_assert (error == NULL || *error == NULL);
 520
 521         hkl_matrix_init_from_euler(&self->U,
 522                                    hkl_parameter_value_get(self->ux, HKL_UNIT_DEFAULT),
 523                                    hkl_parameter_value_get(self->uy, HKL_UNIT_DEFAULT),
 524                                    hkl_parameter_value_get(self->uz, HKL_UNIT_DEFAULT));
 525         hkl_sample_compute_UB(self);
 526
 527         return TRUE;
 528 }
 529
 530 /**
 531  * hkl_sample_uy_set:
 532  * @self: the this ptr
 533  * @uy: the uy parameter to set
 534  * @error: return location for a GError, or NULL
 535  *
 536  * set the uy part of the U matrix.
 537  *
 538  * Returns: TRUE on success, FALSE if an error occurred
 539  **/
 540 int hkl_sample_uy_set(HklSample *self, const HklParameter *uy,
 541                       GError **error)
 542 {
 543         hkl_error (error == NULL || *error == NULL);
 544
 545         if(!hkl_parameter_init_copy(self->uy, uy, error)){
 546                 g_assert (error == NULL || *error != NULL);
 547                 return FALSE;
 548         }
 549         g_assert (error == NULL || *error == NULL);
 550
 551         hkl_matrix_init_from_euler(&self->U,
 552                                    hkl_parameter_value_get(self->ux, HKL_UNIT_DEFAULT),
 553                                    hkl_parameter_value_get(self->uy, HKL_UNIT_DEFAULT),
 554                                    hkl_parameter_value_get(self->uz, HKL_UNIT_DEFAULT));
 555         hkl_sample_compute_UB(self);
 556 }
 557
 558 /**
 559  * hkl_sample_uz_set:
 560  * @self: the this ptr
 561  * @uz: the uz parameter to set
 562  * @error: return location for a GError, or NULL
 563  *
 564  * set the uz part of the U matrix.
 565  *
 566  * Returns: TRUE on success, FALSE if an error occurred
 567  **/
 568 int hkl_sample_uz_set(HklSample *self, const HklParameter *uz,
 569                       GError **error)
 570 {
 571         hkl_error (error == NULL || *error == NULL);
 572
 573         if(!hkl_parameter_init_copy(self->uz, uz, error)){
 574                 g_assert (error == NULL || *error != NULL);
 575                 return FALSE;
 576         }
 577         g_assert (error == NULL || *error == NULL);
 578
 579         hkl_matrix_init_from_euler(&self->U,
 580                                    hkl_parameter_value_get(self->ux, HKL_UNIT_DEFAULT),
 581                                    hkl_parameter_value_get(self->uy, HKL_UNIT_DEFAULT),
 582                                    hkl_parameter_value_get(self->uz, HKL_UNIT_DEFAULT));
 583         hkl_sample_compute_UB(self);
 584 }
 585
 586 /**
 587  * hkl_sample_U_get:
 588  * @self: the this ptr
 589  *
 590  * Return value: the U matrix of the sample
 591  **/
 592 const HklMatrix *hkl_sample_U_get(const HklSample *self)
 593 {
 594         return &self->U;
 595 }
 596
 597 /*
 598  * TODO implemente the error
 599  */
 600 void hkl_sample_U_set(HklSample *self, const HklMatrix *U, GError **error)
 601 {
 602         double x, y, z;
 603
 604         hkl_matrix_matrix_set(&self->U, U);
 605         hkl_sample_compute_UB(self);
 606         hkl_matrix_to_euler(&self->U, &x, &y, &z);
 607         hkl_parameter_value_set(self->ux, x, HKL_UNIT_DEFAULT, NULL);
 608         hkl_parameter_value_set(self->uy, y, HKL_UNIT_DEFAULT, NULL);
 609         hkl_parameter_value_set(self->uz, z, HKL_UNIT_DEFAULT, NULL);
 610 }
 611
 612 /**
 613  * hkl_sample_UB_get:
 614  * @self: the this ptr.
 615  *
 616  * Return value: get the UB matrix of the sample
 617  **/
 618 const HklMatrix *hkl_sample_UB_get(const HklSample *self)
 619 {
 620         return &self->UB;
 621 }
 622
 623 /**
 624  * hkl_sample_UB_set:
 625  * @self: the sample to modify
 626  * @UB: the UB matrix to set
 627  * @error: error set in case of impossibility
 628  *
 629  * Set the UB matrix using an external UB matrix. In fact you give
 630  * the UB matrix but only the U matrix of the sample is affected by
 631  * this operation. We keep the B matrix constant.
 632  * U * B = UB -> U = UB * B^-1
 633  * TODO implemente the error
 634  **/
 635 int hkl_sample_UB_set(HklSample *self, const HklMatrix *UB,
 636                       GError **error)
 637 {
 638         struct set_UB_t params = {
 639                 .sample = self,
 640                 .UB = *UB
 641         };
 642
 643         return minimize(self, set_UB_fitness, &params, error);
 644 }
 645
 646 /**
 647  * hkl_sample_n_reflections_get: (skip)
 648  * @self: the this ptr
 649  *
 650  * return the number of reflections of the sample
 651  *
 652  * Returns:
 653  **/
 654 size_t hkl_sample_n_reflections_get(const HklSample *self)
 655 {
 656         return self->n_reflections;
 657 }
 658
 659 /**
 660  * hkl_sample_reflections_first_get: (skip)
 661  * @self: the this ptr
 662  *
 663  * Return value: the first HklSampleReflection of the sample.
 664  **/
 665 HklSampleReflection *hkl_sample_reflections_first_get(HklSample *self)
 666 {
 667         return list_top(&self->reflections, HklSampleReflection, list);
 668 }
 669
 670 /**
 671  * hkl_sample_reflections_next_get: (skip)
 672  * @self: the this ptr
 673  * @reflection: the current reflection
 674  *
 675  * Return value: (allow-none): the next reflection or NULL if no more reflection
 676  **/
 677 HklSampleReflection *hkl_sample_reflections_next_get(HklSample *self,
 678                                                      HklSampleReflection *reflection)
 679 {
 680         return list_next(&self->reflections, reflection, list);
 681 }
 682
 683 /**
 684  * hkl_sample_add_reflection:
 685  * @self: the this ptr
 686  * @reflection: The reflection to add
 687  *
 688  * add a reflection to the sample, if the reflection is already part
 689  * of the sample reflection list, this method does nothing.
 690  **/
 691 void hkl_sample_add_reflection(HklSample *self,
 692                                HklSampleReflection *reflection)
 693 {
 694         HklSampleReflection *ref;
 695
 696         list_for_each(&self->reflections, ref, list){
 697                 if (ref == reflection)
 698                         return;
 699         }
 700
 701         list_add_tail(&self->reflections, &reflection->list);
 702         self->n_reflections++;
 703 }
 704
 705 /**
 706  * hkl_sample_del_reflection:
 707  * @self: the this ptr
 708  * @reflection: the reflection to remove.
 709  *
 710  * remove an HklSampleRefelction from the reflections list.
 711  **/
 712 void hkl_sample_del_reflection(HklSample *self,
 713                                HklSampleReflection *reflection)
 714 {
 715         list_del(&reflection->list);
 716         hkl_sample_reflection_free(reflection);
 717         self->n_reflections--;
 718 }
 719
 720 /**
 721  * hkl_sample_compute_UB_busing_levy:
 722  * @self: the this ptr
 723  * @r1: the first #HklsampleReflection
 724  * @r2: the second #HklSampleReflection
 725  *
 726  * compute the UB matrix using the Busing and Levy method
 727  * #todo: add ref
 728  *
 729  * Returns: 0 or 1 if it succeed
 730  **/
 731 int hkl_sample_compute_UB_busing_levy(HklSample *self,
 732                                       const HklSampleReflection *r1,
 733                                       const HklSampleReflection *r2,
 734                                       GError **error)
 735
 736 {
 737         hkl_error (error == NULL || *error == NULL);
 738
 739         if (!hkl_vector_is_colinear(&r1->hkl, &r2->hkl)) {
 740                 HklVector h1c;
 741                 HklVector h2c;
 742                 HklMatrix B;
 743                 HklMatrix Tc;
 744
 745                 /* Compute matrix Tc from r1 and r2. */
 746                 h1c = r1->hkl;
 747                 h2c = r2->hkl;
 748                 hkl_lattice_get_B(self->lattice, &B);
 749                 hkl_matrix_times_vector(&B, &h1c);
 750                 hkl_matrix_times_vector(&B, &h2c);
 751                 hkl_matrix_init_from_two_vector(&Tc, &h1c, &h2c);
 752                 hkl_matrix_transpose(&Tc);
 753
 754                 /* compute U */
 755                 hkl_matrix_init_from_two_vector(&self->U,
 756                                                 &r1->_hkl, &r2->_hkl);
 757                 hkl_matrix_times_matrix(&self->U, &Tc);
 758                 hkl_sample_compute_UxUyUz(self);
 759                 hkl_sample_compute_UB(self);
 760         }else{
 761                 g_set_error(error,
 762                             HKL_SAMPLE_ERROR,
 763                             HKL_SAMPLE_ERROR_COMPUTE_UB_BUSING_LEVY,
 764                             "It is not possible to compute the UB matrix when the given reflections are colinear");
 765
 766                 return FALSE;
 767         }
 768         g_assert (error == NULL || *error == NULL);
 769
 770         return TRUE;
 771 }
 772
 773 /**
 774  * hkl_sample_affine:
 775  * @self: the this ptr
 776  *
 777  * affine the sample
 778  *
 779  * Returns: the fitness of the affined #HklSample
 780  **/
 781 int hkl_sample_affine(HklSample *self, GError **error)
 782 {
 783         return minimize(self, mono_crystal_fitness, self, error);
 784 }
 785
 786 /**
 787  * hkl_sample_get_reflection_mesured_angle:
 788  * @self: the this ptr
 789  * @r1: the first #HklSampleReflection
 790  * @r2: the second #HklSampleReflection
 791  *
 792  * get the mesured angles between two #HklSampleReflection
 793  *
 794  * Returns: the mesured angle beetween them
 795  **/
 796 double hkl_sample_get_reflection_mesured_angle(const HklSample *self,
 797                                                const HklSampleReflection *r1,
 798                                                const HklSampleReflection *r2)
 799 {
 800         return hkl_vector_angle(&r1->_hkl,
 801                                 &r2->_hkl);
 802 }
 803
 804 /**
 805  * hkl_sample_get_reflection_theoretical_angle:
 806  * @self: the this ptr
 807  * @r1: the first #HklSampleReflection
 808  * @r2: the second #HklSampleReflection
 809  *
 810  * get the theoretical angles between two #HklSampleReflection
 811  *
 812  * Returns: the theoretical angle beetween them
 813  **/
 814 double hkl_sample_get_reflection_theoretical_angle(const HklSample *self,
 815                                                    const HklSampleReflection *r1,
 816                                                    const HklSampleReflection *r2)
 817 {
 818         HklVector hkl1;
 819         HklVector hkl2;
 820
 821         hkl1 = r1->hkl;
 822         hkl2 = r2->hkl;
 823         hkl_matrix_times_vector(&self->UB, &hkl1);
 824         hkl_matrix_times_vector(&self->UB, &hkl2);
 825
 826         return hkl_vector_angle(&hkl1, &hkl2);
 827 }
 828
 829 /**
 830  * hkl_sample_fprintf: (skip)
 831  * @f:
 832  * @self:
 833  *
 834  * print to a file a sample
 835  **/
 836 void hkl_sample_fprintf(FILE *f, const HklSample *self)
 837 {
 838         if(!self)
 839                 return;
 840
 841         fprintf(f, "\nSample name: \"%s\"", self->name);
 842
 843         fprintf(f, "\nLattice parameters:");
 844         fprintf(f, "\n ");
 845         hkl_parameter_fprintf(f, self->lattice->a);
 846         fprintf(f, "\n ");
 847         hkl_parameter_fprintf(f, self->lattice->b);
 848         fprintf(f, "\n ");
 849         hkl_parameter_fprintf(f, self->lattice->c);
 850         fprintf(f, "\n ");
 851         hkl_parameter_fprintf(f, self->lattice->alpha);
 852         fprintf(f, "\n ");
 853         hkl_parameter_fprintf(f, self->lattice->beta);
 854         fprintf(f, "\n ");
 855         hkl_parameter_fprintf(f, self->lattice->gamma);
 856         fprintf(f, "\n ");
 857         hkl_parameter_fprintf(f, self->ux);
 858         fprintf(f, "\n ");
 859         hkl_parameter_fprintf(f, self->uy);
 860         fprintf(f, "\n ");
 861         hkl_parameter_fprintf(f, self->uz);
 862         fprintf(f, "\nUB:\n");
 863         hkl_matrix_fprintf(f, &self->UB);
 864
 865         if (!list_empty(&self->reflections)){
 866                 HklSampleReflection *reflection;
 867                 HklParameter **axis;
 868
 869                 reflection  = list_top(&self->reflections, HklSampleReflection, list);
 870
 871                 fprintf(f, "Reflections:");
 872                 fprintf(f, "\n");
 873                 fprintf(f, "%-10.6s %-10.6s %-10.6s", "h", "k", "l");
 874                 darray_foreach(axis, reflection->geometry->axes){
 875                         fprintf(f, " %-10.6s", (*axis)->name);
 876                 }
 877
 878                 list_for_each(&self->reflections, reflection, list){
 879                         fprintf(f, "\n%-10.6f %-10.6f %-10.6f",
 880                                 reflection->hkl.data[0],
 881                                 reflection->hkl.data[1],
 882                                 reflection->hkl.data[2]);
 883                         darray_foreach(axis, reflection->geometry->axes){
 884                                 fprintf(f, " %-10.6f", hkl_parameter_value_get(*axis, HKL_UNIT_USER));
 885                         }
 886                 }
 887         }
 888 }
 889
 890 /***********************/
 891 /* HklSampleReflection */
 892 /***********************/
 893
 894 /**
 895  * hkl_sample_reflection_new:
 896  * @geometry:
 897  * @detector:
 898  * @h:
 899  * @k:
 900  * @l:
 901  *
 902  * constructeur
 903  *
 904  * Returns:
 905  **/
 906 HklSampleReflection *hkl_sample_reflection_new(const HklGeometry *geometry,
 907                                                const HklDetector *detector,
 908                                                double h, double k, double l,
 909                                                GError **error)
 910 {
 911         HklSampleReflection *self;
 912
 913         if (!geometry || !detector)
 914                 return NULL;
 915
 916         self = HKL_MALLOC(HklSampleReflection);
 917
 918         self->geometry = hkl_geometry_new_copy(geometry);
 919         self->detector = hkl_detector_new_copy(detector);
 920         self->hkl.data[0] = h;
 921         self->hkl.data[1] = k;
 922         self->hkl.data[2] = l;
 923         self->flag = TRUE;
 924
 925         hkl_sample_reflection_update(self);
 926
 927         return self;
 928 }
 929
 930 /**
 931  * hkl_sample_reflection_new_copy: (skip)
 932  * @self:
 933  *
 934  * copy constructor
 935  *
 936  * Returns:
 937  **/
 938 HklSampleReflection *hkl_sample_reflection_new_copy(const HklSampleReflection *self)
 939 {
 940         HklSampleReflection *dup = NULL;
 941
 942         dup = HKL_MALLOC(HklSampleReflection);
 943
 944         dup->geometry = hkl_geometry_new_copy(self->geometry);
 945         dup->detector = hkl_detector_new_copy(self->detector);
 946         dup->hkl = self->hkl;
 947         dup->_hkl = self->_hkl;
 948         dup->flag = self->flag;
 949
 950         return dup;
 951 }
 952
 953 /**
 954  * hkl_sample_reflection_free: (skip)
 955  * @self:
 956  *
 957  * destructor
 958  **/
 959 void hkl_sample_reflection_free(HklSampleReflection *self)
 960 {
 961         hkl_geometry_free(self->geometry);
 962         hkl_detector_free(self->detector);
 963         free(self);
 964 }
 965
 966 /**
 967  * hkl_sample_reflection_hkl_get:
 968  * @self: the this ptr
 969  * @h: (out caller-allocates): the h-coordinate of the #HklSampleReflection
 970  * @k: (out caller-allocates): the k-coordinate of the #HklSampleReflection
 971  * @l: (out caller-allocates): the l-coordinate of the #HklSampleReflection
 972  *
 973  * get the hkl coordinates of the #HklSampleReflection
 974  **/
 975 void hkl_sample_reflection_hkl_get(const HklSampleReflection *self,
 976                                    double *h, double *k, double *l)
 977 {
 978         *h = self->hkl.data[0];
 979         *k = self->hkl.data[1];
 980         *l = self->hkl.data[2];
 981 }
 982
 983 /**
 984  * hkl_sample_reflection_hkl_set:
 985  * @self: the this ptr
 986  * @h: the h-coordinate of the #HklSampleReflection
 987  * @k: the k-coordinate of the #HklSampleReflection
 988  * @l: the l-coordinate of the #HklSampleReflection
 989  * @error: return location for a GError, or NULL
 990  *
 991  * set the hkl coordinates of the #HklSampleReflection
 992  *
 993  * Returns: TRUE on success, FALSE if an error occurred
 994  **/
 995 int hkl_sample_reflection_hkl_set(HklSampleReflection *self,
 996                                   double h, double k, double l,
 997                                   GError **error)
 998 {
 999         hkl_error (error == NULL || *error == NULL);
1000
1001         if((fabs(h) + fabs(k) + fabs(l) < HKL_EPSILON)){
1002                 g_set_error(error,
1003                             HKL_SAMPLE_REFLECTION_ERROR,
1004                             HKL_SAMPLE_REFLECTION_ERROR_HKL_SET,
1005                             "it is not allow to set a null hkl reflection\n");
1006                 return FALSE;
1007         }
1008
1009         self->hkl.data[0] = h;
1010         self->hkl.data[1] = k;
1011         self->hkl.data[2] = l;
1012
1013         return TRUE;
1014 }
1015
1016 /**
1017  * hkl_sample_reflection_flag_get: (skip)
1018  * @self:
1019  *
1020  * get the flag of the reflection
1021  **/
1022 int hkl_sample_reflection_flag_get(const HklSampleReflection *self)
1023 {
1024         return self->flag;
1025 }
1026
1027 /**
1028  * hkl_sample_reflection_flag_set: (skip)
1029  * @self:
1030  * @flag:
1031  *
1032  * set the flag of the reglection
1033  **/
1034 void hkl_sample_reflection_flag_set(HklSampleReflection *self, int flag)
1035 {
1036         self->flag = flag;
1037 }
1038
1039 /**
1040  * hkl_sample_reflection_geometry_get: (skip)
1041  * @self:
1042  *
1043  * set the geometry of the reflection
1044  **/
1045 const HklGeometry *hkl_sample_reflection_geometry_get(HklSampleReflection *self)
1046 {
1047         return self->geometry;
1048 }
1049
1050 /**
1051  * hkl_sample_reflection_geometry_set: (skip)
1052  * @self:
1053  * @geometry:
1054  *
1055  * set the geometry of the reflection
1056  **/
1057 void hkl_sample_reflection_geometry_set(HklSampleReflection *self,
1058                                         const HklGeometry *geometry)
1059 {
1060         if(self->geometry){
1061                 if(self->geometry != geometry){
1062                         hkl_geometry_free(self->geometry);
1063                         self->geometry = hkl_geometry_new_copy(geometry);
1064                 }
1065         }else
1066                 self->geometry = hkl_geometry_new_copy(geometry);
1067
1068         hkl_sample_reflection_update(self);
1069 }