hkl/hkl-sample.c

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