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