hkl/hkl-geometry.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-2010 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 #include <math.h>
  23 #include <string.h>
  24 #include <stdarg.h>
  25 #ifndef _MSC_VER
  26 # include <alloca.h>
  27 #endif
  28
  29 #include <gsl/gsl_math.h>
  30 #include <gsl/gsl_sf_trig.h>
  31
  32 #include <hkl/hkl-geometry.h>
  33
  34 /*
  35  * Try to add a axis to the axes list,
  36  * if a identical axis is present in the list return it
  37  * else create a new on and add it to the list.
  38  * die if try to add an axis with the same name but a different axis_v
  39  */
  40 static size_t hkl_geometry_add_rotation(HklGeometry *self,
  41                                         const char *name, const HklVector *axis_v)
  42 {
  43         size_t i, len;
  44         HklAxis axis;
  45
  46         /* check if an axis with the same name is on the axis list */
  47         for(i=0; i<self->len; ++i){
  48                 HklAxis *axis;
  49
  50                 axis = &self->axes[i];
  51                 if(!strcmp(hkl_axis_get_name(axis), name)){
  52                         if (hkl_vector_cmp(&axis->axis_v, axis_v)){
  53                                 fprintf(stderr, "can not add two axis with the same name \"%s\" but different axes <%f, %f, %f> != <%f, %f, %f> into an HklAxes.",
  54                                         name,
  55                                         axis->axis_v.data[0], axis->axis_v.data[1], axis->axis_v.data[2],
  56                                         axis_v->data[0], axis_v->data[1], axis_v->data[2]);
  57                                 exit(128);
  58                         }else{
  59                                 return i;
  60                         }
  61                 }
  62         }
  63
  64         /* no so create and add it to the list */
  65         self->axes = realloc(self->axes, sizeof(*self->axes) * (self->len + 1));
  66         hkl_axis_init(&self->axes[self->len], name, axis_v);
  67
  68         return self->len++;
  69 }
  70
  71 /*******************/
  72 /* HklHolderConfig */
  73 /*******************/
  74
  75 static struct HklHolderConfig *hkl_holder_config_new(void)
  76 {
  77         struct HklHolderConfig *self;
  78
  79         self = HKL_MALLOC(struct HklHolderConfig);
  80
  81         self->gc = 1;
  82         self->idx = NULL;
  83         self->len = 0;
  84
  85         return self;
  86 }
  87
  88 static struct HklHolderConfig *hkl_holder_config_ref(struct HklHolderConfig *self)
  89 {
  90         if(!self)
  91                 return NULL;
  92
  93         self->gc++;
  94
  95         return self;
  96 }
  97
  98 static void hkl_holder_config_unref(struct HklHolderConfig *self)
  99 {
 100         if(!self)
 101                 return;
 102
 103         if(--self->gc)
 104                 return;
 105
 106         free(self->idx);
 107         free(self);
 108 }
 109
 110 /*************/
 111 /* HklHolder */
 112 /*************/
 113
 114 static void hkl_holder_init(HklHolder *self, HklGeometry *geometry)
 115 {
 116         static HklQuaternion q0 = {{0, 0, 0, 0}};
 117         self->config = hkl_holder_config_new();
 118         self->geometry = geometry;
 119         self->q = q0;
 120 }
 121
 122 static int hkl_holder_init_copy(HklHolder *self, HklGeometry *geometry,
 123                                 const HklHolder *holder)
 124 {
 125         /* check axes compatibility */
 126         if (geometry->len != holder->geometry->len)
 127                 return HKL_FALSE;
 128
 129         self->config = hkl_holder_config_ref(holder->config);
 130         self->geometry = geometry;
 131         self->q = holder->q;
 132
 133         return HKL_TRUE;
 134 }
 135
 136 static void hkl_holder_release_memory(HklHolder *self)
 137 {
 138         hkl_holder_config_unref(self->config);
 139 }
 140
 141 static void hkl_holder_update(HklHolder *self)
 142 {
 143         static HklQuaternion q0 = {{1, 0, 0, 0}};
 144         size_t i;
 145         HklAxis *axes;
 146         size_t *idx;
 147
 148         self->q = q0;
 149         axes = self->geometry->axes;
 150         idx = self->config->idx;
 151         for(i=0; i<self->config->len; ++i)
 152                 hkl_quaternion_times_quaternion(&self->q, &axes[idx[i]].q);
 153 }
 154
 155 HklAxis *hkl_holder_add_rotation_axis(HklHolder *self,
 156                                       const char *name, double x, double y, double z)
 157 {
 158         HklAxis *axis = NULL;
 159         size_t i, idx;
 160         HklVector axis_v;
 161
 162         axis_v.data[0] = x;
 163         axis_v.data[1] = y;
 164         axis_v.data[2] = z;
 165
 166         idx = hkl_geometry_add_rotation(self->geometry, name, &axis_v);
 167
 168         /* check that the axis is not already in the holder */
 169         for(i=0; i<self->config->len; i++)
 170                 if (idx == self->config->idx[i])
 171                         return NULL;
 172
 173         axis = &self->geometry->axes[idx];
 174         self->config->idx = realloc(self->config->idx, sizeof(*self->config->idx) * (self->config->len + 1));
 175         self->config->idx[self->config->len++] = idx;
 176
 177         return axis;
 178 }
 179
 180 /***************/
 181 /* HklGeometry */
 182 /***************/
 183
 184 /**
 185  * hkl_geometry_new: (skip)
 186  *
 187  * constructor
 188  *
 189  * Returns:
 190  **/
 191 HklGeometry *hkl_geometry_new(void)
 192 {
 193         HklGeometry *g = NULL;
 194
 195         g = HKL_MALLOC(HklGeometry);
 196
 197         g->config = NULL;
 198         hkl_source_init(&g->source, 1.54, 1, 0, 0);
 199         g->axes = NULL;
 200         g->len = 0;
 201         g->holders = NULL;
 202         g->holders_len = 0;
 203
 204         return g;
 205 }
 206
 207 /**
 208  * hkl_geometry_new_copy: (skip)
 209  * @self:
 210  *
 211  * copy constructor
 212  *
 213  * Returns:
 214  **/
 215 HklGeometry *hkl_geometry_new_copy(const HklGeometry *self)
 216 {
 217         HklGeometry *dup = NULL;
 218         unsigned int i;
 219         size_t len;
 220
 221         if(!self)
 222                 return dup;
 223
 224         dup = HKL_MALLOC(HklGeometry);
 225
 226         dup->config = self->config;
 227         dup->source = self->source;
 228         dup->len = self->len;
 229         dup->holders_len = self->holders_len;
 230
 231         /* copy the axes */
 232         dup->axes = malloc(sizeof(*dup->axes) * dup->len);
 233         memcpy(dup->axes, self->axes, sizeof(*dup->axes) * dup->len);
 234
 235         /* copy the holders */
 236         dup->holders = malloc(sizeof(*dup->holders) * dup->holders_len);
 237         for(i=0; i<dup->holders_len; ++i)
 238                 hkl_holder_init_copy(&dup->holders[i], dup,
 239                                      &self->holders[i]);
 240
 241         return dup;
 242 }
 243
 244 /**
 245  * hkl_geometry_free: (skip)
 246  * @self:
 247  *
 248  * destructor
 249  **/
 250 void hkl_geometry_free(HklGeometry *self)
 251 {
 252         size_t i;
 253
 254         free(self->axes);
 255         self->len = 0;
 256
 257         if(self->holders_len) {
 258                 for(i=0; i<self->holders_len; ++i)
 259                         hkl_holder_release_memory(&self->holders[i]);
 260                 free(self->holders);
 261         }
 262
 263         free(self);
 264 }
 265
 266 /**
 267  * hkl_geometry_init_geometry: (skip)
 268  * @self:
 269  * @src:
 270  *
 271  * initilize an HklGeometry
 272  **/
 273 void hkl_geometry_init_geometry(HklGeometry *self, const HklGeometry *src)
 274 {
 275         size_t i;
 276
 277         if(!self || !src
 278            || self->config->type != src->config->type)
 279                 return;
 280
 281         self->source = src->source;
 282
 283         /* copy the axes configuration and mark it as dirty */
 284         memcpy(self->axes, src->axes, sizeof(*self->axes) * self->len);
 285         for(i=0; i<src->holders_len; ++i)
 286                 self->holders[i].q = src->holders[i].q;
 287 }
 288
 289 /**
 290  * hkl_geometry_add_holder: (skip)
 291  * @self:
 292  *
 293  * add an Holder to the #HklGeometry
 294  *
 295  * Returns:
 296  **/
 297 HklHolder *hkl_geometry_add_holder(HklGeometry *self)
 298 {
 299         self->holders = realloc(self->holders, sizeof(*self->holders) * (self->holders_len + 1));
 300         hkl_holder_init(&self->holders[self->holders_len], self);
 301
 302         return &self->holders[self->holders_len++];
 303 }
 304
 305 /**
 306  * hkl_geometry_update: (skip)
 307  * @self:
 308  *
 309  * update the geometry internal once an Axis values changed
 310  **/
 311 void hkl_geometry_update(HklGeometry *self)
 312 {
 313         size_t i;
 314         int ko = 0;
 315
 316         for(i=0; i<self->len; ++i)
 317                 if (hkl_axis_get_changed(&self->axes[i])) {
 318                         ko = 1;
 319                         break;
 320                 }
 321
 322         if (ko) {
 323                 for(i=0; i<self->holders_len; i++)
 324                         hkl_holder_update(&self->holders[i]);
 325
 326                 for(i=0; i<self->len; i++)
 327                         hkl_axis_set_changed(&self->axes[i], HKL_FALSE);
 328         }
 329 }
 330
 331 /**
 332  * hkl_geometry_get_axis_idx_by_name: (skip)
 333  * @self:
 334  * @name:
 335  *
 336  * get the index of the axes named @name in the geometry
 337  *
 338  * Returns: -1 if the axis was not found
 339  **/
 340 int hkl_geometry_get_axis_idx_by_name(HklGeometry *self, const char *name)
 341 {
 342         uint i;
 343         HklAxis *axis;
 344
 345         if (!self || !name)
 346                 return -1;
 347
 348         for(i=0; i<self->len; ++i){
 349                 axis = &self->axes[i];
 350                 if (!strcmp(hkl_axis_get_name(axis), name))
 351                         return i;
 352         }
 353
 354         return -1;
 355 }
 356
 357 /**
 358  * hkl_geometry_get_axis_by_name:
 359  * @self:
 360  * @name:
 361  *
 362  * get an #HklAxis using its name
 363  *
 364  * Returns: (transfer none):
 365  **/
 366 HklAxis *hkl_geometry_get_axis_by_name(HklGeometry *self, const char *name)
 367 {
 368         size_t i;
 369         HklAxis *axis;
 370         for(i=0; i<self->len; ++i) {
 371                 axis = &self->axes[i];
 372                 if (!strcmp(hkl_axis_get_name(axis), name))
 373                         return axis;
 374         }
 375         return NULL;
 376 }
 377
 378 /**
 379  * hkl_geometry_randomize: (skip)
 380  * @self:
 381  *
 382  * randomize the #HklGeometry
 383  **/
 384 void hkl_geometry_randomize(HklGeometry *self)
 385 {
 386         size_t i;
 387
 388         for(i=0; i<self->len; ++i)
 389                 hkl_axis_randomize(&self->axes[i]);
 390         hkl_geometry_update(self);
 391 }
 392
 393 /**
 394  * hkl_geometry_set_values_v: (skip)
 395  * @self:
 396  * @len:
 397  * "...:
 398  *
 399  * set the axes values
 400  *
 401  * Returns:
 402  **/
 403 int hkl_geometry_set_values_v(HklGeometry *self, size_t len, ...)
 404 {
 405         va_list ap;
 406         size_t i;
 407
 408         if (!self || len != self->len)
 409                 return HKL_FALSE;
 410
 411         va_start(ap, len);
 412         for(i=0; i<len; ++i)
 413                 hkl_axis_set_value(&self->axes[i], va_arg(ap, double));
 414
 415         va_end(ap);
 416         hkl_geometry_update(self);
 417
 418         return HKL_TRUE;
 419 }
 420
 421 int hkl_geometry_set_values_unit_v(HklGeometry *self, ...)
 422 {
 423         va_list ap;
 424         size_t i;
 425
 426         if (!self)
 427                 return HKL_FALSE;
 428
 429         va_start(ap, self);
 430         for(i=0; i<self->len; ++i)
 431                 hkl_axis_set_value_unit(&self->axes[i], va_arg(ap, double));
 432
 433         va_end(ap);
 434         hkl_geometry_update(self);
 435
 436         return HKL_TRUE;
 437 }
 438
 439 /**
 440  * hkl_geometry_distance: (skip)
 441  * @self:
 442  * @geom:
 443  *
 444  * compute the distance between two #HklGeometries
 445  *
 446  * Returns:
 447  **/
 448 double hkl_geometry_distance(HklGeometry *self, HklGeometry *geom)
 449 {
 450         size_t i;
 451         double value1, value2;
 452         double distance = 0.;
 453
 454         if (!self || !geom)
 455                 return 0.;
 456
 457         for(i=0; i<self->len; ++i){
 458                 value1 = hkl_axis_get_value(&self->axes[i]);
 459                 value2 = hkl_axis_get_value(&geom->axes[i]);
 460                 distance += fabs(value2 - value1);
 461         }
 462
 463         return distance;
 464 }
 465
 466 /**
 467  * hkl_geometry_distance_orthodromic: (skip)
 468  * @self:
 469  * @geom:
 470  *
 471  * compute the orthodromique distance between two #HklGeometry
 472  *
 473  * Returns:
 474  **/
 475 double hkl_geometry_distance_orthodromic(HklGeometry *self, HklGeometry *geom)
 476 {
 477         size_t i;
 478         double value1, value2;
 479         double distance = 0.;
 480
 481         if (!self || !geom)
 482                 return 0.;
 483
 484         for(i=0; i<self->len; ++i){
 485                 double d;
 486
 487                 value1 = hkl_axis_get_value(&self->axes[i]);
 488                 value2 = hkl_axis_get_value(&geom->axes[i]);
 489                 d = fabs(gsl_sf_angle_restrict_symm(value2) - gsl_sf_angle_restrict_symm(value1));
 490                 /* as M_PI and -M_PI are included in the GSL restriction */
 491                 if (d > M_PI)
 492                         d = 2*M_PI - d;
 493                 distance += d;
 494         }
 495
 496         return distance;
 497 }
 498
 499 /**
 500  * hkl_geometry_is_valid: (skip)
 501  * @self:
 502  *
 503  * check if all axes of the #HklGeometry are valid.
 504  *
 505  * Returns:
 506  **/
 507 int hkl_geometry_is_valid(const HklGeometry *self)
 508 {
 509         size_t i;
 510
 511         for(i=0; i<self->len; ++i)
 512                 if(hkl_axis_is_valid(&self->axes[i]) == HKL_FALSE)
 513                         return HKL_FALSE;
 514
 515         return HKL_TRUE;
 516 }
 517
 518 /**
 519  * hkl_geometry_closest_from_geometry_with_range: (skip)
 520  * @self:
 521  * @ref:
 522  *
 523  * get the closest axes values in the HklInterval compatible with the
 524  * current axes values
 525  *
 526  * Returns:
 527  **/
 528 int hkl_geometry_closest_from_geometry_with_range(HklGeometry *self, HklGeometry *ref)
 529 {
 530         size_t i;
 531         double *values = alloca(self->len * sizeof(*values));
 532         int ko = HKL_FALSE;
 533
 534         for(i=0;i<self->len;++i){
 535                 values[i] = hkl_axis_get_value_closest(&self->axes[i], &ref->axes[i]);
 536                 if(gsl_isnan(values[i])){
 537                         ko = HKL_TRUE;
 538                         break;
 539                 }
 540         }
 541         if(!ko){
 542                 for(i=0;i<self->len;++i)
 543                         hkl_axis_set_value(&self->axes[i], values[i]);
 544                 hkl_geometry_update(self);
 545         }
 546         return ko;
 547 }
 548
 549 /**
 550  * hkl_geometry_get_axes_values_unit:
 551  * @self:
 552  * @len: (out caller-allocates)
 553  *
 554  * return all the axes values (must be free by the user)
 555  *
 556  * Returns: (array length=len) (transfer full):
 557  **/
 558 double *hkl_geometry_get_axes_values_unit(const HklGeometry *self, unsigned int *len)
 559 {
 560         double *values;
 561         uint i;
 562
 563         if(!self || !len || self->len == 0)
 564                 return NULL;
 565
 566         *len = self->len;
 567         values = malloc(self->len * sizeof(*values));
 568         if(!values)
 569                 return NULL;
 570         for(i=0; i<self->len; ++i)
 571                 values[i] = hkl_axis_get_value_unit(&self->axes[i]);
 572
 573         return values;
 574 }
 575
 576 /**
 577  * hkl_geometry_set_axes_values_unit:
 578  * @self:
 579  * @values: (array length=len):
 580  * @len:
 581  *
 582  * set the axes values
 583  **/
 584 void hkl_geometry_set_axes_values_unit(HklGeometry *self, double *values, unsigned int len)
 585 {
 586         uint i;
 587
 588         if (!self || !values || len != self->len)
 589                 return;
 590
 591         for(i=0; i<self->len; ++i)
 592                 hkl_axis_set_value_unit(&self->axes[i], values[i]);
 593         hkl_geometry_update(self);
 594 }
 595
 596 /**
 597  * hkl_geometry_fprintf: (skip)
 598  * @file:
 599  * @self:
 600  *
 601  * print into a file the #HklGeometry
 602  **/
 603 void hkl_geometry_fprintf(FILE *file, const HklGeometry *self)
 604 {
 605         size_t i;
 606
 607         if(!self)
 608                 return;
 609
 610         for(i=0; i<self->len; ++i){
 611                 if(i)
 612                         fprintf(file, "\n");
 613                 hkl_parameter_fprintf(file, (HklParameter *)(&self->axes[i]));
 614         }
 615 }
 616
 617 /*******************/
 618 /* HklGeometryList */
 619 /*******************/
 620
 621 /**
 622  * hkl_geometry_list_new: (skip)
 623  *
 624  * constructor
 625  *
 626  * Returns:
 627  **/
 628 HklGeometryList *hkl_geometry_list_new(void)
 629 {
 630         HklGeometryList *self;
 631
 632         self = HKL_MALLOC(HklGeometryList);
 633
 634         list_head_init(&self->items);
 635         self->len = 0;
 636         self->multiply = NULL;
 637
 638         return self;
 639 }
 640
 641 /**
 642  * hkl_geometry_list_new_copy: (skip)
 643  * @self:
 644  *
 645  * copy constructor
 646  *
 647  * Returns:
 648  **/
 649 HklGeometryList *hkl_geometry_list_new_copy(const HklGeometryList *self)
 650 {
 651         HklGeometryList *dup;
 652         HklGeometryListItem *item;
 653         int i;
 654
 655         if (!self)
 656                 return NULL;
 657
 658         dup = HKL_MALLOC(HklGeometryList);
 659
 660         *dup = *self;
 661         list_head_init(&dup->items);
 662
 663         /* now copy the item arrays */
 664         list_for_each(&self->items, item, node){
 665                 HklGeometryListItem *dup_item = hkl_geometry_list_item_new_copy(item);
 666                 list_add_tail(&dup->items, &dup_item->node);
 667         }
 668
 669         return dup;
 670 }
 671
 672 /**
 673  * hkl_geometry_list_free: (skip)
 674  * @self:
 675  *
 676  * destructor
 677  **/
 678 void hkl_geometry_list_free(HklGeometryList *self)
 679 {
 680         hkl_geometry_list_reset(self);
 681         free(self);
 682 }
 683
 684 /**
 685  * @brief: (skip)
 686  * @self: The current #HklGeometryList
 687  * @geometry: the #HklGeometry to add
 688  *
 689  * this method Add a geometry to the geometries
 690  *
 691  * This method try to be clever by allocating memory only if the
 692  * current length of the geometries is not large enought. Then it just
 693  * set the geometry axes and copy it to the right geometries. We do
 694  * not gives the x len as it is equal to the self->axes_len.
 695  **/
 696 void hkl_geometry_list_add(HklGeometryList *self, HklGeometry *geometry)
 697 {
 698         int ko;
 699         HklGeometryListItem *item;
 700
 701         /* now check if the geometry is already in the geometry list */
 702         ko = HKL_FALSE;
 703         list_for_each(&self->items, item, node)
 704                 if (hkl_geometry_distance_orthodromic(geometry,
 705                                                       item->geometry) < HKL_EPSILON)
 706                         ko = HKL_TRUE;
 707
 708         if(ko == HKL_FALSE){
 709                 item = hkl_geometry_list_item_new(geometry);
 710                 list_add_tail(&self->items, &item->node);
 711                 self->len++;
 712         }
 713 }
 714
 715 /**
 716  * hkl_geometry_list_reset: (skip)
 717  * @self:
 718  *
 719  * reset the HklGeometry, in fact it is a sort of clean method remove
 720  * all the items of the list.
 721  **/
 722 void hkl_geometry_list_reset(HklGeometryList *self)
 723 {
 724         HklGeometryListItem *item;
 725         HklGeometryListItem *next;
 726
 727         list_for_each_safe(&self->items, item, next, node)
 728                 hkl_geometry_list_item_free(item);
 729
 730         list_head_init(&self->items);
 731         self->len = 0;
 732 }
 733
 734 /**
 735  * hkl_geometry_list_sort: (skip)
 736  * @self:
 737  * @ref:
 738  *
 739  * sort the #HklGeometryList compare to the distance of the given
 740  * #HklGeometry
 741  **/
 742 void hkl_geometry_list_sort(HklGeometryList *self, HklGeometry *ref)
 743 {
 744         double *distances = alloca(self->len * sizeof(*distances));
 745         size_t *idx = alloca(self->len * sizeof(*idx));
 746         int i = 0;
 747         size_t x;
 748         int j, p;
 749         HklGeometryListItem **items;
 750         HklGeometryListItem *item;
 751         HklGeometryListItem *next;
 752
 753         /* prepare a vector to sort the items */
 754         /* it should be better to sort directly the list instead of
 755          * this hack */
 756         items = malloc(self->len * sizeof(*items));
 757
 758         /* compute the distances once for all */
 759         list_for_each(&self->items, item, node){
 760                 distances[i] = hkl_geometry_distance(ref, item->geometry);
 761                 idx[i] = i;
 762                 items[i] = item;
 763                 i++;
 764         }
 765
 766         /* insertion sorting */
 767         for(i=1; i<self->len; ++i){
 768                 x = idx[i];
 769                 /* find the smallest idx p lower than i with distance[idx[p]] >= distance[x] */
 770                 for(p = 0; distances[idx[p]] < distances[x] && fabs(distances[idx[p]] - distances[x]) > HKL_EPSILON; p++);
 771
 772                 /* move everythings in between p and i */
 773                 for(j=i-1; j>=p; j--)
 774                         idx[j+1] = idx[j];
 775
 776                 idx[p] = x; /* insert the saved idx */
 777         }
 778
 779         list_head_init(&self->items);
 780         for(i=0; i<self->len; ++i)
 781                 list_add_tail(&self->items, &items[idx[i]]->node);
 782         free(items);
 783 }
 784
 785 /**
 786  * hkl_geometry_list_fprintf: (skip)
 787  * @f:
 788  * @self:
 789  *
 790  * print to a file the #HklGeometryList
 791  **/
 792 void hkl_geometry_list_fprintf(FILE *f, const HklGeometryList *self)
 793 {
 794         HklParameter *parameter;
 795         HklGeometry *g;
 796         double value;
 797         size_t axes_len;
 798         size_t i, j;
 799
 800         if(!self)
 801                 return;
 802
 803         fprintf(f, "multiply method: %p \n", self->multiply);
 804         if(self->len){
 805                 HklGeometryListItem *item;
 806
 807                 item = list_top(&self->items, HklGeometryListItem, node);
 808                 axes_len = item->geometry->len;
 809                 fprintf(f, "    ");
 810                 for(i=0; i<axes_len; ++i)
 811                         fprintf(f, "%19s", hkl_axis_get_name(&item->geometry->axes[i]));
 812
 813                 /* geometries */
 814                 list_for_each(&self->items, item, node){
 815                         fprintf(f, "\n%d :", i);
 816                         for(j=0; j<axes_len; ++j) {
 817                                 parameter = (HklParameter *)(&item->geometry->axes[j]);
 818                                 value = hkl_parameter_get_value_unit(parameter);
 819                                 if (parameter->punit)
 820                                         fprintf(f, " % 18.15f %s", value, parameter->punit->repr);
 821                                 else
 822                                         fprintf(f, " % 18.15f", value);
 823
 824                         }
 825                         fprintf(f, "\n   ");
 826                         for(j=0; j<axes_len; ++j) {
 827                                 parameter = (HklParameter *)(&item->geometry->axes[j]);
 828                                 value = gsl_sf_angle_restrict_symm(parameter->value) * hkl_unit_factor(parameter->unit, parameter->punit);
 829                                 if (parameter->punit)
 830                                         fprintf(f, " % 18.15f %s", value, parameter->punit->repr);
 831                                 else
 832                                         fprintf(f, " % 18.15f", value);
 833                         }
 834                         fprintf(f, "\n");
 835                 }
 836         }
 837 }
 838
 839 /**
 840  * hkl_geometry_list_multiply: (skip)
 841  * @self:
 842  *
 843  * apply the multiply lenthod to the #HklGeometry
 844  **/
 845 void hkl_geometry_list_multiply(HklGeometryList *self)
 846 {
 847         HklGeometryListItem *item;
 848         HklGeometryListItem *last;
 849
 850         if(!self || !self->multiply)
 851                 return;
 852
 853         /*
 854          * warning this method change the self->len so we need to save it
 855          * before using the recursive perm_r calls
 856          */
 857         last = list_tail(&self->items, HklGeometryListItem, node);
 858         list_for_each(&self->items, item, node){
 859                 self->multiply(self, item);
 860                 if (item == last)
 861                         break;
 862         }
 863 }
 864
 865 static void perm_r(HklGeometryList *self, HklGeometry *ref, HklGeometry *geometry,
 866                    int perm[], size_t axis_idx)
 867 {
 868         if (axis_idx == geometry->len){
 869                 if(hkl_geometry_distance(ref, geometry) > HKL_EPSILON){
 870                         HklGeometryListItem *item;
 871
 872                         item = hkl_geometry_list_item_new(geometry);
 873                         list_add_tail(&self->items, &item->node);
 874                         self->len++;
 875                 }
 876         }else{
 877                 if(perm[axis_idx]){
 878                         HklAxis *axis;
 879                         double max;
 880                         double value;
 881                         double value0;
 882
 883                         axis = &geometry->axes[axis_idx];
 884                         max = hkl_axis_get_max(axis);
 885                         value = hkl_axis_get_value(axis);
 886                         value0 = value;
 887                         do{
 888                                 /* fprintf(stdout, "\n%d %s, %f", axis_idx, hkl_axis_get_name(axis), value * HKL_RADTODEG); */
 889                                 perm_r(self, ref, geometry, perm, axis_idx + 1);
 890                                 value +=  2*M_PI;
 891                                 if(value <= (max + HKL_EPSILON))
 892                                         hkl_axis_set_value(axis, value);
 893                         }while(value <= (max + HKL_EPSILON));
 894                         hkl_axis_set_value(axis, value0);
 895                 } else
 896                         perm_r(self, ref, geometry, perm, axis_idx + 1);
 897         }
 898 }
 899
 900 void hkl_geometry_list_multiply_from_range(HklGeometryList *self)
 901 {
 902         size_t j;
 903         HklGeometryListItem *item;
 904         HklGeometryListItem *last;
 905
 906         if(!self)
 907                 return;
 908
 909         /*
 910          * warning this method change the self->len so we need to save it
 911          * before using the recursive perm_r calls
 912          */
 913         last = list_tail(&self->items, HklGeometryListItem, node);
 914         list_for_each(&self->items, item, node){
 915                 HklGeometry *geometry;
 916                 HklGeometry *ref;
 917                 int *perm;
 918
 919                 ref = item->geometry;
 920                 geometry = hkl_geometry_new_copy(ref);
 921                 perm = alloca(geometry->len * sizeof(*perm));
 922
 923                 /* find axes to permute and the first solution of thoses axes */
 924                 for(j=0; j<geometry->len; ++j){
 925                         HklAxis *axis = &geometry->axes[j];
 926                         perm[j] = hkl_axis_is_value_compatible_with_range(axis);
 927                         /* fprintf(stdout, "%d %d\n", j, perm[j]); */
 928                         if (perm[j])
 929                                 hkl_axis_set_value_smallest_in_range(axis);
 930                 }
 931                 /*
 932                  * fprintf(stdout, "FIRST SOLUTION\n");
 933                  * hkl_geometry_fprintf(stdout, geometry);
 934                  */
 935
 936                 perm_r(self, ref, geometry, perm, 0);
 937                 hkl_geometry_free(geometry);
 938                 if (item == last)
 939                         break;
 940         }
 941 }
 942
 943 /**
 944  * hkl_geometry_list_remove_invalid: (skip)
 945  * @self:
 946  *
 947  * remove all invalid #HklGeometry from the #HklGeometryList
 948  **/
 949 void hkl_geometry_list_remove_invalid(HklGeometryList *self)
 950 {
 951         HklGeometryListItem *item;
 952         HklGeometryListItem *next;
 953
 954         if(!self)
 955                 return;
 956
 957         list_for_each_safe(&self->items, item, next, node)
 958                 if(!hkl_geometry_is_valid(item->geometry)){
 959                         list_del(&item->node);
 960                         hkl_geometry_list_item_free(item);
 961                         self->len--;
 962                 }
 963 }
 964
 965 /***********************/
 966 /* HklGeometryListItem */
 967 /***********************/
 968
 969 /**
 970  * hkl_geometry_list_item_new: (skip)
 971  * @geometry:
 972  *
 973  * constructor
 974  *
 975  * Returns:
 976  **/
 977 HklGeometryListItem *hkl_geometry_list_item_new(HklGeometry *geometry)
 978 {
 979         HklGeometryListItem *self;
 980
 981         if(!geometry)
 982                 return NULL;
 983
 984         self = HKL_MALLOC(HklGeometryListItem);
 985
 986         self->geometry = hkl_geometry_new_copy(geometry);
 987
 988         return self;
 989 }
 990
 991 /**
 992  * hkl_geometry_list_item_new_copy: (skip)
 993  * @self:
 994  *
 995  * copy constructor
 996  *
 997  * Returns:
 998  **/
 999 HklGeometryListItem *hkl_geometry_list_item_new_copy(const HklGeometryListItem *self)
1000 {
1001         HklGeometryListItem *dup;
1002
1003         if(!self)
1004                 return NULL;
1005
1006         dup = HKL_MALLOC(HklGeometryListItem);
1007
1008         dup->geometry = hkl_geometry_new_copy(self->geometry);
1009
1010         return dup;
1011 }
1012
1013 /**
1014  * hkl_geometry_list_item_free: (skip)
1015  * @self:
1016  *
1017  * destructor
1018  **/
1019 void hkl_geometry_list_item_free(HklGeometryListItem *self)
1020 {
1021         if(!self)
1022                 return;
1023
1024         hkl_geometry_free(self->geometry);
1025         free(self);
1026 }