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