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