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-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_axes_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_axes_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 \"\"\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: (skip)
 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         size_t i;
 471         int ko = 0;
 472
 473         darray_foreach(axis, self->axes){
 474                 if ((*axis)->changed) {
 475                         ko = 1;
 476                         break;
 477                 }
 478         }
 479
 480         if (ko) {
 481                 HklHolder **holder;
 482
 483                 darray_foreach(holder, self->holders){
 484                         hkl_holder_update(*holder);
 485                 }
 486
 487                 darray_foreach(axis, self->axes){
 488                         (*axis)->changed = FALSE;
 489                 }
 490         }
 491 }
 492
 493 const char *hkl_geometry_name_get(const HklGeometry *self)
 494 {
 495         return hkl_factory_name_get(self->factory);
 496 }
 497
 498 /**
 499  * hkl_geometry_get_axis_idx_by_name: (skip)
 500  * @self:
 501  * @name:
 502  *
 503  * get the index of the axes named @name in the geometry
 504  *
 505  * Returns: -1 if the axis was not found
 506  **/
 507 int hkl_geometry_get_axis_idx_by_name(const HklGeometry *self, const char *name)
 508 {
 509         uint i = 0;
 510         HklParameter **axis;
 511
 512         if (!self || !name)
 513                 return -1;
 514
 515         darray_foreach(axis, self->axes){
 516                 if (!strcmp((*axis)->name, name))
 517                         return i;
 518                 ++i;
 519         }
 520
 521         return -1;
 522 }
 523
 524 /**
 525  * hkl_geometry_get_axis_by_name:
 526  * @self:
 527  * @name:
 528  *
 529  * get an #HklAxis using its name
 530  *
 531  * Returns: (transfer none):
 532  **/
 533 HklParameter *hkl_geometry_get_axis_by_name(HklGeometry *self, const char *name)
 534 {
 535         HklParameter **axis;
 536
 537         darray_foreach(axis, self->axes){
 538                 if (!strcmp((*axis)->name, name))
 539                         return (*axis);
 540         }
 541         return NULL;
 542 }
 543
 544 /**
 545  * hkl_geometry_axes_values_get:
 546  * @self: the this ptr
 547  * @values: (array length=n_values): the values to get
 548  * @n_values: the size of the values array.
 549  * @unit_type: the unit type (default or user) of the returned value
 550  *
 551  * fill the values array with the #HklGeometry axes.
 552  **/
 553 void hkl_geometry_axes_values_get(const HklGeometry *self,
 554                                   double values[], size_t n_values,
 555                                   HklUnitEnum unit_type)
 556 {
 557         size_t i = 0;
 558         HklParameter **axis;
 559
 560         g_return_if_fail (n_values == darray_size(self->axes));
 561
 562         darray_foreach(axis, self->axes){
 563                 values[i++] = hkl_parameter_value_get(*axis, unit_type);
 564         }
 565 }
 566
 567 /**
 568  * hkl_geometry_axes_values_set:
 569  * @self: the this ptr
 570  * @values: (array length=n_values): the values to set.
 571  * @n_values: the length of the values array.
 572  * @unit_type: the unit type (default or user) of the returned value
 573  * @error: return location for a GError, or NULL
 574  *
 575  * Set the #HklGeometry axes values
 576  *
 577  * Returns: TRUE on success, FALSE if an error occurred
 578  **/
 579 int hkl_geometry_axes_values_set(HklGeometry *self,
 580                                  double values[], size_t n_values,
 581                                  HklUnitEnum unit_type,
 582                                  GError **error)
 583 {
 584         uint i = 0;
 585         HklParameter **axis;
 586
 587         hkl_error (error == NULL || *error == NULL);
 588         g_assert(n_values == darray_size(self->axes));
 589
 590         darray_foreach(axis, self->axes){
 591                 if(!hkl_parameter_value_set(*axis, values[i++], unit_type, error)){
 592                         g_assert (error == NULL || *error != NULL);
 593                         return FALSE;
 594                 }
 595         }
 596         g_assert (error == NULL || *error == NULL);
 597
 598         hkl_geometry_update(self);
 599
 600         return TRUE;
 601 }
 602
 603 /**
 604  * hkl_geometry_randomize: (skip)
 605  * @self:
 606  *
 607  * randomize the #HklGeometry
 608  **/
 609 void hkl_geometry_randomize(HklGeometry *self)
 610 {
 611         HklParameter **axis;
 612
 613         darray_foreach(axis, self->axes){
 614                 hkl_parameter_randomize(*axis);
 615         }
 616         hkl_geometry_update(self);
 617 }
 618
 619 /**
 620  * hkl_geometry_set_values_v: (skip)
 621  * @self:
 622  * @unit_type: the unit type (default or user) of the returned value
 623  * @error:
 624  * "...:
 625  *
 626  * set the axes values
 627  *
 628  * Returns:
 629  **/
 630 int hkl_geometry_set_values_v(HklGeometry *self, HklUnitEnum unit_type, GError **error, ...)
 631 {
 632         va_list ap;
 633         HklParameter **axis;
 634
 635         hkl_error (error == NULL || *error == NULL);
 636
 637         va_start(ap, error);
 638         darray_foreach(axis, self->axes){
 639                 if(!hkl_parameter_value_set(*axis,
 640                                             va_arg(ap, double),
 641                                             unit_type, error)){
 642                         g_assert (error == NULL || *error != NULL);
 643                         va_end(ap);
 644                         hkl_geometry_update(self);
 645                         return FALSE;
 646                 }
 647         }
 648         g_assert (error == NULL || *error == NULL);
 649
 650         va_end(ap);
 651
 652         hkl_geometry_update(self);
 653
 654         return TRUE;
 655 }
 656
 657 /**
 658  * hkl_geometry_distance:
 659  * @self: the this ptr
 660  * @ref: the #HklGeometry to compare with
 661  *
 662  * compute the distance between two #HklGeometries
 663  *
 664  * Returns: the distance between the two geometries
 665  **/
 666 double hkl_geometry_distance(const HklGeometry *self,
 667                              const HklGeometry *ref)
 668 {
 669         size_t i;
 670         double value1, value2;
 671         double distance = 0.;
 672
 673         if (!self || !ref)
 674                 return 0.;
 675
 676         for(i=0; i<darray_size(self->axes); ++i){
 677                 value1 = darray_item(self->axes, i)->_value;
 678                 value2 = darray_item(ref->axes, i)->_value;
 679                 distance += fabs(value2 - value1);
 680         }
 681
 682         return distance;
 683 }
 684
 685 /**
 686  * hkl_geometry_distance_orthodromic: (skip)
 687  * @self: the this ptr
 688  * @ref: the reference #HklGeometry to compare with.
 689  *
 690  * Returns: the orthodromique distance
 691  **/
 692 double hkl_geometry_distance_orthodromic(const HklGeometry *self,
 693                                          const HklGeometry *ref)
 694 {
 695         size_t i;
 696         double value1, value2;
 697         double distance = 0.;
 698
 699         if (!self || !ref)
 700                 return 0.;
 701
 702         for(i=0; i<darray_size(self->axes); ++i){
 703                 double d;
 704
 705                 value1 = darray_item(self->axes, i)->_value;
 706                 value2 = darray_item(ref->axes, i)->_value;
 707                 d = fabs(gsl_sf_angle_restrict_symm(value2) - gsl_sf_angle_restrict_symm(value1));
 708                 /* as M_PI and -M_PI are included in the GSL restriction */
 709                 if (d > M_PI)
 710                         d = 2*M_PI - d;
 711                 distance += d;
 712         }
 713
 714         return distance;
 715 }
 716
 717 /**
 718  * hkl_geometry_is_valid: (skip)
 719  * @self:
 720  *
 721  * check if all axes of the #HklGeometry are valid.
 722  *
 723  * Returns:
 724  **/
 725 int hkl_geometry_is_valid(const HklGeometry *self)
 726 {
 727         HklParameter **axis;
 728
 729         darray_foreach(axis, self->axes){
 730                 if(!hkl_parameter_is_valid(*axis))
 731                         return FALSE;
 732         }
 733
 734         return TRUE;
 735 }
 736
 737 /**
 738  * hkl_geometry_closest_from_geometry_with_range: (skip)
 739  * @self:
 740  * @ref:
 741  *
 742  * get the closest axes values in the HklInterval compatible with the
 743  * current axes values
 744  *
 745  * Returns:
 746  **/
 747 int hkl_geometry_closest_from_geometry_with_range(HklGeometry *self,
 748                                                   const HklGeometry *ref)
 749 {
 750         size_t i;
 751         uint len = darray_size(self->axes);
 752         double *values = alloca(len * sizeof(*values));
 753         int ko = FALSE;
 754
 755         for(i=0;i<len;++i){
 756                 values[i] = hkl_parameter_value_get_closest(darray_item(self->axes, i),
 757                                                             darray_item(ref->axes, i));
 758                 if(gsl_isnan(values[i])){
 759                         ko = TRUE;
 760                         break;
 761                 }
 762         }
 763         if(!ko){
 764                 for(i=0;i<len;++i)
 765                         hkl_parameter_value_set(darray_item(self->axes, i),
 766                                                 values[i],
 767                                                 HKL_UNIT_DEFAULT, NULL);
 768                 hkl_geometry_update(self);
 769         }
 770         return ko;
 771 }
 772
 773 /**
 774  * hkl_geometry_fprintf: (skip)
 775  * @file:
 776  * @self:
 777  *
 778  * print into a file the #HklGeometry
 779  **/
 780 void hkl_geometry_fprintf(FILE *file, const HklGeometry *self)
 781 {
 782         fprintf(file, " HklGeometry type: \"%s\" wavelength: %f",
 783                 self->factory->name,
 784                 self->source.wave_length);
 785         for(unsigned int i=0; i<darray_size(self->axes); ++i){
 786                 fprintf(file, " ");
 787                 hkl_parameter_fprintf(file, darray_item(self->axes, i));
 788         }
 789 }
 790
 791 /*******************/
 792 /* HklGeometryList */
 793 /*******************/
 794
 795 /**
 796  * hkl_geometry_list_new: (skip)
 797  *
 798  * constructor
 799  *
 800  * Returns:
 801  **/
 802 HklGeometryList *hkl_geometry_list_new(void)
 803 {
 804         HklGeometryList *self;
 805
 806         self = HKL_MALLOC(HklGeometryList);
 807
 808         list_head_init(&self->items);
 809         self->n_items = 0;
 810         self->multiply = NULL;
 811
 812         return self;
 813 }
 814
 815 /**
 816  * hkl_geometry_list_new_copy: (skip)
 817  * @self:
 818  *
 819  * copy constructor
 820  *
 821  * Returns:
 822  **/
 823 HklGeometryList *hkl_geometry_list_new_copy(const HklGeometryList *self)
 824 {
 825         HklGeometryList *dup;
 826         HklGeometryListItem *item;
 827
 828         if (!self)
 829                 return NULL;
 830
 831         dup = HKL_MALLOC(HklGeometryList);
 832
 833         list_head_init(&dup->items);
 834         /* now copy the item arrays */
 835         list_for_each(&self->items, item, list){
 836                 list_add_tail(&dup->items,
 837                               &hkl_geometry_list_item_new_copy(item)->list);
 838         }
 839         dup->n_items = self->n_items;
 840         dup->multiply = self->multiply;
 841
 842         return dup;
 843 }
 844
 845 /**
 846  * hkl_geometry_list_free: (skip)
 847  * @self:
 848  *
 849  * destructor
 850  **/
 851 void hkl_geometry_list_free(HklGeometryList *self)
 852 {
 853         hkl_geometry_list_reset(self);
 854         free(self);
 855 }
 856
 857 /**
 858  * hkl_geometry_list_add: (skip)
 859  * @self: The current #HklGeometryList
 860  * @geometry: the #HklGeometry to add
 861  *
 862  * this method Add a geometry to the geometries
 863  *
 864  * This method try to be clever by allocating memory only if the
 865  * current length of the geometries is not large enought. Then it just
 866  * set the geometry axes and copy it to the right geometries. We do
 867  * not gives the x len as it is equal to the self->axes_len.
 868  **/
 869 void hkl_geometry_list_add(HklGeometryList *self, HklGeometry *geometry)
 870 {
 871         HklGeometryListItem *item;
 872
 873         /* now check if the geometry is already in the geometry list */
 874         list_for_each(&self->items, item, list){
 875                 if (hkl_geometry_distance_orthodromic(geometry,
 876                                                       item->geometry) < HKL_EPSILON)
 877                         return;
 878         }
 879
 880         list_add_tail(&self->items,
 881                       &hkl_geometry_list_item_new(geometry)->list);
 882         self->n_items += 1;
 883 }
 884
 885 /**
 886  * hkl_geometry_list_n_items_get: (skip)
 887  * @self: the this ptr
 888  *
 889  * get the number of items in the #HklGeometryList
 890  *
 891  * Returns: the number of items in the list
 892  **/
 893 size_t hkl_geometry_list_n_items_get(const HklGeometryList *self)
 894 {
 895         return self->n_items;
 896 }
 897
 898 /**
 899  * hkl_geometry_list_items_first_get: (skip)
 900  * @self: the this ptr
 901  *
 902  * get the first solution of the #HklGeometryList
 903  *
 904  * Returns: the first solution of the list
 905  **/
 906 const HklGeometryListItem *hkl_geometry_list_items_first_get(const HklGeometryList *self)
 907 {
 908         return list_top(&self->items, HklGeometryListItem, list);
 909 }
 910
 911 /**
 912  * hkl_geometry_list_items_next_get: (skip)
 913  * @self: the this ptr
 914  * @item: the current #HklGeometryListItem solution of the #HklGeometryList
 915  *
 916  * get the next solution of the #HklGeometryList from the current item location.
 917  *
 918  * Returns: the next solution of the list
 919  **/
 920 const HklGeometryListItem *hkl_geometry_list_items_next_get(const HklGeometryList *self,
 921                                                             const HklGeometryListItem *item)
 922 {
 923         return list_next(&self->items, item, list);
 924 }
 925
 926 /**
 927  * hkl_geometry_list_reset: (skip)
 928  * @self: the this ptr
 929  *
 930  * reset the HklGeometry, in fact it is a sort of clean method remove
 931  * all the items of the list.
 932  **/
 933 void hkl_geometry_list_reset(HklGeometryList *self)
 934 {
 935         HklGeometryListItem *item;
 936         HklGeometryListItem *next;
 937
 938         list_for_each_safe(&self->items, item, next, list)
 939                 hkl_geometry_list_item_free(item);
 940
 941         list_head_init(&self->items);
 942         self->n_items = 0;
 943 }
 944
 945 /**
 946  * hkl_geometry_list_sort: (skip)
 947  * @self:
 948  * @ref:
 949  *
 950  * sort the #HklGeometryList compare to the distance of the given
 951  * #HklGeometry
 952  **/
 953 void hkl_geometry_list_sort(HklGeometryList *self, HklGeometry *ref)
 954 {
 955         double *distances = alloca(self->n_items * sizeof(*distances));
 956         size_t *idx = alloca(self->n_items * sizeof(*idx));
 957         HklGeometryListItem **items = alloca(self->n_items * sizeof(*items));
 958         HklGeometryListItem *item;
 959         HklGeometryListItem *next;
 960         int i = 0;
 961         size_t x;
 962         int j, p;
 963
 964         /* compute the distances once for all */
 965         list_for_each(&self->items, item, list){
 966                 distances[i] = hkl_geometry_distance(ref, item->geometry);
 967                 idx[i] = i;
 968                 items[i] = item;
 969                 i++;
 970         }
 971
 972         /* insertion sorting */
 973         for(i=1; i<self->n_items; ++i){
 974                 x = idx[i];
 975                 /* find the smallest idx p lower than i with distance[idx[p]] >= distance[x] */
 976                 for(p = 0; distances[idx[p]] < distances[x] && fabs(distances[idx[p]] - distances[x]) > HKL_EPSILON; p++);
 977
 978                 /* move everythings in between p and i */
 979                 for(j=i-1; j>=p; j--)
 980                         idx[j+1] = idx[j];
 981
 982                 idx[p] = x; /* insert the saved idx */
 983         }
 984
 985         list_head_init(&self->items);
 986
 987         for(i=0; i<self->n_items; ++i)
 988                 list_add_tail(&self->items, &items[idx[i]]->list);
 989 }
 990
 991 /**
 992  * hkl_geometry_list_fprintf: (skip)
 993  * @f:
 994  * @self:
 995  *
 996  * print to a file the #HklGeometryList
 997  **/
 998 void hkl_geometry_list_fprintf(FILE *f, const HklGeometryList *self)
 999 {
1000         HklGeometry *g;
1001         uint i = 0;
1002         double value;
1003
1004         if(!self)
1005                 return;
1006
1007         fprintf(f, "multiply method: %p \n", self->multiply);
1008         if(self->n_items){
1009                 HklGeometryListItem *item;
1010                 HklParameter **axis;
1011
1012                 fprintf(f, "    ");
1013                 darray_foreach(axis, list_top(&self->items, HklGeometryListItem, list)->geometry->axes){
1014                         fprintf(f, "%19s", (*axis)->name);
1015                 }
1016
1017                 /* geometries */
1018                 list_for_each(&self->items, item, list){
1019                         fprintf(f, "\n%d :", i++);
1020                         darray_foreach(axis, item->geometry->axes){
1021                                 value = hkl_parameter_value_get(*axis, HKL_UNIT_DEFAULT);
1022                                 if ((*axis)->punit)
1023                                         fprintf(f, " % 18.15f %s", value, (*axis)->punit->repr);
1024                                 else
1025                                         fprintf(f, " % 18.15f", value);
1026
1027                         }
1028                         fprintf(f, "\n   ");
1029                         darray_foreach(axis, item->geometry->axes){
1030                                 value = hkl_parameter_value_get(*axis, HKL_UNIT_DEFAULT);
1031                                 value = gsl_sf_angle_restrict_symm(value);
1032                                 value *= hkl_unit_factor((*axis)->unit,
1033                                                          (*axis)->punit);
1034                                 if ((*axis)->punit)
1035                                         fprintf(f, " % 18.15f %s", value, (*axis)->punit->repr);
1036                                 else
1037                                         fprintf(f, " % 18.15f", value);
1038                         }
1039                         fprintf(f, "\n");
1040                 }
1041         }
1042 }
1043
1044 /**
1045  * hkl_geometry_list_multiply: (skip)
1046  * @self:
1047  *
1048  * apply the multiply lenthod to the #HklGeometry
1049  **/
1050 void hkl_geometry_list_multiply(HklGeometryList *self)
1051 {
1052         uint i = 0;
1053         uint len = self->n_items;
1054         HklGeometryListItem *item;
1055
1056         if(!self || !self->multiply)
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         for(i=0, item=list_top(&self->items, HklGeometryListItem, list);
1064             i<len;
1065             ++i, item=list_next(&self->items, item, list))
1066                 self->multiply(self, item);
1067 }
1068
1069 static void perm_r(HklGeometryList *self, const HklGeometry *ref,
1070                    const HklGeometry *geometry, const int perm[],
1071                    const unsigned int axis_idx)
1072 {
1073         if (axis_idx == darray_size(geometry->axes)){
1074                 if(hkl_geometry_distance(geometry, ref) > HKL_EPSILON){
1075                         list_add_tail(&self->items,
1076                                       &hkl_geometry_list_item_new(geometry)->list);
1077                         self->n_items++;
1078                 }
1079         }else{
1080                 if(perm[axis_idx]){
1081                         HklParameter *axis = darray_item(geometry->axes, axis_idx);
1082                         const double max = axis->range.max;;
1083                         const double value0 = axis->_value;
1084                         double value;
1085
1086                         value = value0;
1087                         do{
1088                                 /* fprintf(stdout, "\n%d %s, %f", axis_idx, hkl_axis_get_name(axis), value * HKL_RADTODEG); */
1089                                 perm_r(self, ref, geometry, perm, axis_idx + 1);
1090                                 value +=  2*M_PI;
1091                                 if(value <= (max + HKL_EPSILON)){
1092                                         /* optimisation here: */
1093                                         /* instead of using set_value
1094                                          * we directly write the
1095                                          * HklParameter value, BEWARE
1096                                          * that it require that
1097                                          * HklParameter is a rotation
1098                                          * (for now it is always
1099                                          * true */
1100                                         axis->_value = value;
1101                                 }
1102                         }while(value <= (max + HKL_EPSILON));
1103                         /* restore the initial value */
1104                         axis->_value = value0;
1105                 } else
1106                         perm_r(self, ref, geometry, perm, axis_idx + 1);
1107         }
1108 }
1109
1110 void hkl_geometry_list_multiply_from_range(HklGeometryList *self)
1111 {
1112         uint i;
1113         uint len = self->n_items;
1114         size_t j = 0;
1115         const HklGeometryListItem *item;
1116
1117         if(!self)
1118                 return;
1119
1120         /*
1121          * warning this method change the self->len so we need to save it
1122          * before using the recursive perm_r calls
1123          */
1124
1125         for(i=0, item=list_top(&self->items, HklGeometryListItem, list);
1126             i<len;
1127             ++i, item=list_next(&self->items, item, list)){
1128                 HklGeometry *geometry;
1129                 HklParameter **axis;
1130                 int *perm;
1131
1132                 geometry = hkl_geometry_new_copy(item->geometry);
1133                 perm = alloca(darray_size(geometry->axes) * sizeof(*perm));
1134
1135                 /* find axes to permute and the first solution of thoses axes */
1136                 darray_foreach(axis, geometry->axes){
1137                         perm[j] = hkl_parameter_is_valid(*axis);
1138                         /* fprintf(stdout, "%d %d\n", j, perm[j]); */
1139                         if (perm[j])
1140                                 hkl_parameter_value_set_smallest_in_range(*axis);
1141                         ++j;
1142                 }
1143                 /*
1144                  * fprintf(stdout, "FIRST SOLUTION\n");
1145                  * hkl_geometry_fprintf(stdout, geometry);
1146                  */
1147
1148                 perm_r(self, item->geometry, geometry, perm, 0);
1149                 hkl_geometry_free(geometry);
1150         }
1151 }
1152
1153 /**
1154  * hkl_geometry_list_remove_invalid: (skip)
1155  * @self:
1156  *
1157  * remove all invalid #HklGeometry from the #HklGeometryList
1158  **/
1159 void hkl_geometry_list_remove_invalid(HklGeometryList *self)
1160 {
1161         HklGeometryListItem *item, *next;
1162
1163         list_for_each_safe(&self->items, item, next, list)
1164                 if(!hkl_geometry_is_valid(item->geometry)){
1165                         list_del(&item->list);
1166                         self->n_items--;
1167                         hkl_geometry_list_item_free(item);
1168                 }
1169 }
1170
1171 /***********************/
1172 /* HklGeometryListItem */
1173 /***********************/
1174
1175 /**
1176  * hkl_geometry_list_item_new: (skip)
1177  * @geometry:
1178  *
1179  * constructor
1180  *
1181  * Returns:
1182  **/
1183 HklGeometryListItem *hkl_geometry_list_item_new(const HklGeometry *geometry)
1184 {
1185         HklGeometryListItem *self;
1186
1187         if(!geometry)
1188                 return NULL;
1189
1190         self = HKL_MALLOC(HklGeometryListItem);
1191
1192         self->geometry = hkl_geometry_new_copy(geometry);
1193
1194         return self;
1195 }
1196
1197 /**
1198  * hkl_geometry_list_item_new_copy: (skip)
1199  * @self:
1200  *
1201  * copy constructor
1202  *
1203  * Returns:
1204  **/
1205 HklGeometryListItem *hkl_geometry_list_item_new_copy(const HklGeometryListItem *self)
1206 {
1207         HklGeometryListItem *dup;
1208
1209         if(!self)
1210                 return NULL;
1211
1212         dup = HKL_MALLOC(HklGeometryListItem);
1213
1214         dup->geometry = hkl_geometry_new_copy(self->geometry);
1215
1216         return dup;
1217 }
1218
1219 /**
1220  * hkl_geometry_list_item_free: (skip)
1221  * @self:
1222  *
1223  * destructor
1224  **/
1225 void hkl_geometry_list_item_free(HklGeometryListItem *self)
1226 {
1227         if(!self)
1228                 return;
1229
1230         hkl_geometry_free(self->geometry);
1231         free(self);
1232 }
1233
1234 /**
1235  * hkl_geometry_list_item_geometry_get:
1236  * @self: the this ptr
1237  *
1238  * Return value: The geometry contain inside the HklGeometryListItem
1239  **/
1240 const HklGeometry *hkl_geometry_list_item_geometry_get(const HklGeometryListItem *self)
1241 {
1242         return self->geometry;
1243 }