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-2010 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 <stdlib.h>
  23 #include <math.h>
  24
  25 #include "hkl/hkl-lattice.h"
  26 #include "hkl/hkl-unit.h"
  27
  28 /* private */
  29
  30 static int check_lattice_param(double a, double b, double c,
  31                                double alpha, double beta, double gamma)
  32 {
  33         double D = 1. - cos(alpha)*cos(alpha) - cos(beta)*cos(beta)
  34                 - cos(gamma)*cos(gamma) + 2. * cos(alpha)*cos(beta)*cos(gamma);
  35
  36         if (D < 0.)
  37                 return HKL_FAIL;
  38         else
  39                 return HKL_SUCCESS;
  40 }
  41
  42 /* public */
  43
  44 HklLattice *hkl_lattice_new(double a, double b, double c,
  45                             double alpha, double beta, double gamma)
  46 {
  47         HklLattice *self = NULL;
  48         if(!check_lattice_param(a, b, c, alpha, beta, gamma)) {
  49                 self = HKL_MALLOC(HklLattice);
  50
  51                 self->a = hkl_parameter_new("a", 0, a, a+10,
  52                                             HKL_TRUE, HKL_TRUE,
  53                                             &hkl_unit_length_nm,
  54                                             &hkl_unit_length_nm);
  55                 self->b = hkl_parameter_new("b", 0, b, b+10,
  56                                             HKL_TRUE, HKL_TRUE,
  57                                             &hkl_unit_length_nm,
  58                                             &hkl_unit_length_nm);
  59                 self->c = hkl_parameter_new("c", 0, c, c+10,
  60                                             HKL_TRUE, HKL_TRUE,
  61                                             &hkl_unit_length_nm,
  62                                             &hkl_unit_length_nm);
  63                 self->alpha = hkl_parameter_new("alpha", -M_PI, alpha, M_PI,
  64                                                 HKL_TRUE, HKL_TRUE,
  65                                                 &hkl_unit_angle_rad,
  66                                                 &hkl_unit_angle_deg);
  67                 self->beta = hkl_parameter_new("beta", -M_PI, beta, M_PI,
  68                                                HKL_TRUE, HKL_TRUE,
  69                                                &hkl_unit_angle_rad,
  70                                                &hkl_unit_angle_deg);
  71                 self->gamma = hkl_parameter_new("gamma", -M_PI, gamma, M_PI,
  72                                                 HKL_TRUE, HKL_TRUE,
  73                                                 &hkl_unit_angle_rad,
  74                                                 &hkl_unit_angle_deg);
  75         }
  76
  77         return self;
  78 }
  79
  80 HklLattice *hkl_lattice_new_copy(HklLattice const *self)
  81 {
  82         HklLattice *copy = NULL;
  83
  84         copy = HKL_MALLOC(HklLattice);
  85
  86         copy->a = hkl_parameter_new_copy(self->a);
  87         copy->b = hkl_parameter_new_copy(self->b);
  88         copy->c = hkl_parameter_new_copy(self->c);
  89         copy->alpha = hkl_parameter_new_copy(self->alpha);
  90         copy->beta = hkl_parameter_new_copy(self->beta);
  91         copy->gamma = hkl_parameter_new_copy(self->gamma);
  92
  93         return copy;
  94 }
  95
  96 HklLattice* hkl_lattice_new_default(void)
  97 {
  98         return hkl_lattice_new(1.54, 1.54, 1.54,
  99                                90*HKL_DEGTORAD, 90*HKL_DEGTORAD, 90*HKL_DEGTORAD);
 100 }
 101
 102 void hkl_lattice_free(HklLattice *self)
 103 {
 104         hkl_parameter_free(self->a);
 105         hkl_parameter_free(self->b);
 106         hkl_parameter_free(self->c);
 107         hkl_parameter_free(self->alpha);
 108         hkl_parameter_free(self->beta);
 109         hkl_parameter_free(self->gamma);
 110         free(self);
 111 }
 112
 113 int hkl_lattice_set(HklLattice *self,
 114                     double a, double b, double c,
 115                     double alpha, double beta, double gamma)
 116 {
 117         int res = HKL_FAIL;
 118
 119         if(!check_lattice_param(a, b, c, alpha, beta, gamma)) {
 120                 self->a->value = a;
 121                 self->b->value = b;
 122                 self->c->value = c;
 123                 self->alpha->value = alpha;
 124                 self->beta->value = beta;
 125                 self->gamma->value = gamma;
 126                 res = HKL_SUCCESS;
 127         }
 128         return res;
 129 }
 130
 131 /*
 132  * Get the B matrix from the l parameters
 133  */
 134 int hkl_lattice_get_B(HklLattice const *self, HklMatrix *B)
 135 {
 136         double D;
 137         double c_alpha, s_alpha;
 138         double c_beta, s_beta;
 139         double c_gamma, s_gamma;
 140         double b11, b22, tmp;
 141
 142         c_alpha = cos(self->alpha->value);
 143         c_beta = cos(self->beta->value);
 144         c_gamma = cos(self->gamma->value);
 145         D = 1 - c_alpha*c_alpha - c_beta*c_beta - c_gamma*c_gamma
 146                 + 2*c_alpha*c_beta*c_gamma;
 147
 148         if (D > 0.)
 149                 D = sqrt(D);
 150         else
 151                 return HKL_FAIL;
 152
 153         s_alpha = sin(self->alpha->value);
 154         s_beta  = sin(self->beta->value);
 155         s_gamma = sin(self->gamma->value);
 156
 157         b11 = HKL_TAU / (self->b->value * s_alpha);
 158         b22 = HKL_TAU / self->c->value;
 159         tmp = b22 / s_alpha;
 160
 161         B->data[0][0] = HKL_TAU * s_alpha / (self->a->value * D);
 162         B->data[0][1] = b11 / D * (c_alpha*c_beta - c_gamma);
 163         B->data[0][2] = tmp / D * (c_gamma*c_alpha - c_beta);
 164
 165         B->data[1][0] = 0;
 166         B->data[1][1] = b11;
 167         B->data[1][2] = tmp / (s_beta*s_gamma) * (c_beta*c_gamma - c_alpha);
 168
 169         B->data[2][0] = 0;
 170         B->data[2][1] = 0;
 171         B->data[2][2] = b22;
 172
 173         return HKL_SUCCESS;
 174 }
 175
 176 /**
 177  * hkl_lattice_get_1_B:
 178  * @self: the @HklLattice
 179  * @B: the @HklMatrix returned
 180  *
 181  * Compute the invert of B (needed by the hkl_sample_set_UB method)
 182  * should be optimized
 183  *
 184  * Returns: HKL_SUCCESS or HKL_FAIL depending of the success of the
 185  * computation.
 186  **/
 187 int hkl_lattice_get_1_B(const HklLattice *self, HklMatrix *B)
 188 {
 189         HklMatrix tmp;
 190         double a;
 191         double b;
 192         double c;
 193         double d;
 194         double e;
 195         double f;
 196
 197         if(!self || !B)
 198                 return HKL_FAIL;
 199
 200         /*
 201          * first compute the B matrix
 202          * | a b c |
 203          * | 0 d e |
 204          * | 0 0 f |
 205          */
 206         hkl_lattice_get_B(self, &tmp);
 207
 208         /*
 209          * now invert this triangular matrix
 210          */
 211         a = tmp.data[0][0];
 212         b = tmp.data[0][1];
 213         c = tmp.data[0][2];
 214         d = tmp.data[1][1];
 215         e = tmp.data[1][2];
 216         f = tmp.data[2][2];
 217
 218         B->data[0][0] = 1 / a;
 219         B->data[0][1] = -b / a / d;
 220         B->data[0][2] = (b * e - d * c) / a / d / f;
 221
 222         B->data[1][0] = 0;
 223         B->data[1][1] = 1 / d;
 224         B->data[1][2] = -e / d / f;
 225
 226         B->data[2][0] = 0;
 227         B->data[2][1] = 0;
 228         B->data[2][2] = 1 / f;
 229
 230         return HKL_SUCCESS;
 231 }
 232
 233 int hkl_lattice_reciprocal(HklLattice const *self, HklLattice *reciprocal)
 234 {
 235         double c_alpha, c_beta, c_gamma;
 236         double s_alpha, s_beta, s_gamma;
 237         double c_beta1, c_beta2, c_beta3;
 238         double s_beta1, s_beta2, s_beta3;
 239         double s_beta_s_gamma, s_gamma_s_alpha, s_alpha_s_beta;
 240         double D;
 241
 242         c_alpha = cos(self->alpha->value);
 243         c_beta  = cos(self->beta->value);
 244         c_gamma = cos(self->gamma->value);
 245         D = 1 - c_alpha*c_alpha - c_beta*c_beta - c_gamma*c_gamma
 246                 + 2*c_alpha*c_beta*c_gamma;
 247
 248         if (D > 0.)
 249                 D = sqrt(D);
 250         else
 251                 return HKL_FAIL;
 252
 253         s_alpha = sin(self->alpha->value);
 254         s_beta  = sin(self->beta->value);
 255         s_gamma = sin(self->gamma->value);
 256
 257         s_beta_s_gamma  = s_beta  * s_gamma;
 258         s_gamma_s_alpha = s_gamma * s_alpha;
 259         s_alpha_s_beta  = s_alpha * s_beta;
 260
 261         c_beta1 = (c_beta  * c_gamma - c_alpha) / s_beta_s_gamma;
 262         c_beta2 = (c_gamma * c_alpha - c_beta)  / s_gamma_s_alpha;
 263         c_beta3 = (c_alpha * c_beta  - c_gamma) / s_alpha_s_beta;
 264         s_beta1 = D / s_beta_s_gamma;
 265         s_beta2 = D / s_gamma_s_alpha;
 266         s_beta3 = D / s_alpha_s_beta;
 267
 268         reciprocal->a->value = HKL_TAU * s_alpha / (self->a->value * D);
 269         reciprocal->b->value = HKL_TAU * s_beta  / (self->b->value * D);
 270         reciprocal->c->value = HKL_TAU * s_gamma / (self->c->value * D);
 271         reciprocal->alpha->value = atan2(s_beta1, c_beta1);
 272         reciprocal->beta->value  = atan2(s_beta2, c_beta2);
 273         reciprocal->gamma->value = atan2(s_beta3, c_beta3);
 274
 275         return HKL_SUCCESS;
 276 }
 277
 278 void hkl_lattice_randomize(HklLattice *self)
 279 {
 280         static HklVector vector_x = {{1, 0, 0}};
 281         HklVector a, b, c;
 282         HklVector axe;
 283         unsigned int angles_to_randomize;
 284
 285         /* La valeur des angles alpha, beta et gamma ne sont pas indépendant. */
 286         /* Il faut donc gérer les différents cas. */
 287         hkl_parameter_randomize(self->a);
 288         hkl_parameter_randomize(self->b);
 289         hkl_parameter_randomize(self->c);
 290
 291         angles_to_randomize = self->alpha->fit
 292                 + self->beta->fit
 293                 + self->gamma->fit;
 294         switch (angles_to_randomize) {
 295         case 0:
 296                 break;
 297         case 1:
 298                 if (self->alpha->fit) {
 299                         /* alpha */
 300                         a = b = c = vector_x;
 301
 302                         /* randomize b */
 303                         hkl_vector_randomize_vector(&axe, &a);
 304                         hkl_vector_rotated_around_vector(&b, &axe, self->gamma->value);
 305
 306                         /* randomize c */
 307                         hkl_vector_randomize_vector(&axe, &a);
 308                         hkl_vector_rotated_around_vector(&c, &axe, self->beta->value);
 309
 310                         /* compute the alpha angle. */
 311                         self->alpha->value = hkl_vector_angle(&b, &c);
 312                 } else if (self->beta->fit) {
 313                         /* beta */
 314                         a = b = vector_x;
 315
 316                         /* randomize b */
 317                         hkl_vector_randomize_vector(&axe, &a);
 318                         hkl_vector_rotated_around_vector(&b, &axe, self->gamma->value);
 319
 320                         /* randomize c */
 321                         c = b;
 322                         hkl_vector_randomize_vector(&axe, &b);
 323                         hkl_vector_rotated_around_vector(&c, &axe, self->alpha->value);
 324
 325                         /* compute beta */
 326                         self->beta->value = hkl_vector_angle(&a, &c);
 327                 } else {
 328                         /* gamma */
 329                         a = c = vector_x;
 330
 331                         /* randomize c */
 332                         hkl_vector_randomize_vector(&axe, &a);
 333                         hkl_vector_rotated_around_vector(&c, &axe, self->beta->value);
 334
 335                         /* randomize b */
 336                         b = c;
 337                         hkl_vector_randomize_vector(&axe, &c);
 338                         hkl_vector_rotated_around_vector(&b, &axe, self->alpha->value);
 339
 340                         /* compute gamma */
 341                         self->gamma->value = hkl_vector_angle(&a, &b);
 342                 }
 343                 break;
 344         case 2:
 345                 if (self->alpha->fit) {
 346                         if (self->beta->fit) {
 347                                 /* alpha + beta */
 348                                 a = b = vector_x;
 349
 350                                 /* randomize b */
 351                                 hkl_vector_randomize_vector(&axe, &a);
 352                                 hkl_vector_rotated_around_vector(&b, &axe, self->gamma->value);
 353
 354                                 /* randomize c */
 355                                 hkl_vector_randomize_vector_vector(&c, &a, &b);
 356
 357                                 self->alpha->value = hkl_vector_angle(&b, &c);
 358                                 self->beta->value = hkl_vector_angle(&a, &c);
 359                         } else {
 360                                 /* alpha + gamma */
 361                                 a = c = vector_x;
 362
 363                                 /* randomize c */
 364                                 hkl_vector_randomize_vector(&axe, &a);
 365                                 hkl_vector_rotated_around_vector(&c, &axe, self->beta->value);
 366
 367                                 /* randomize c */
 368                                 hkl_vector_randomize_vector_vector(&b, &a, &c);
 369
 370                                 self->alpha->value = hkl_vector_angle(&b, &c);
 371                                 self->gamma->value = hkl_vector_angle(&a, &b);
 372                         }
 373                 } else {
 374                         /* beta + gamma */
 375                         b = c = vector_x;
 376
 377                         /* randomize c */
 378                         hkl_vector_randomize_vector(&axe, &b);
 379                         hkl_vector_rotated_around_vector(&c, &axe, self->alpha->value);
 380
 381                         /* randomize c */
 382                         hkl_vector_randomize_vector_vector(&a, &b, &c);
 383
 384                         self->beta->value = hkl_vector_angle(&a, &c);
 385                         self->gamma->value = hkl_vector_angle(&a, &b);
 386                 }
 387                 break;
 388         case 3:
 389                 hkl_vector_randomize(&a);
 390                 hkl_vector_randomize_vector(&b, &a);
 391                 hkl_vector_randomize_vector_vector(&c, &b, &a);
 392
 393                 self->alpha->value = hkl_vector_angle(&b, &c);
 394                 self->beta->value = hkl_vector_angle(&a, &c);
 395                 self->gamma->value = hkl_vector_angle(&a, &b);
 396                 break;
 397         }
 398 }
 399
 400 void hkl_lattice_fprintf(FILE *f, HklLattice const *self)
 401 {
 402         fprintf(f, "\n");
 403         hkl_parameter_fprintf(f, self->a);
 404         fprintf(f, "\n");
 405         hkl_parameter_fprintf(f, self->b);
 406         fprintf(f, "\n");
 407         hkl_parameter_fprintf(f, self->c);
 408         fprintf(f, "\n");
 409         hkl_parameter_fprintf(f, self->alpha);
 410         fprintf(f, "\n");
 411         hkl_parameter_fprintf(f, self->beta);
 412         fprintf(f, "\n");
 413         hkl_parameter_fprintf(f, self->gamma);
 414 }