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