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-2019, 2022, 2023, 2024 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 **err)
  88 {
  89         HklLattice *self = NULL;
  90         HklParameter *parameter = NULL;
  91         double volume;
  92
  93         g_return_val_if_fail (err == NULL || *err == NULL, NULL);
  94
  95         if(!check_lattice_param(a, b, c, alpha, beta, gamma, &volume, err))
  96         {
  97                 hkl_assert (err == NULL || *err != NULL);
  98                 return FALSE;
  99         }
 100         hkl_assert (err == NULL || *err == NULL);
 101
 102         self = g_new(HklLattice, 1);
 103
 104         /* a */
 105         parameter = hkl_parameter_new("a", "The length of the first lattice vector",
 106                                       0, a, a+10,
 107                                       TRUE, TRUE,
 108                                       &hkl_unit_length_nm,
 109                                       &hkl_unit_length_nm, err);
 110         if (!parameter) {
 111                 hkl_assert (err == NULL || *err != NULL);
 112                 goto failed;
 113         }else{
 114                 self->a = parameter;
 115         }
 116
 117         /* b */
 118         parameter = hkl_parameter_new("b", "The length of the second lattice vector",
 119                                       0, b, b+10,
 120                                       TRUE, TRUE,
 121                                       &hkl_unit_length_nm,
 122                                       &hkl_unit_length_nm,
 123                                       err);
 124         if (!parameter) {
 125                 hkl_assert (err == NULL || *err != NULL);
 126                 goto free_a;
 127         }else{
 128                 self->b = parameter;
 129         }
 130
 131         /* c */
 132         parameter = hkl_parameter_new("c", "The length of the third lattice vector",
 133                                       0, c, c+10,
 134                                       TRUE, TRUE,
 135                                       &hkl_unit_length_nm,
 136                                       &hkl_unit_length_nm,
 137                                       err);
 138         if (!parameter) {
 139                 hkl_assert (err == NULL || *err != NULL);
 140                 goto free_b;
 141         }else{
 142                 self->c = parameter;
 143         }
 144
 145         /* alpha */
 146         parameter = hkl_parameter_new("alpha",
 147                                       "The angle between the second and third lattice vector",
 148                                       -M_PI, alpha, M_PI,
 149                                       TRUE, TRUE,
 150                                       &hkl_unit_angle_rad,
 151                                       &hkl_unit_angle_deg,
 152                                       err);
 153         if (!parameter) {
 154                 hkl_assert (err == NULL || *err != NULL);
 155                 goto free_c;
 156         }else{
 157                 self->alpha = parameter;
 158         }
 159
 160         /* beta */
 161         parameter = hkl_parameter_new("beta",
 162                                       "The angle between the first and third lattice vector",
 163                                       -M_PI, beta, M_PI,
 164                                       TRUE, TRUE,
 165                                       &hkl_unit_angle_rad,
 166                                       &hkl_unit_angle_deg,
 167                                       err);
 168         if (!parameter) {
 169                 hkl_assert (err == NULL || *err != NULL);
 170                 goto free_alpha;
 171         }else{
 172                 self->beta = parameter;
 173         }
 174
 175         /* gamma */
 176         parameter = hkl_parameter_new("gamma",
 177                                       "The angle between the first and second lattice vector",
 178                                       -M_PI, gamma, M_PI,
 179                                       TRUE, TRUE,
 180                                       &hkl_unit_angle_rad,
 181                                       &hkl_unit_angle_deg,
 182                                       err);
 183         if (!parameter) {
 184                 hkl_assert (err == NULL || *err != NULL);
 185                 goto free_beta;
 186         }else{
 187                 self->gamma = parameter;
 188         }
 189
 190         /* volume */
 191         parameter = hkl_parameter_new("volume",
 192                                       "The volume of the lattice",
 193                                       0, volume, a*b*c,
 194                                       FALSE, FALSE,
 195                                       &hkl_unit_length_nm,
 196                                       &hkl_unit_length_nm,
 197                                       err);
 198         if (!parameter) {
 199                 hkl_assert (err == NULL || *err != NULL);
 200                 goto free_gamma;
 201         }else{
 202                 self->volume = parameter;
 203         }
 204
 205         return self;
 206
 207 free_gamma:
 208         hkl_parameter_free(self->gamma);
 209 free_beta:
 210         hkl_parameter_free(self->beta);
 211 free_alpha:
 212         hkl_parameter_free(self->alpha);
 213 free_c:
 214         hkl_parameter_free(self->c);
 215 free_b:
 216         hkl_parameter_free(self->b);
 217 free_a:
 218         hkl_parameter_free(self->a);
 219 failed:
 220         return FALSE;
 221 }
 222
 223 /**
 224  * hkl_lattice_new_copy: (skip)
 225  * @self:
 226  *
 227  * copy constructor
 228  *
 229  * Returns:
 230  **/
 231 HklLattice *hkl_lattice_new_copy(const HklLattice *self)
 232 {
 233         HklLattice *copy = g_new(HklLattice, 1);
 234
 235         copy->a = hkl_parameter_new_copy(self->a);
 236         copy->b = hkl_parameter_new_copy(self->b);
 237         copy->c = hkl_parameter_new_copy(self->c);
 238         copy->alpha = hkl_parameter_new_copy(self->alpha);
 239         copy->beta = hkl_parameter_new_copy(self->beta);
 240         copy->gamma = hkl_parameter_new_copy(self->gamma);
 241         copy->volume = hkl_parameter_new_copy(self->volume);
 242
 243         return copy;
 244 }
 245
 246 /**
 247  * hkl_lattice_new_default: (skip)
 248  *
 249  * default constructor
 250  *
 251  * Returns:
 252  **/
 253 HklLattice* hkl_lattice_new_default(void)
 254 {
 255         return hkl_lattice_new(1.54, 1.54, 1.54,
 256                                90*HKL_DEGTORAD, 90*HKL_DEGTORAD, 90*HKL_DEGTORAD,
 257                                NULL);
 258 }
 259
 260 /**
 261  * hkl_lattice_free: (skip)
 262  * @self:
 263  *
 264  * destructor
 265  **/
 266 void hkl_lattice_free(HklLattice *self)
 267 {
 268         hkl_parameter_free(self->a);
 269         hkl_parameter_free(self->b);
 270         hkl_parameter_free(self->c);
 271         hkl_parameter_free(self->alpha);
 272         hkl_parameter_free(self->beta);
 273         hkl_parameter_free(self->gamma);
 274         hkl_parameter_free(self->volume);
 275         free(self);
 276 }
 277
 278 #define HKL_LATTICE_X_SET(_p, _parameter, _error) do{                   \
 279                 hkl_error ((_error) == NULL || *(_error) == NULL);      \
 280                 double a, b, c, alpha, beta, gamma;                     \
 281                 /* check if the combinaison of parameters is ok */      \
 282                 hkl_lattice_get(self, &a, &b, &c, &alpha, &beta, &gamma, HKL_UNIT_DEFAULT); \
 283                 _p = hkl_parameter_value_get((_parameter), HKL_UNIT_DEFAULT); \
 284                                                                         \
 285                 if(!hkl_lattice_set(self, a, b, c, alpha, beta, gamma, HKL_UNIT_DEFAULT, (_error))){ \
 286                         hkl_assert ((_error) == NULL || *(_error) != NULL); \
 287                         return FALSE;                                   \
 288                 }                                                       \
 289                 hkl_assert ((_error) == NULL || *(_error) == NULL);     \
 290                 return hkl_parameter_init_copy(self->_p, (_parameter), (_error)); \
 291         }while(0)
 292
 293 /**
 294  * hkl_lattice_a_get:
 295  * @self: the this ptr
 296  **/
 297 const HklParameter *hkl_lattice_a_get(const HklLattice *self)
 298 {
 299         return self->a;
 300 }
 301
 302 /**
 303  * hkl_lattice_a_set:
 304  * @self: the this ptr
 305  * @parameter: the parameter to set
 306  * @error: return location for a GError, or NULL
 307  *
 308  * Returns: TRUE on success, FALSE if an error occurred
 309  **/
 310 int hkl_lattice_a_set(HklLattice *self, const HklParameter *parameter,
 311                       GError **error)
 312 {
 313         HKL_LATTICE_X_SET(a, parameter, error);
 314 }
 315
 316 /**
 317  * hkl_lattice_b_get:
 318  * @self: the this ptr
 319  **/
 320 const HklParameter *hkl_lattice_b_get(const HklLattice *self)
 321 {
 322         return self->b;
 323 }
 324
 325 /**
 326  * hkl_lattice_b_set:
 327  * @self: the this ptr
 328  * @parameter: the parameter to set
 329  * @error: return location for a GError, or NULL
 330  *
 331  * Returns: TRUE on success, FALSE if an error occurred
 332  **/
 333 int hkl_lattice_b_set(HklLattice *self, const HklParameter *parameter,
 334                       GError **error)
 335 {
 336         HKL_LATTICE_X_SET(b, parameter, error);
 337 }
 338
 339 /**
 340  * hkl_lattice_c_get:
 341  * @self: the this ptr
 342  **/
 343 const HklParameter *hkl_lattice_c_get(const HklLattice *self)
 344 {
 345         return self->c;
 346 }
 347
 348 /**
 349  * hkl_lattice_c_set:
 350  * @self: the this ptr
 351  * @parameter: the parameter to set
 352  * @error: return location for a GError, or NULL
 353  *
 354  * Returns: TRUE on success, FALSE if an error occurred
 355  **/
 356 int hkl_lattice_c_set(HklLattice *self, const HklParameter *parameter,
 357                       GError **error)
 358 {
 359         HKL_LATTICE_X_SET(c, parameter, error);
 360 }
 361
 362 /**
 363  * hkl_lattice_alpha_get:
 364  * @self: the this ptr
 365  **/
 366 const HklParameter *hkl_lattice_alpha_get(const HklLattice *self)
 367 {
 368         return self->alpha;
 369 }
 370
 371 /**
 372  * hkl_lattice_alpha_set:
 373  * @self: the this ptr
 374  * @parameter: the parameter to set
 375  * @error: return location for a GError, or NULL
 376  *
 377  * Returns: TRUE on success, FALSE if an error occurred
 378  **/
 379 int hkl_lattice_alpha_set(HklLattice *self, const HklParameter *parameter,
 380                           GError **error)
 381 {
 382         HKL_LATTICE_X_SET(alpha, parameter, error);
 383 }
 384
 385 /**
 386  * hkl_lattice_beta_get:
 387  * @self: the this ptr
 388  **/
 389 const HklParameter *hkl_lattice_beta_get(const HklLattice *self)
 390 {
 391         return self->beta;
 392 }
 393
 394 /**
 395  * hkl_lattice_beta_set:
 396  * @self: the this ptr
 397  * @parameter: the parameter to set
 398  * @error: return location for a GError, or NULL
 399  *
 400  * Returns: TRUE on success, FALSE if an error occurred
 401  **/
 402 int hkl_lattice_beta_set(HklLattice *self, const HklParameter *parameter,
 403                          GError **error)
 404 {
 405         HKL_LATTICE_X_SET(beta, parameter, error);
 406 }
 407
 408 /**
 409  * hkl_lattice_gamma_get:
 410  * @self: the this ptr
 411  **/
 412 const HklParameter *hkl_lattice_gamma_get(const HklLattice *self)
 413 {
 414         return self->gamma;
 415 }
 416
 417 /**
 418  * hkl_lattice_gamma_set:
 419  * @self: the this ptr
 420  * @parameter: the parameter to set
 421  * @error: return location for a GError, or NULL
 422  *
 423  * Returns: TRUE on success, FALSE if an error occurred
 424  **/
 425 int hkl_lattice_gamma_set(HklLattice *self, const HklParameter *parameter,
 426                           GError **error)
 427 {
 428         HKL_LATTICE_X_SET(gamma, parameter, error);
 429 }
 430
 431 /**
 432  * hkl_lattice_volume_get:
 433  * @self: the this ptr
 434  **/
 435 const HklParameter *hkl_lattice_volume_get(const HklLattice *self)
 436 {
 437         return self->volume;
 438 }
 439
 440 /**
 441  * hkl_lattice_lattice_set: (skip)
 442  * @self: the this ptr
 443  * @lattice: the lattice to set from.
 444  **/
 445 void hkl_lattice_lattice_set(HklLattice *self, const HklLattice *lattice)
 446 {
 447         if (self == lattice)
 448                 return;
 449
 450         hkl_parameter_init_copy(self->a, lattice->a, NULL);
 451         hkl_parameter_init_copy(self->b, lattice->b, NULL);
 452         hkl_parameter_init_copy(self->c, lattice->c, NULL);
 453         hkl_parameter_init_copy(self->alpha, lattice->alpha, NULL);
 454         hkl_parameter_init_copy(self->beta, lattice->beta, NULL);
 455         hkl_parameter_init_copy(self->gamma, lattice->gamma, NULL);
 456         hkl_parameter_init_copy(self->volume, lattice->volume, NULL);
 457 }
 458
 459 /**
 460  * hkl_lattice_set:
 461  * @self:
 462  * @a:
 463  * @b:
 464  * @c:
 465  * @alpha:
 466  * @beta:
 467  * @gamma:
 468  *
 469  * set the lattice parameters
 470  *
 471  * Returns:
 472  **/
 473 int hkl_lattice_set(HklLattice *self,
 474                     double a, double b, double c,
 475                     double alpha, double beta, double gamma,
 476                     HklUnitEnum unit_type, GError **error)
 477 {
 478         hkl_error (error == NULL || *error == NULL);
 479
 480         double new_a, previous_a;
 481         double new_b, previous_b;
 482         double new_c, previous_c;
 483         double new_alpha, previous_alpha;
 484         double new_beta, previous_beta;
 485         double new_gamma, previous_gamma;
 486         double new_volume;
 487
 488         new_a = convert_to_default(self->a, a, unit_type);
 489         new_b = convert_to_default(self->b, b, unit_type);
 490         new_c = convert_to_default(self->c, c, unit_type);
 491         new_alpha = convert_to_default(self->alpha, alpha, unit_type);
 492         new_beta = convert_to_default(self->beta, beta, unit_type);
 493         new_gamma = convert_to_default(self->gamma, gamma, unit_type);
 494
 495         /* need to do the conversion before the check */
 496         if(!check_lattice_param(new_a, new_b, new_c, new_alpha, new_beta, new_gamma, &new_volume, error)){
 497                 hkl_assert (error == NULL || *error != NULL);
 498                 return FALSE;
 499         }
 500         hkl_assert (error == NULL || *error == NULL);
 501
 502         /* check a */
 503         if(!hkl_parameter_value_set(self->a, new_a, HKL_UNIT_DEFAULT, NULL)){
 504                 hkl_assert (error == NULL || *error != NULL);
 505                 goto failed;
 506         }
 507         hkl_assert (error == NULL || *error == NULL);
 508
 509         if(!hkl_parameter_value_set(self->b, new_b, HKL_UNIT_DEFAULT, NULL)){
 510                 hkl_assert (error == NULL || *error != NULL);
 511                 goto restore_a;
 512         }
 513         hkl_assert (error == NULL || *error == NULL);
 514
 515         if(!hkl_parameter_value_set(self->c, new_c, HKL_UNIT_DEFAULT, NULL)){
 516                 hkl_assert (error == NULL || *error != NULL);
 517                 goto restore_b;
 518         }
 519         hkl_assert (error == NULL || *error == NULL);
 520
 521         if(!hkl_parameter_value_set(self->alpha, new_alpha, HKL_UNIT_DEFAULT, NULL)){
 522                 hkl_assert (error == NULL || *error != NULL);
 523                 goto restore_c;
 524         }
 525         hkl_assert (error == NULL || *error == NULL);
 526
 527         if(!hkl_parameter_value_set(self->beta, new_beta, HKL_UNIT_DEFAULT, NULL)){
 528                 hkl_assert (error == NULL || *error != NULL);
 529                 goto restore_alpha;
 530         }
 531         hkl_assert (error == NULL || *error == NULL);
 532
 533         if(!hkl_parameter_value_set(self->gamma, new_gamma, HKL_UNIT_DEFAULT, NULL)){
 534                 hkl_assert (error == NULL || *error != NULL);
 535                 goto restore_beta;
 536         }
 537         hkl_assert (error == NULL || *error == NULL);
 538
 539         if(!hkl_parameter_value_set(self->volume, new_volume, HKL_UNIT_DEFAULT, NULL)){
 540                 hkl_assert (error == NULL || *error != NULL);
 541                 goto restore_gamma;
 542         }
 543         hkl_assert (error == NULL || *error == NULL);
 544
 545         return TRUE;
 546
 547 restore_gamma:
 548         IGNORE(hkl_parameter_value_set(self->gamma, previous_gamma,
 549                                        HKL_UNIT_DEFAULT, NULL));
 550 restore_beta:
 551         IGNORE(hkl_parameter_value_set(self->beta, previous_beta,
 552                                        HKL_UNIT_DEFAULT, NULL));
 553 restore_alpha:
 554         IGNORE(hkl_parameter_value_set(self->alpha, previous_alpha,
 555                                        HKL_UNIT_DEFAULT, NULL));
 556 restore_c:
 557         IGNORE(hkl_parameter_value_set(self->c, previous_c,
 558                                        HKL_UNIT_DEFAULT, NULL));
 559 restore_b:
 560         IGNORE(hkl_parameter_value_set(self->b, previous_b,
 561                                        HKL_UNIT_DEFAULT, NULL));
 562 restore_a:
 563         IGNORE(hkl_parameter_value_set(self->a, previous_a,
 564                                        HKL_UNIT_DEFAULT, NULL));
 565 failed:
 566         return FALSE;
 567 }
 568
 569 /**
 570  * hkl_lattice_get:
 571  * @self:
 572  * @a: (out caller-allocates):
 573  * @b: (out caller-allocates):
 574  * @c: (out caller-allocates):
 575  * @alpha: (out caller-allocates):
 576  * @beta: (out caller-allocates):
 577  * @gamma: (out caller-allocates):
 578  *
 579  * get the lattice parameters
 580  **/
 581 void hkl_lattice_get(const HklLattice *self,
 582                      double *a, double *b, double *c,
 583                      double *alpha, double *beta, double *gamma,
 584                      HklUnitEnum unit_type)
 585 {
 586         *a = hkl_parameter_value_get(self->a, unit_type);
 587         *b = hkl_parameter_value_get(self->b, unit_type);
 588         *c = hkl_parameter_value_get(self->c, unit_type);
 589         *alpha = hkl_parameter_value_get(self->alpha, unit_type);
 590         *beta = hkl_parameter_value_get(self->beta, unit_type);
 591         *gamma = hkl_parameter_value_get(self->gamma, unit_type);
 592 }
 593
 594 /**
 595  * hkl_lattice_get_B: (skip)
 596  * @self:
 597  * @B: (out): where to store the B matrix
 598  *
 599  * Get the B matrix from the lattice parameters
 600  *
 601  * Returns:
 602  **/
 603 int hkl_lattice_get_B(const HklLattice *self, HklMatrix *B)
 604 {
 605         double D;
 606         double c_alpha, s_alpha;
 607         double c_beta, s_beta;
 608         double c_gamma, s_gamma;
 609         double b11, b22, tmp;
 610
 611         c_alpha = cos(hkl_parameter_value_get(self->alpha, HKL_UNIT_DEFAULT));
 612         c_beta = cos(hkl_parameter_value_get(self->beta, HKL_UNIT_DEFAULT));
 613         c_gamma = cos(hkl_parameter_value_get(self->gamma, HKL_UNIT_DEFAULT));
 614         D = 1 - c_alpha*c_alpha - c_beta*c_beta - c_gamma*c_gamma
 615                 + 2*c_alpha*c_beta*c_gamma;
 616
 617         if (D > 0.)
 618                 D = sqrt(D);
 619         else
 620                 return FALSE;
 621
 622         s_alpha = sin(hkl_parameter_value_get(self->alpha, HKL_UNIT_DEFAULT));
 623         s_beta  = sin(hkl_parameter_value_get(self->beta, HKL_UNIT_DEFAULT));
 624         s_gamma = sin(hkl_parameter_value_get(self->gamma, HKL_UNIT_DEFAULT));
 625
 626         b11 = HKL_TAU / (hkl_parameter_value_get(self->b, HKL_UNIT_DEFAULT) * s_alpha);
 627         b22 = HKL_TAU / hkl_parameter_value_get(self->c, HKL_UNIT_DEFAULT);
 628         tmp = b22 / s_alpha;
 629
 630         B->data[0][0] = HKL_TAU * s_alpha / (hkl_parameter_value_get(self->a, HKL_UNIT_DEFAULT) * D);
 631         B->data[0][1] = b11 / D * (c_alpha*c_beta - c_gamma);
 632         B->data[0][2] = tmp / D * (c_gamma*c_alpha - c_beta);
 633
 634         B->data[1][0] = 0;
 635         B->data[1][1] = b11;
 636         B->data[1][2] = tmp / (s_beta*s_gamma) * (c_beta*c_gamma - c_alpha);
 637
 638         B->data[2][0] = 0;
 639         B->data[2][1] = 0;
 640         B->data[2][2] = b22;
 641
 642         return TRUE;
 643 }
 644
 645 /**
 646  * hkl_lattice_get_1_B: (skip)
 647  * @self: the @HklLattice
 648  * @B: (out): where to store the 1/B matrix
 649  *
 650  * Compute the invert of B (needed by the hkl_sample_UB_set method)
 651  * should be optimized
 652  *
 653  * Returns: TRUE or FALSE depending of the success of the
 654  * computation.
 655  **/
 656 int hkl_lattice_get_1_B(const HklLattice *self, HklMatrix *B)
 657 {
 658         HklMatrix tmp;
 659         double a;
 660         double b;
 661         double c;
 662         double d;
 663         double e;
 664         double f;
 665
 666         /*
 667          * first compute the B matrix
 668          * | a b c |
 669          * | 0 d e |
 670          * | 0 0 f |
 671          */
 672         hkl_lattice_get_B(self, &tmp);
 673
 674         /*
 675          * now invert this triangular matrix
 676          */
 677         a = tmp.data[0][0];
 678         b = tmp.data[0][1];
 679         c = tmp.data[0][2];
 680         d = tmp.data[1][1];
 681         e = tmp.data[1][2];
 682         f = tmp.data[2][2];
 683
 684         B->data[0][0] = 1 / a;
 685         B->data[0][1] = -b / a / d;
 686         B->data[0][2] = (b * e - d * c) / a / d / f;
 687
 688         B->data[1][0] = 0;
 689         B->data[1][1] = 1 / d;
 690         B->data[1][2] = -e / d / f;
 691
 692         B->data[2][0] = 0;
 693         B->data[2][1] = 0;
 694         B->data[2][2] = 1 / f;
 695
 696         return TRUE;
 697 }
 698
 699 /**
 700  * hkl_lattice_reciprocal:
 701  * @self: the this ptr
 702  * @reciprocal: the lattice where the result will be computed
 703  *
 704  * compute the reciprocal #HklLattice and put the result id the
 705  * provided @reciprocal parameter
 706  *
 707  * Returns: 0 or 1 if it succeed.
 708  **/
 709 int hkl_lattice_reciprocal(const HklLattice *self, HklLattice *reciprocal)
 710 {
 711         double c_alpha, c_beta, c_gamma;
 712         double s_alpha, s_beta, s_gamma;
 713         double c_beta1, c_beta2, c_beta3;
 714         double s_beta1, s_beta2, s_beta3;
 715         double s_beta_s_gamma, s_gamma_s_alpha, s_alpha_s_beta;
 716         double D;
 717
 718         c_alpha = cos(hkl_parameter_value_get(self->alpha, HKL_UNIT_DEFAULT));
 719         c_beta  = cos(hkl_parameter_value_get(self->beta, HKL_UNIT_DEFAULT));
 720         c_gamma = cos(hkl_parameter_value_get(self->gamma, HKL_UNIT_DEFAULT));
 721         D = 1 - c_alpha*c_alpha - c_beta*c_beta - c_gamma*c_gamma
 722                 + 2*c_alpha*c_beta*c_gamma;
 723
 724         if (D > 0.)
 725                 D = sqrt(D);
 726         else
 727                 return FALSE;
 728
 729         s_alpha = sin(hkl_parameter_value_get(self->alpha, HKL_UNIT_DEFAULT));
 730         s_beta  = sin(hkl_parameter_value_get(self->beta, HKL_UNIT_DEFAULT));
 731         s_gamma = sin(hkl_parameter_value_get(self->gamma, HKL_UNIT_DEFAULT));
 732
 733         s_beta_s_gamma  = s_beta  * s_gamma;
 734         s_gamma_s_alpha = s_gamma * s_alpha;
 735         s_alpha_s_beta  = s_alpha * s_beta;
 736
 737         c_beta1 = (c_beta  * c_gamma - c_alpha) / s_beta_s_gamma;
 738         c_beta2 = (c_gamma * c_alpha - c_beta)  / s_gamma_s_alpha;
 739         c_beta3 = (c_alpha * c_beta  - c_gamma) / s_alpha_s_beta;
 740         s_beta1 = D / s_beta_s_gamma;
 741         s_beta2 = D / s_gamma_s_alpha;
 742         s_beta3 = D / s_alpha_s_beta;
 743
 744         IGNORE(hkl_lattice_set(reciprocal,
 745                                HKL_TAU * s_alpha / (hkl_parameter_value_get(self->a, HKL_UNIT_DEFAULT) * D),
 746                                HKL_TAU * s_beta  / (hkl_parameter_value_get(self->b, HKL_UNIT_DEFAULT) * D),
 747                                HKL_TAU * s_gamma / (hkl_parameter_value_get(self->c, HKL_UNIT_DEFAULT) * D),
 748                                atan2(s_beta1, c_beta1),
 749                                atan2(s_beta2, c_beta2),
 750                                atan2(s_beta3, c_beta3),
 751                                HKL_UNIT_DEFAULT, NULL));
 752
 753         return TRUE;
 754 }
 755
 756 /**
 757  * hkl_lattice_randomize: (skip)
 758  * @self:
 759  *
 760  * randomize the lattice
 761  **/
 762 void hkl_lattice_randomize(HklLattice *self)
 763 {
 764         static HklVector vector_x = {{1, 0, 0}};
 765         HklVector a, b, c;
 766         HklVector axe;
 767         unsigned int angles_to_randomize;
 768
 769         /* La valeur des angles alpha, beta et gamma ne sont pas indépendant. */
 770         /* Il faut donc gérer les différents cas. */
 771         hkl_parameter_randomize(self->a);
 772         hkl_parameter_randomize(self->b);
 773         hkl_parameter_randomize(self->c);
 774
 775         angles_to_randomize = self->alpha->fit
 776                 + self->beta->fit
 777                 + self->gamma->fit;
 778         switch (angles_to_randomize) {
 779         case 0:
 780                 break;
 781         case 1:
 782                 if (self->alpha->fit) {
 783                         /* alpha */
 784                         a = b = c = vector_x;
 785
 786                         /* randomize b */
 787                         hkl_vector_randomize_vector(&axe, &a);
 788                         hkl_vector_rotated_around_vector(&b, &axe,
 789                                                          hkl_parameter_value_get(self->gamma,
 790                                                                                  HKL_UNIT_DEFAULT));
 791
 792                         /* randomize c */
 793                         hkl_vector_randomize_vector(&axe, &a);
 794                         hkl_vector_rotated_around_vector(&c, &axe,
 795                                                          hkl_parameter_value_get(self->beta,
 796                                                                                  HKL_UNIT_DEFAULT));
 797
 798                         /* compute the alpha angle. */
 799                         IGNORE(hkl_parameter_value_set(self->alpha, hkl_vector_angle(&b, &c),
 800                                                        HKL_UNIT_DEFAULT, NULL));
 801                 } else if (self->beta->fit) {
 802                         /* beta */
 803                         a = b = vector_x;
 804
 805                         /* randomize b */
 806                         hkl_vector_randomize_vector(&axe, &a);
 807                         hkl_vector_rotated_around_vector(&b, &axe,
 808                                                          hkl_parameter_value_get(self->gamma,
 809                                                                                  HKL_UNIT_DEFAULT));
 810
 811                         /* randomize c */
 812                         c = b;
 813                         hkl_vector_randomize_vector(&axe, &b);
 814                         hkl_vector_rotated_around_vector(&c, &axe,
 815                                                          hkl_parameter_value_get(self->alpha,
 816                                                                                  HKL_UNIT_DEFAULT));
 817
 818                         /* compute beta */
 819                         IGNORE(hkl_parameter_value_set(self->beta, hkl_vector_angle(&a, &c),
 820                                                        HKL_UNIT_DEFAULT, NULL));
 821                 } else {
 822                         /* gamma */
 823                         a = c = vector_x;
 824
 825                         /* randomize c */
 826                         hkl_vector_randomize_vector(&axe, &a);
 827                         hkl_vector_rotated_around_vector(&c, &axe,
 828                                                          hkl_parameter_value_get(self->beta,
 829                                                                                  HKL_UNIT_DEFAULT));
 830
 831                         /* randomize b */
 832                         b = c;
 833                         hkl_vector_randomize_vector(&axe, &c);
 834                         hkl_vector_rotated_around_vector(&b, &axe,
 835                                                          hkl_parameter_value_get(self->alpha,
 836                                                                                  HKL_UNIT_DEFAULT));
 837
 838                         /* compute gamma */
 839                         IGNORE(hkl_parameter_value_set(self->gamma, hkl_vector_angle(&a, &b),
 840                                                        HKL_UNIT_DEFAULT, NULL));
 841                 }
 842                 break;
 843         case 2:
 844                 if (self->alpha->fit) {
 845                         if (self->beta->fit) {
 846                                 /* alpha + beta */
 847                                 a = b = vector_x;
 848
 849                                 /* randomize b */
 850                                 hkl_vector_randomize_vector(&axe, &a);
 851                                 hkl_vector_rotated_around_vector(&b, &axe,
 852                                                                  hkl_parameter_value_get(self->gamma,
 853                                                                                          HKL_UNIT_DEFAULT));
 854
 855                                 /* randomize c */
 856                                 hkl_vector_randomize_vector_vector(&c, &a, &b);
 857
 858                                 IGNORE(hkl_parameter_value_set(self->alpha, hkl_vector_angle(&b, &c),
 859                                                                HKL_UNIT_DEFAULT, NULL));
 860                                 IGNORE(hkl_parameter_value_set(self->beta, hkl_vector_angle(&a, &c),
 861                                                                HKL_UNIT_DEFAULT, NULL));
 862                         } else {
 863                                 /* alpha + gamma */
 864                                 a = c = vector_x;
 865
 866                                 /* randomize c */
 867                                 hkl_vector_randomize_vector(&axe, &a);
 868                                 hkl_vector_rotated_around_vector(&c, &axe,
 869                                                                  hkl_parameter_value_get(self->beta,
 870                                                                                          HKL_UNIT_DEFAULT));
 871
 872                                 /* randomize c */
 873                                 hkl_vector_randomize_vector_vector(&b, &a, &c);
 874
 875                                 IGNORE(hkl_parameter_value_set(self->alpha, hkl_vector_angle(&b, &c),
 876                                                                HKL_UNIT_DEFAULT, NULL));
 877                                 IGNORE(hkl_parameter_value_set(self->gamma, hkl_vector_angle(&a, &b),
 878                                                                HKL_UNIT_DEFAULT, NULL));
 879                         }
 880                 } else {
 881                         /* beta + gamma */
 882                         b = c = vector_x;
 883
 884                         /* randomize c */
 885                         hkl_vector_randomize_vector(&axe, &b);
 886                         hkl_vector_rotated_around_vector(&c, &axe,
 887                                                          hkl_parameter_value_get(self->alpha,
 888                                                                                  HKL_UNIT_DEFAULT));
 889
 890                         /* randomize c */
 891                         hkl_vector_randomize_vector_vector(&a, &b, &c);
 892
 893                         IGNORE(hkl_parameter_value_set(self->beta, hkl_vector_angle(&a, &c),
 894                                                        HKL_UNIT_DEFAULT, NULL));
 895                         IGNORE(hkl_parameter_value_set(self->gamma, hkl_vector_angle(&a, &b),
 896                                                        HKL_UNIT_DEFAULT, NULL));
 897                 }
 898                 break;
 899         case 3:
 900                 hkl_vector_randomize(&a);
 901                 hkl_vector_randomize_vector(&b, &a);
 902                 hkl_vector_randomize_vector_vector(&c, &b, &a);
 903
 904                 IGNORE(hkl_parameter_value_set(self->alpha, hkl_vector_angle(&b, &c),
 905                                                HKL_UNIT_DEFAULT, NULL));
 906                 IGNORE(hkl_parameter_value_set(self->beta, hkl_vector_angle(&a, &c),
 907                                                HKL_UNIT_DEFAULT, NULL));
 908                 IGNORE(hkl_parameter_value_set(self->gamma, hkl_vector_angle(&a, &b),
 909                                                HKL_UNIT_DEFAULT, NULL));
 910                 break;
 911         }
 912 }
 913
 914 /**
 915  * hkl_lattice_fprintf: (skip)
 916  * @f:
 917  * @self:
 918  *
 919  * print into a file the lattice.
 920  **/
 921 void hkl_lattice_fprintf(FILE *f, HklLattice const *self)
 922 {
 923         fprintf(f, "\n");
 924         hkl_parameter_fprintf(f, self->a);
 925         fprintf(f, "\n");
 926         hkl_parameter_fprintf(f, self->b);
 927         fprintf(f, "\n");
 928         hkl_parameter_fprintf(f, self->c);
 929         fprintf(f, "\n");
 930         hkl_parameter_fprintf(f, self->alpha);
 931         fprintf(f, "\n");
 932         hkl_parameter_fprintf(f, self->beta);
 933         fprintf(f, "\n");
 934         hkl_parameter_fprintf(f, self->gamma);
 935 }