hkl/hkl-lattice.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 <math.h>                       // for cos, sin, M_PI, atan2, sqrt
  23 #include <stdio.h>                      // for fprintf, FILE
  24 #include <stdlib.h>                     // for NULL, free
  25 #include "hkl-lattice-private.h"        // for _HklLattice
  26 #include "hkl-macros-private.h"         // for HKL_MALLOC
  27 #include "hkl-matrix-private.h"         // for _HklMatrix
  28 #include "hkl-parameter-private.h"      // for hkl_parameter_init_copy, etc
  29 #include "hkl-unit-private.h"           // for hkl_unit_length_nm, etc
  30 #include "hkl-vector-private.h"         // for hkl_vector_angle, etc
  31 #include "hkl.h"                        // for HklLattice, etc
  32
  33 /* private */
  34
  35 static double convert_to_default(const HklParameter *p, double value, HklUnitEnum unit_type)
  36 {
  37         switch(unit_type){
  38         case HKL_UNIT_DEFAULT:
  39                 return value;
  40         case HKL_UNIT_USER:
  41                 return value / hkl_unit_factor(p->unit, p->punit);
  42         default:
  43                 return NAN;
  44         }
  45 }
  46
  47 static int check_lattice_param(double a, double b, double c,
  48                                double alpha, double beta, double gamma,
  49                                double *volume,
  50                                GError **error)
  51 {
  52         hkl_error (error == NULL || *error == NULL);
  53
  54         double D = 1. - cos(alpha)*cos(alpha) - cos(beta)*cos(beta)
  55                 - cos(gamma)*cos(gamma) + 2. * cos(alpha)*cos(beta)*cos(gamma);
  56
  57         if (D < 0.){
  58                 g_set_error(error,
  59                             HKL_LATTICE_ERROR,
  60                             HKL_LATTICE_CHECK_LATTICE,
  61                             "these lattice parameters are not valid, check alpha, beta and gamma");
  62                 return FALSE;
  63         }else{
  64                 *volume = a * b * c * sqrt(D);
  65                 return TRUE;
  66         }
  67 }
  68
  69 /* public */
  70
  71 /**
  72  * hkl_lattice_new:
  73  * @a: the length of the a parameter
  74  * @b: the length of the b parameter
  75  * @c: the length of the c parameter
  76  * @alpha: the angle between b and c (radian)
  77  * @beta: the angle between a and c (radian)
  78  * @gamma: the angle between a and b (radian)
  79  * @error: return location for a GError, or NULL
  80  *
  81  * constructor
  82  *
  83  * Returns: a new HklLattice
  84  **/
  85 HklLattice *hkl_lattice_new(double a, double b, double c,
  86                             double alpha, double beta, double gamma,
  87                             GError **error)
  88 {
  89         HklLattice *self = NULL;
  90         double volume;
  91
  92         hkl_error (error == NULL || *error == NULL);
  93
  94         if(!check_lattice_param(a, b, c, alpha, beta, gamma, &volume, error))
  95         {
  96                 g_assert (error == NULL || *error != NULL);
  97                 return FALSE;
  98         }
  99         g_assert (error == NULL || *error == NULL);
 100
 101         self = HKL_MALLOC(HklLattice);
 102
 103         self->a = hkl_parameter_new("a", "The length of the first lattice vector",
 104                                     0, a, a+10,
 105                                     TRUE, TRUE,
 106                                     &hkl_unit_length_nm,
 107                                     &hkl_unit_length_nm);
 108         self->b = hkl_parameter_new("b", "The length of the second lattice vector",
 109                                     0, b, b+10,
 110                                     TRUE, TRUE,
 111                                     &hkl_unit_length_nm,
 112                                     &hkl_unit_length_nm);
 113         self->c = hkl_parameter_new("c", "The length of the third lattice vector",
 114                                     0, c, c+10,
 115                                     TRUE, TRUE,
 116                                     &hkl_unit_length_nm,
 117                                     &hkl_unit_length_nm);
 118         self->alpha = hkl_parameter_new("alpha",
 119                                         "The angle between the second and third lattice vector",
 120                                         -M_PI, alpha, M_PI,
 121                                         TRUE, TRUE,
 122                                         &hkl_unit_angle_rad,
 123                                         &hkl_unit_angle_deg);
 124         self->beta = hkl_parameter_new("beta",
 125                                        "The angle between the first and third lattice vector",
 126                                        -M_PI, beta, M_PI,
 127                                        TRUE, TRUE,
 128                                        &hkl_unit_angle_rad,
 129                                        &hkl_unit_angle_deg);
 130         self->gamma = hkl_parameter_new("gamma",
 131                                         "The angle between the first and second lattice vector",
 132                                         -M_PI, gamma, M_PI,
 133                                         TRUE, TRUE,
 134                                         &hkl_unit_angle_rad,
 135                                         &hkl_unit_angle_deg);
 136         self->volume = hkl_parameter_new("volume",
 137                                          "The volume of the lattice",
 138                                          0, volume, a*b*c,
 139                                          FALSE, FALSE,
 140                                          &hkl_unit_length_nm,
 141                                          &hkl_unit_length_nm);
 142
 143         return self;
 144 }
 145
 146 /**
 147  * hkl_lattice_new_copy: (skip)
 148  * @self:
 149  *
 150  * copy constructor
 151  *
 152  * Returns:
 153  **/
 154 HklLattice *hkl_lattice_new_copy(const HklLattice *self)
 155 {
 156         HklLattice *copy = NULL;
 157
 158         copy = HKL_MALLOC(HklLattice);
 159
 160         copy->a = hkl_parameter_new_copy(self->a);
 161         copy->b = hkl_parameter_new_copy(self->b);
 162         copy->c = hkl_parameter_new_copy(self->c);
 163         copy->alpha = hkl_parameter_new_copy(self->alpha);
 164         copy->beta = hkl_parameter_new_copy(self->beta);
 165         copy->gamma = hkl_parameter_new_copy(self->gamma);
 166         copy->volume = hkl_parameter_new_copy(self->volume);
 167
 168         return copy;
 169 }
 170
 171 /**
 172  * hkl_lattice_new_default: (skip)
 173  *
 174  * default constructor
 175  *
 176  * Returns:
 177  **/
 178 HklLattice* hkl_lattice_new_default(void)
 179 {
 180         return hkl_lattice_new(1.54, 1.54, 1.54,
 181                                90*HKL_DEGTORAD, 90*HKL_DEGTORAD, 90*HKL_DEGTORAD,
 182                                NULL);
 183 }
 184
 185 /**
 186  * hkl_lattice_free: (skip)
 187  * @self:
 188  *
 189  * destructor
 190  **/
 191 void hkl_lattice_free(HklLattice *self)
 192 {
 193         hkl_parameter_free(self->a);
 194         hkl_parameter_free(self->b);
 195         hkl_parameter_free(self->c);
 196         hkl_parameter_free(self->alpha);
 197         hkl_parameter_free(self->beta);
 198         hkl_parameter_free(self->gamma);
 199         hkl_parameter_free(self->volume);
 200         free(self);
 201 }
 202
 203 #define HKL_LATTICE_X_SET(_p, _parameter, _error) do{                   \
 204                 hkl_error ((_error) == NULL || *(_error) == NULL);      \
 205                 double a, b, c, alpha, beta, gamma;                     \
 206                 /* check if the combinaison of parameters is ok */      \
 207                 hkl_lattice_get(self, &a, &b, &c, &alpha, &beta, &gamma, HKL_UNIT_DEFAULT); \
 208                 _p = hkl_parameter_value_get((_parameter), HKL_UNIT_DEFAULT); \
 209                                                                         \
 210                 if(!hkl_lattice_set(self, a, b, c, alpha, beta, gamma, HKL_UNIT_DEFAULT, (_error))){ \
 211                         g_assert ((_error) == NULL || *(_error) != NULL); \
 212                         return FALSE;                                   \
 213                 }                                                       \
 214                 g_assert ((_error) == NULL || *(_error) == NULL);       \
 215                 return hkl_parameter_init_copy(self->_p, (_parameter), (_error)); \
 216         }while(0)
 217
 218 /**
 219  * hkl_lattice_a_get:
 220  * @self: the this ptr
 221  **/
 222 const HklParameter *hkl_lattice_a_get(const HklLattice *self)
 223 {
 224         return self->a;
 225 }
 226
 227 /**
 228  * hkl_lattice_a_set:
 229  * @self: the this ptr
 230  * @parameter: the parameter to set
 231  * @error: return location for a GError, or NULL
 232  *
 233  * Returns: TRUE on success, FALSE if an error occurred
 234  **/
 235 int hkl_lattice_a_set(HklLattice *self, const HklParameter *parameter,
 236                       GError **error)
 237 {
 238         HKL_LATTICE_X_SET(a, parameter, error);
 239 }
 240
 241 /**
 242  * hkl_lattice_b_get:
 243  * @self: the this ptr
 244  **/
 245 const HklParameter *hkl_lattice_b_get(const HklLattice *self)
 246 {
 247         return self->b;
 248 }
 249
 250 /**
 251  * hkl_lattice_b_set:
 252  * @self: the this ptr
 253  * @parameter: the parameter to set
 254  * @error: return location for a GError, or NULL
 255  *
 256  * Returns: TRUE on success, FALSE if an error occurred
 257  **/
 258 int hkl_lattice_b_set(HklLattice *self, const HklParameter *parameter,
 259                       GError **error)
 260 {
 261         HKL_LATTICE_X_SET(b, parameter, error);
 262 }
 263
 264 /**
 265  * hkl_lattice_c_get:
 266  * @self: the this ptr
 267  **/
 268 const HklParameter *hkl_lattice_c_get(const HklLattice *self)
 269 {
 270         return self->c;
 271 }
 272
 273 /**
 274  * hkl_lattice_c_set:
 275  * @self: the this ptr
 276  * @parameter: the parameter to set
 277  * @error: return location for a GError, or NULL
 278  *
 279  * Returns: TRUE on success, FALSE if an error occurred
 280  **/
 281 int hkl_lattice_c_set(HklLattice *self, const HklParameter *parameter,
 282                       GError **error)
 283 {
 284         HKL_LATTICE_X_SET(c, parameter, error);
 285 }
 286
 287 /**
 288  * hkl_lattice_alpha_get:
 289  * @self: the this ptr
 290  **/
 291 const HklParameter *hkl_lattice_alpha_get(const HklLattice *self)
 292 {
 293         return self->alpha;
 294 }
 295
 296 /**
 297  * hkl_lattice_alpha_set:
 298  * @self: the this ptr
 299  * @parameter: the parameter to set
 300  * @error: return location for a GError, or NULL
 301  *
 302  * Returns: TRUE on success, FALSE if an error occurred
 303  **/
 304 int hkl_lattice_alpha_set(HklLattice *self, const HklParameter *parameter,
 305                           GError **error)
 306 {
 307         HKL_LATTICE_X_SET(alpha, parameter, error);
 308 }
 309
 310 /**
 311  * hkl_lattice_beta_get:
 312  * @self: the this ptr
 313  **/
 314 const HklParameter *hkl_lattice_beta_get(const HklLattice *self)
 315 {
 316         return self->beta;
 317 }
 318
 319 /**
 320  * hkl_lattice_beta_set:
 321  * @self: the this ptr
 322  * @parameter: the parameter to set
 323  * @error: return location for a GError, or NULL
 324  *
 325  * Returns: TRUE on success, FALSE if an error occurred
 326  **/
 327 int hkl_lattice_beta_set(HklLattice *self, const HklParameter *parameter,
 328                          GError **error)
 329 {
 330         HKL_LATTICE_X_SET(beta, parameter, error);
 331 }
 332
 333 /**
 334  * hkl_lattice_gamma_get:
 335  * @self: the this ptr
 336  **/
 337 const HklParameter *hkl_lattice_gamma_get(const HklLattice *self)
 338 {
 339         return self->gamma;
 340 }
 341
 342 /**
 343  * hkl_lattice_gamma_set:
 344  * @self: the this ptr
 345  * @parameter: the parameter to set
 346  * @error: return location for a GError, or NULL
 347  *
 348  * Returns: TRUE on success, FALSE if an error occurred
 349  **/
 350 int hkl_lattice_gamma_set(HklLattice *self, const HklParameter *parameter,
 351                           GError **error)
 352 {
 353         HKL_LATTICE_X_SET(gamma, parameter, error);
 354 }
 355
 356 /**
 357  * hkl_lattice_volume_get:
 358  * @self: the this ptr
 359  **/
 360 const HklParameter *hkl_lattice_volume_get(const HklLattice *self)
 361 {
 362         return self->volume;
 363 }
 364
 365 /**
 366  * hkl_lattice_lattice_set: (skip)
 367  * @self: the this ptr
 368  * @lattice: the lattice to set from.
 369  **/
 370 void hkl_lattice_lattice_set(HklLattice *self, const HklLattice *lattice)
 371 {
 372         if (self == lattice)
 373                 return;
 374
 375         hkl_parameter_init_copy(self->a, lattice->a, NULL);
 376         hkl_parameter_init_copy(self->b, lattice->b, NULL);
 377         hkl_parameter_init_copy(self->c, lattice->c, NULL);
 378         hkl_parameter_init_copy(self->alpha, lattice->alpha, NULL);
 379         hkl_parameter_init_copy(self->beta, lattice->beta, NULL);
 380         hkl_parameter_init_copy(self->gamma, lattice->gamma, NULL);
 381         hkl_parameter_init_copy(self->volume, lattice->volume, NULL);
 382 }
 383
 384 /**
 385  * hkl_lattice_set:
 386  * @self:
 387  * @a:
 388  * @b:
 389  * @c:
 390  * @alpha:
 391  * @beta:
 392  * @gamma:
 393  *
 394  * set the lattice parameters
 395  *
 396  * Returns:
 397  **/
 398 int hkl_lattice_set(HklLattice *self,
 399                     double a, double b, double c,
 400                     double alpha, double beta, double gamma,
 401                     HklUnitEnum unit_type, GError **error)
 402 {
 403         hkl_error (error == NULL || *error == NULL);
 404
 405         double _a, _b, _c, _alpha, _beta, _gamma;
 406         double _volume;
 407
 408         _a = convert_to_default(self->a, a, unit_type);
 409         _b = convert_to_default(self->b, b, unit_type);
 410         _c = convert_to_default(self->c, c, unit_type);
 411         _alpha = convert_to_default(self->alpha, alpha, unit_type);
 412         _beta = convert_to_default(self->beta, beta, unit_type);
 413         _gamma = convert_to_default(self->gamma, gamma, unit_type);
 414
 415         /* need to do the conversion before the check */
 416         if(!check_lattice_param(_a, _b, _c, _alpha, _beta, _gamma, &_volume, error)){
 417                 g_assert (error == NULL || *error != NULL);
 418                 return FALSE;
 419         }
 420         g_assert (error == NULL || *error == NULL);
 421
 422         hkl_parameter_value_set(self->a, _a, HKL_UNIT_DEFAULT, NULL);
 423         hkl_parameter_value_set(self->b, _b, HKL_UNIT_DEFAULT, NULL);
 424         hkl_parameter_value_set(self->c, _c, HKL_UNIT_DEFAULT, NULL);
 425         hkl_parameter_value_set(self->alpha, _alpha, HKL_UNIT_DEFAULT, NULL);
 426         hkl_parameter_value_set(self->beta, _beta, HKL_UNIT_DEFAULT, NULL);
 427         hkl_parameter_value_set(self->gamma, _gamma, HKL_UNIT_DEFAULT, NULL);
 428
 429         hkl_parameter_value_set(self->volume, _volume, HKL_UNIT_DEFAULT, NULL);
 430
 431         return TRUE;
 432 }
 433
 434 /**
 435  * hkl_lattice_get:
 436  * @self:
 437  * @a: (out caller-allocates):
 438  * @b: (out caller-allocates):
 439  * @c: (out caller-allocates):
 440  * @alpha: (out caller-allocates):
 441  * @beta: (out caller-allocates):
 442  * @gamma: (out caller-allocates):
 443  *
 444  * get the lattice parameters
 445  * Return value: all the parameters
 446  **/
 447 void hkl_lattice_get(const HklLattice *self,
 448                      double *a, double *b, double *c,
 449                      double *alpha, double *beta, double *gamma,
 450                      HklUnitEnum unit_type)
 451 {
 452         *a = hkl_parameter_value_get(self->a, unit_type);
 453         *b = hkl_parameter_value_get(self->b, unit_type);
 454         *c = hkl_parameter_value_get(self->c, unit_type);
 455         *alpha = hkl_parameter_value_get(self->alpha, unit_type);
 456         *beta = hkl_parameter_value_get(self->beta, unit_type);
 457         *gamma = hkl_parameter_value_get(self->gamma, unit_type);
 458 }
 459
 460 /**
 461  * hkl_lattice_get_B: (skip)
 462  * @self:
 463  * @B: (out): where to store the B matrix
 464  *
 465  * Get the B matrix from the lattice parameters
 466  *
 467  * Returns:
 468  **/
 469 int hkl_lattice_get_B(const HklLattice *self, HklMatrix *B)
 470 {
 471         double D;
 472         double c_alpha, s_alpha;
 473         double c_beta, s_beta;
 474         double c_gamma, s_gamma;
 475         double b11, b22, tmp;
 476
 477         c_alpha = cos(hkl_parameter_value_get(self->alpha, HKL_UNIT_DEFAULT));
 478         c_beta = cos(hkl_parameter_value_get(self->beta, HKL_UNIT_DEFAULT));
 479         c_gamma = cos(hkl_parameter_value_get(self->gamma, HKL_UNIT_DEFAULT));
 480         D = 1 - c_alpha*c_alpha - c_beta*c_beta - c_gamma*c_gamma
 481                 + 2*c_alpha*c_beta*c_gamma;
 482
 483         if (D > 0.)
 484                 D = sqrt(D);
 485         else
 486                 return FALSE;
 487
 488         s_alpha = sin(hkl_parameter_value_get(self->alpha, HKL_UNIT_DEFAULT));
 489         s_beta  = sin(hkl_parameter_value_get(self->beta, HKL_UNIT_DEFAULT));
 490         s_gamma = sin(hkl_parameter_value_get(self->gamma, HKL_UNIT_DEFAULT));
 491
 492         b11 = HKL_TAU / (hkl_parameter_value_get(self->b, HKL_UNIT_DEFAULT) * s_alpha);
 493         b22 = HKL_TAU / hkl_parameter_value_get(self->c, HKL_UNIT_DEFAULT);
 494         tmp = b22 / s_alpha;
 495
 496         B->data[0][0] = HKL_TAU * s_alpha / (hkl_parameter_value_get(self->a, HKL_UNIT_DEFAULT) * D);
 497         B->data[0][1] = b11 / D * (c_alpha*c_beta - c_gamma);
 498         B->data[0][2] = tmp / D * (c_gamma*c_alpha - c_beta);
 499
 500         B->data[1][0] = 0;
 501         B->data[1][1] = b11;
 502         B->data[1][2] = tmp / (s_beta*s_gamma) * (c_beta*c_gamma - c_alpha);
 503
 504         B->data[2][0] = 0;
 505         B->data[2][1] = 0;
 506         B->data[2][2] = b22;
 507
 508         return TRUE;
 509 }
 510
 511 /**
 512  * hkl_lattice_get_1_B: (skip)
 513  * @self: the @HklLattice
 514  * @B: (out): where to store the 1/B matrix
 515  *
 516  * Compute the invert of B (needed by the hkl_sample_UB_set method)
 517  * should be optimized
 518  *
 519  * Returns: TRUE or FALSE depending of the success of the
 520  * computation.
 521  **/
 522 int hkl_lattice_get_1_B(const HklLattice *self, HklMatrix *B)
 523 {
 524         HklMatrix tmp;
 525         double a;
 526         double b;
 527         double c;
 528         double d;
 529         double e;
 530         double f;
 531
 532         if(!self || !B)
 533                 return FALSE;
 534
 535         /*
 536          * first compute the B matrix
 537          * | a b c |
 538          * | 0 d e |
 539          * | 0 0 f |
 540          */
 541         hkl_lattice_get_B(self, &tmp);
 542
 543         /*
 544          * now invert this triangular matrix
 545          */
 546         a = tmp.data[0][0];
 547         b = tmp.data[0][1];
 548         c = tmp.data[0][2];
 549         d = tmp.data[1][1];
 550         e = tmp.data[1][2];
 551         f = tmp.data[2][2];
 552
 553         B->data[0][0] = 1 / a;
 554         B->data[0][1] = -b / a / d;
 555         B->data[0][2] = (b * e - d * c) / a / d / f;
 556
 557         B->data[1][0] = 0;
 558         B->data[1][1] = 1 / d;
 559         B->data[1][2] = -e / d / f;
 560
 561         B->data[2][0] = 0;
 562         B->data[2][1] = 0;
 563         B->data[2][2] = 1 / f;
 564
 565         return TRUE;
 566 }
 567
 568 /**
 569  * hkl_lattice_reciprocal:
 570  * @self: the this ptr
 571  * @reciprocal: the lattice where the result will be computed
 572  *
 573  * compute the reciprocal #HklLattice and put the result id the
 574  * provided @reciprocal parameter
 575  *
 576  * Returns: 0 or 1 if it succeed.
 577  **/
 578 int hkl_lattice_reciprocal(const HklLattice *self, HklLattice *reciprocal)
 579 {
 580         double c_alpha, c_beta, c_gamma;
 581         double s_alpha, s_beta, s_gamma;
 582         double c_beta1, c_beta2, c_beta3;
 583         double s_beta1, s_beta2, s_beta3;
 584         double s_beta_s_gamma, s_gamma_s_alpha, s_alpha_s_beta;
 585         double D;
 586
 587         c_alpha = cos(hkl_parameter_value_get(self->alpha, HKL_UNIT_DEFAULT));
 588         c_beta  = cos(hkl_parameter_value_get(self->beta, HKL_UNIT_DEFAULT));
 589         c_gamma = cos(hkl_parameter_value_get(self->gamma, HKL_UNIT_DEFAULT));
 590         D = 1 - c_alpha*c_alpha - c_beta*c_beta - c_gamma*c_gamma
 591                 + 2*c_alpha*c_beta*c_gamma;
 592
 593         if (D > 0.)
 594                 D = sqrt(D);
 595         else
 596                 return FALSE;
 597
 598         s_alpha = sin(hkl_parameter_value_get(self->alpha, HKL_UNIT_DEFAULT));
 599         s_beta  = sin(hkl_parameter_value_get(self->beta, HKL_UNIT_DEFAULT));
 600         s_gamma = sin(hkl_parameter_value_get(self->gamma, HKL_UNIT_DEFAULT));
 601
 602         s_beta_s_gamma  = s_beta  * s_gamma;
 603         s_gamma_s_alpha = s_gamma * s_alpha;
 604         s_alpha_s_beta  = s_alpha * s_beta;
 605
 606         c_beta1 = (c_beta  * c_gamma - c_alpha) / s_beta_s_gamma;
 607         c_beta2 = (c_gamma * c_alpha - c_beta)  / s_gamma_s_alpha;
 608         c_beta3 = (c_alpha * c_beta  - c_gamma) / s_alpha_s_beta;
 609         s_beta1 = D / s_beta_s_gamma;
 610         s_beta2 = D / s_gamma_s_alpha;
 611         s_beta3 = D / s_alpha_s_beta;
 612
 613         hkl_lattice_set(reciprocal,
 614                         HKL_TAU * s_alpha / (hkl_parameter_value_get(self->a, HKL_UNIT_DEFAULT) * D),
 615                         HKL_TAU * s_beta  / (hkl_parameter_value_get(self->b, HKL_UNIT_DEFAULT) * D),
 616                         HKL_TAU * s_gamma / (hkl_parameter_value_get(self->c, HKL_UNIT_DEFAULT) * D),
 617                         atan2(s_beta1, c_beta1),
 618                         atan2(s_beta2, c_beta2),
 619                         atan2(s_beta3, c_beta3),
 620                         HKL_UNIT_DEFAULT, NULL);
 621
 622         return TRUE;
 623 }
 624
 625 /**
 626  * hkl_lattice_randomize: (skip)
 627  * @self:
 628  *
 629  * randomize the lattice
 630  **/
 631 void hkl_lattice_randomize(HklLattice *self)
 632 {
 633         static HklVector vector_x = {{1, 0, 0}};
 634         HklVector a, b, c;
 635         HklVector axe;
 636         unsigned int angles_to_randomize;
 637
 638         /* La valeur des angles alpha, beta et gamma ne sont pas indépendant. */
 639         /* Il faut donc gérer les différents cas. */
 640         hkl_parameter_randomize(self->a);
 641         hkl_parameter_randomize(self->b);
 642         hkl_parameter_randomize(self->c);
 643
 644         angles_to_randomize = self->alpha->fit
 645                 + self->beta->fit
 646                 + self->gamma->fit;
 647         switch (angles_to_randomize) {
 648         case 0:
 649                 break;
 650         case 1:
 651                 if (self->alpha->fit) {
 652                         /* alpha */
 653                         a = b = c = vector_x;
 654
 655                         /* randomize b */
 656                         hkl_vector_randomize_vector(&axe, &a);
 657                         hkl_vector_rotated_around_vector(&b, &axe,
 658                                                          hkl_parameter_value_get(self->gamma,
 659                                                                                  HKL_UNIT_DEFAULT));
 660
 661                         /* randomize c */
 662                         hkl_vector_randomize_vector(&axe, &a);
 663                         hkl_vector_rotated_around_vector(&c, &axe,
 664                                                          hkl_parameter_value_get(self->beta,
 665                                                                                  HKL_UNIT_DEFAULT));
 666
 667                         /* compute the alpha angle. */
 668                         hkl_parameter_value_set(self->alpha, hkl_vector_angle(&b, &c),
 669                                                 HKL_UNIT_DEFAULT, NULL);
 670                 } else if (self->beta->fit) {
 671                         /* beta */
 672                         a = b = vector_x;
 673
 674                         /* randomize b */
 675                         hkl_vector_randomize_vector(&axe, &a);
 676                         hkl_vector_rotated_around_vector(&b, &axe,
 677                                                          hkl_parameter_value_get(self->gamma,
 678                                                                                  HKL_UNIT_DEFAULT));
 679
 680                         /* randomize c */
 681                         c = b;
 682                         hkl_vector_randomize_vector(&axe, &b);
 683                         hkl_vector_rotated_around_vector(&c, &axe,
 684                                                          hkl_parameter_value_get(self->alpha,
 685                                                                                  HKL_UNIT_DEFAULT));
 686
 687                         /* compute beta */
 688                         hkl_parameter_value_set(self->beta, hkl_vector_angle(&a, &c),
 689                                                 HKL_UNIT_DEFAULT, NULL);
 690                 } else {
 691                         /* gamma */
 692                         a = c = vector_x;
 693
 694                         /* randomize c */
 695                         hkl_vector_randomize_vector(&axe, &a);
 696                         hkl_vector_rotated_around_vector(&c, &axe,
 697                                                          hkl_parameter_value_get(self->beta,
 698                                                                                  HKL_UNIT_DEFAULT));
 699
 700                         /* randomize b */
 701                         b = c;
 702                         hkl_vector_randomize_vector(&axe, &c);
 703                         hkl_vector_rotated_around_vector(&b, &axe,
 704                                                          hkl_parameter_value_get(self->alpha,
 705                                                                                  HKL_UNIT_DEFAULT));
 706
 707                         /* compute gamma */
 708                         hkl_parameter_value_set(self->gamma, hkl_vector_angle(&a, &b),
 709                                                 HKL_UNIT_DEFAULT, NULL);
 710                 }
 711                 break;
 712         case 2:
 713                 if (self->alpha->fit) {
 714                         if (self->beta->fit) {
 715                                 /* alpha + beta */
 716                                 a = b = vector_x;
 717
 718                                 /* randomize b */
 719                                 hkl_vector_randomize_vector(&axe, &a);
 720                                 hkl_vector_rotated_around_vector(&b, &axe,
 721                                                                  hkl_parameter_value_get(self->gamma,
 722                                                                                          HKL_UNIT_DEFAULT));
 723
 724                                 /* randomize c */
 725                                 hkl_vector_randomize_vector_vector(&c, &a, &b);
 726
 727                                 hkl_parameter_value_set(self->alpha, hkl_vector_angle(&b, &c),
 728                                                         HKL_UNIT_DEFAULT, NULL);
 729                                 hkl_parameter_value_set(self->beta, hkl_vector_angle(&a, &c),
 730                                                         HKL_UNIT_DEFAULT, NULL);
 731                         } else {
 732                                 /* alpha + gamma */
 733                                 a = c = vector_x;
 734
 735                                 /* randomize c */
 736                                 hkl_vector_randomize_vector(&axe, &a);
 737                                 hkl_vector_rotated_around_vector(&c, &axe,
 738                                                                  hkl_parameter_value_get(self->beta,
 739                                                                                          HKL_UNIT_DEFAULT));
 740
 741                                 /* randomize c */
 742                                 hkl_vector_randomize_vector_vector(&b, &a, &c);
 743
 744                                 hkl_parameter_value_set(self->alpha, hkl_vector_angle(&b, &c),
 745                                                         HKL_UNIT_DEFAULT, NULL);
 746                                 hkl_parameter_value_set(self->gamma, hkl_vector_angle(&a, &b),
 747                                                         HKL_UNIT_DEFAULT, NULL);
 748                         }
 749                 } else {
 750                         /* beta + gamma */
 751                         b = c = vector_x;
 752
 753                         /* randomize c */
 754                         hkl_vector_randomize_vector(&axe, &b);
 755                         hkl_vector_rotated_around_vector(&c, &axe,
 756                                                          hkl_parameter_value_get(self->alpha,
 757                                                                                  HKL_UNIT_DEFAULT));
 758
 759                         /* randomize c */
 760                         hkl_vector_randomize_vector_vector(&a, &b, &c);
 761
 762                         hkl_parameter_value_set(self->beta, hkl_vector_angle(&a, &c),
 763                                                 HKL_UNIT_DEFAULT, NULL);
 764                         hkl_parameter_value_set(self->gamma, hkl_vector_angle(&a, &b),
 765                                                 HKL_UNIT_DEFAULT, NULL);
 766                 }
 767                 break;
 768         case 3:
 769                 hkl_vector_randomize(&a);
 770                 hkl_vector_randomize_vector(&b, &a);
 771                 hkl_vector_randomize_vector_vector(&c, &b, &a);
 772
 773                 hkl_parameter_value_set(self->alpha, hkl_vector_angle(&b, &c),
 774                                         HKL_UNIT_DEFAULT, NULL);
 775                 hkl_parameter_value_set(self->beta, hkl_vector_angle(&a, &c),
 776                                         HKL_UNIT_DEFAULT, NULL);
 777                 hkl_parameter_value_set(self->gamma, hkl_vector_angle(&a, &b),
 778                                         HKL_UNIT_DEFAULT, NULL);
 779                 break;
 780         }
 781 }
 782
 783 /**
 784  * hkl_lattice_fprintf: (skip)
 785  * @f:
 786  * @self:
 787  *
 788  * print into a file the lattice.
 789  **/
 790 void hkl_lattice_fprintf(FILE *f, HklLattice const *self)
 791 {
 792         fprintf(f, "\n");
 793         hkl_parameter_fprintf(f, self->a);
 794         fprintf(f, "\n");
 795         hkl_parameter_fprintf(f, self->b);
 796         fprintf(f, "\n");
 797         hkl_parameter_fprintf(f, self->c);
 798         fprintf(f, "\n");
 799         hkl_parameter_fprintf(f, self->alpha);
 800         fprintf(f, "\n");
 801         hkl_parameter_fprintf(f, self->beta);
 802         fprintf(f, "\n");
 803         hkl_parameter_fprintf(f, self->gamma);
 804 }