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