src/lattice.cpp

   1 #include "config.h"
   2 #include "svector.h"
   3 #include "lattice.h"
   4
   5 extern struct hkl_svector hkl_svector_X;
   6
   7 namespace hkl
   8   {
   9
  10   /**
  11    * @brief The default constructor.
  12    */
  13   Lattice::Lattice()
  14   {
  15     _a = new FitParameter("a", "The a parameter of the crystal",
  16                           0., 1.54, 1000,
  17                           true, HKL_EPSILON);
  18     _b = new FitParameter("b", "The b parameter of the crystal",
  19                           0., 1.54, 1000,
  20                           true, HKL_EPSILON);
  21     _c = new FitParameter("c", "The c parameter of the crystal",
  22                           0., 1.54, 1000,
  23                           true, HKL_EPSILON);
  24     _alpha = new FitParameter("alpha", "The alpha parameter of the crystal",
  25                               0. * HKL_DEGTORAD, 90. * HKL_DEGTORAD, 180. * HKL_DEGTORAD,
  26                               true, HKL_EPSILON);
  27     _beta = new FitParameter("beta", "The beta parameter of the crystal",
  28                              0. * HKL_DEGTORAD, 90. * HKL_DEGTORAD, 180. * HKL_DEGTORAD,
  29                              true, HKL_EPSILON);
  30     _gamma = new FitParameter("gamma", "The gamma parameter of the cell",
  31                               0. * HKL_DEGTORAD, 90. * HKL_DEGTORAD, 180. * HKL_DEGTORAD,
  32                               true, HKL_EPSILON);
  33
  34     // set a old values different than current values to force _B computation.
  35     _old_a = 0;
  36     _old_b = 0;
  37     _old_c = 0;
  38     _old_alpha = 0;
  39     _old_beta = 0;
  40     _old_gamma = 0;
  41     _computeB();
  42   }
  43
  44   /**
  45    * @brief Another constructor.
  46    * @param a the a parameter of the Lattice
  47    * @param b the b parameter of the Lattice
  48    * @param c the c parameter of the Lattice
  49    * @param alpha the alpha parameter of the Lattice
  50    * @param beta the beta parameter of the Lattice
  51    * @param gamma the gamma parameter of the Lattice
  52    */
  53   Lattice::Lattice(const hkl::Value & a, const hkl::Value & b, const hkl::Value & c, const hkl::Value & alpha, const hkl::Value & beta, const hkl::Value & gamma)
  54   {
  55     _a = new FitParameter("a", "The a parameter of the crystal",
  56                           0., a, 1000,
  57                           true, HKL_EPSILON);
  58     _b = new FitParameter("b", "The b parameter of the crystal",
  59                           0., b, 1000,
  60                           true, HKL_EPSILON);
  61     _c = new FitParameter("c", "The c parameter of the crystal",
  62                           0., c, 1000,
  63                           true, HKL_EPSILON);
  64     _alpha = new FitParameter("alpha", "The alpha parameter of the crystal",
  65                               0. * HKL_DEGTORAD, alpha, 180. * HKL_DEGTORAD,
  66                               true, HKL_EPSILON);
  67     _beta = new FitParameter("beta", "The beta parameter of the crystal",
  68                              0. * HKL_DEGTORAD, beta, 180. * HKL_DEGTORAD,
  69                              true, HKL_EPSILON);
  70     _gamma = new FitParameter("gamma", "The gamma parameter of the cell",
  71                               0. * HKL_DEGTORAD, gamma, 180. * HKL_DEGTORAD,
  72                               true, HKL_EPSILON);
  73
  74     // set a old values different than current values to force _B computation.
  75     _old_a = 0;
  76     _old_b = 0;
  77     _old_c = 0;
  78     _old_alpha = 0;
  79     _old_beta = 0;
  80     _old_gamma = 0;
  81     _computeB();
  82   }
  83
  84   /**
  85    * @brief The copy constructor.
  86    * @param source The Lattice to copy.
  87    */
  88   Lattice::Lattice(const hkl::Lattice & source)
  89   {
  90     _a = new FitParameter(*(source._a));
  91     _b = new FitParameter(*(source._b));
  92     _c = new FitParameter(*(source._c));
  93     _alpha = new FitParameter(*(source._alpha));
  94     _beta = new FitParameter(*(source._beta));
  95     _gamma = new FitParameter(*(source._gamma));
  96
  97     // update the old value to compute the B matrix
  98     _old_a = source._old_a;
  99     _old_b = source._old_b;
 100     _old_c = source._old_c;
 101     _old_alpha = source._old_alpha;
 102     _old_beta = source._old_beta;
 103     _old_gamma = source._old_gamma;
 104
 105     _B = source._B;
 106   }
 107
 108   /**
 109    * @brief The default destructor.
 110    */
 111   Lattice::~Lattice()
 112   {
 113     delete _a;
 114     delete _b;
 115     delete _c;
 116     delete _alpha;
 117     delete _beta;
 118     delete _gamma;
 119   }
 120
 121   /**
 122    * @brief Get the a FitParameter of the Lattice.
 123    * @return A reference on the a FitParameter.
 124    * @todo return fitparameter * instead of fitParameter &.
 125    */
 126   hkl::FitParameter & Lattice::a()
 127   {
 128     return *_a;
 129   }
 130
 131   /**
 132    * @brief Get the b FitParameter of the Lattice.
 133    * @return A reference on the b FitParameter.
 134    * @todo return fitparameter * instead of fitParameter &.
 135    */
 136   hkl::FitParameter & Lattice::b()
 137   {
 138     return *_b;
 139   }
 140
 141   /**
 142    * @brief Get the c FitParameter of the Lattice.
 143    * @return A reference on the c FitParameter.
 144    * @todo return fitparameter * instead of fitParameter &.
 145    */
 146   hkl::FitParameter & Lattice::c()
 147   {
 148     return *_c;
 149   }
 150
 151   /**
 152    * @brief Get the alpha FitParameter of the Lattice.
 153    * @return A reference on the alpha FitParameter.
 154    * @todo return fitparameter * instead of fitParameter &.
 155    */
 156   hkl::FitParameter & Lattice::alpha()
 157   {
 158     return *_alpha;
 159   }
 160
 161   /**
 162    * @brief Get the beta FitParameter of the Lattice.
 163    * @return A reference on the beta FitParameter.
 164    * @todo return fitparameter * instead of fitParameter &.
 165    */
 166   hkl::FitParameter & Lattice::beta()
 167   {
 168     return *_beta;
 169   }
 170
 171   /**
 172    * @brief Get the gamma FitParameter of the Lattice.
 173    * @return A reference on the gamma FitParameter.
 174    * @todo return fitparameter * instead of fitParameter &.
 175    */
 176   hkl::FitParameter & Lattice::gamma()
 177   {
 178     return *_gamma;
 179   }
 180
 181   /**
 182    * @brief Get the a FitParameter of the Lattice.
 183    * @return A reference on the a FitParameter.
 184    * @todo return fitparameter * instead of fitParameter &.
 185    */
 186   const hkl::FitParameter & Lattice::a() const
 187     {
 188       return *_a;
 189     }
 190
 191   /**
 192    * @brief Get the b FitParameter of the Lattice.
 193    * @return A reference on the b FitParameter.
 194    * @todo return fitparameter * instead of fitParameter &.
 195    */
 196   const hkl::FitParameter & Lattice::b() const
 197     {
 198       return *_b;
 199     }
 200
 201   /**
 202    * @brief Get the c FitParameter of the Lattice.
 203    * @return A reference on the c FitParameter.
 204    * @todo return fitparameter * instead of fitParameter &.
 205    */
 206   const hkl::FitParameter & Lattice::c() const
 207     {
 208       return *_c;
 209     }
 210
 211   /**
 212    * @brief Get the alpha FitParameter of the Lattice.
 213    * @return A reference on the alpha FitParameter.
 214    * @todo return fitparameter * instead of fitParameter &.
 215    */
 216   const hkl::FitParameter & Lattice::alpha() const
 217     {
 218       return *_alpha;
 219     }
 220
 221   /**
 222    * @brief Get the beta FitParameter of the Lattice.
 223    * @return A reference on the beta FitParameter.
 224    * @todo return fitparameter * instead of fitParameter &.
 225    */
 226   const hkl::FitParameter & Lattice::beta() const
 227     {
 228       return *_beta;
 229     }
 230
 231   /**
 232    * @brief Get the gamma FitParameter of the Lattice.
 233    * @return A reference on the gamma FitParameter.
 234    * @todo return fitparameter * instead of fitParameter &.
 235    */
 236   const hkl::FitParameter & Lattice::gamma() const
 237     {
 238       return *_gamma;
 239     }
 240
 241   hkl_smatrix const * Lattice::get_B() const throw(hkl::HKLException)
 242   {
 243     bool status = _computeB();
 244     if (status)
 245       return &_B;
 246     else
 247       HKLEXCEPTION("can not compute B", "Check the lattice parameters");
 248   }
 249
 250   hkl_smatrix const * Lattice::get_B(bool & status) const
 251     {
 252       status = _computeB();
 253       return &_B;
 254     }
 255
 256   /**
 257    * @brief Compute the reciprocal Lattice.
 258    * @return The reciprocal Lattice.
 259    * @throw HKLException if the reciprocal Lattice can not be compute.
 260    * @todo See for the consign assignation.
 261    */
 262   hkl::Lattice Lattice::reciprocal() const throw(hkl::HKLException)
 263   {
 264     double a_star, b_star, c_star;
 265     double alpha_star, beta_star, gamma_star;
 266
 267     _compute_reciprocal(a_star, b_star, c_star, alpha_star, beta_star, gamma_star);
 268
 269     hkl::Lattice lattice(a_star, b_star, c_star, alpha_star, beta_star, gamma_star);
 270     lattice._a->set_consign(_a->get_consign());
 271     lattice._b->set_consign(_b->get_consign());
 272     lattice._c->set_consign(_c->get_consign());
 273     lattice._alpha->set_consign(_alpha->get_consign());
 274     lattice._beta->set_consign(_beta->get_consign());
 275     lattice._gamma->set_consign(_gamma->get_consign());
 276
 277     return lattice;
 278   }
 279
 280   /**
 281    * @brief Randomize the Lattice.
 282    */
 283   void Lattice::randomize()
 284   {
 285     static hkl_svector svector_x = {{1, 0, 0}};
 286     hkl_svector a, b, c;
 287     hkl_svector axe;
 288
 289     // La valeur des angles alpha, beta et gamma ne sont pas indépendant.
 290     // Il faut donc gérer les différents cas.
 291
 292     _a->randomize();
 293     _b->randomize();
 294     _c->randomize();
 295     unsigned int angles_to_randomize = _alpha->get_flagFit() + _beta->get_flagFit() + _gamma->get_flagFit();
 296
 297     switch (angles_to_randomize)
 298       {
 299       case 0:
 300         break;
 301       case 1:
 302         if (_alpha->get_flagFit()) // alpha
 303           {
 304             a = b = c = svector_x;
 305
 306             // randomize b
 307             ::hkl_svector_randomize_svector(&axe, &a);
 308             ::hkl_svector_rotated_around_vector(&b, &axe, _gamma->get_current().get_value());
 309
 310             // randomize c
 311             ::hkl_svector_randomize_svector(&axe, &a);
 312             ::hkl_svector_rotated_around_vector(&c, &axe, _beta->get_current().get_value());
 313
 314             //compute the alpha angle.
 315             _alpha->set_current(::hkl_svector_angle(&b, &c));
 316           }
 317         else if (_beta->get_flagFit())
 318           {
 319             // beta
 320             a = b = svector_x;
 321
 322             // randomize b
 323             ::hkl_svector_randomize_svector(&axe, &a);
 324             ::hkl_svector_rotated_around_vector(&b, &axe, _gamma->get_current().get_value());
 325
 326             // randomize c
 327             c = b;
 328             ::hkl_svector_randomize_svector(&axe, &b);
 329             ::hkl_svector_rotated_around_vector(&c, &axe, _alpha->get_current().get_value());
 330
 331             //compute beta
 332             _beta->set_current(::hkl_svector_angle(&a, &c));
 333           }
 334         else
 335           {
 336             // gamma
 337             a = c = svector_x;
 338
 339             // randomize c
 340             ::hkl_svector_randomize_svector(&axe, &a);
 341             ::hkl_svector_rotated_around_vector(&c, &axe, _beta->get_current().get_value());
 342
 343             // randomize b
 344             b = c;
 345             ::hkl_svector_randomize_svector(&axe, &c);
 346             ::hkl_svector_rotated_around_vector(&b, &axe, _alpha->get_current().get_value());
 347
 348             //compute beta
 349             _gamma->set_current(::hkl_svector_angle(&a, &b));
 350           }
 351         break;
 352       case 2:
 353         if (_alpha->get_flagFit())
 354           {
 355             if (_beta->get_flagFit()) // alpha + beta
 356               {
 357                 a = b = svector_x;
 358
 359                 // randomize b
 360                 ::hkl_svector_randomize_svector(&axe, &a);
 361                 ::hkl_svector_rotated_around_vector(&b, &axe, _gamma->get_current().get_value());
 362
 363                 //randomize c
 364                 ::hkl_svector_randomize_svector_svector(&c, &a, &b);
 365
 366                 _alpha->set_current(::hkl_svector_angle(&b, &c));
 367                 _beta->set_current(::hkl_svector_angle(&a, &c));
 368               }
 369             else
 370               {
 371                 // alpha + gamma
 372                 a = c = svector_x;
 373
 374                 // randomize c
 375                 ::hkl_svector_randomize_svector(&axe, &a);
 376                 ::hkl_svector_rotated_around_vector(&c, &axe, _beta->get_current().get_value());
 377
 378                 //randomize c
 379                 ::hkl_svector_randomize_svector_svector(&b, &a, &c);
 380
 381                 _alpha->set_current(::hkl_svector_angle(&b, &c));
 382                 _gamma->set_current(::hkl_svector_angle(&a, &b));
 383               }
 384           }
 385         else
 386           {
 387             // beta + gamma
 388             b = c = svector_x;
 389
 390             // randomize c
 391             ::hkl_svector_randomize_svector(&axe, &b);
 392             ::hkl_svector_rotated_around_vector(&c, &axe, _alpha->get_current().get_value());
 393
 394             //randomize c
 395             ::hkl_svector_randomize_svector_svector(&a, &b, &c);
 396
 397             _beta->set_current(::hkl_svector_angle(&a, &c));
 398             _gamma->set_current(::hkl_svector_angle(&a, &b));
 399           }
 400         break;
 401       case 3:
 402         ::hkl_svector_randomize(&a);
 403         ::hkl_svector_randomize_svector(&b, &a);
 404         ::hkl_svector_randomize_svector_svector(&c, &b, &a);
 405
 406         _alpha->set_current(::hkl_svector_angle(&b, &c));
 407         _beta->set_current(::hkl_svector_angle(&a, &c));
 408         _gamma->set_current(::hkl_svector_angle(&a, &b));
 409         break;
 410       }
 411     // no exception the lattice is always valid.
 412     _computeB();
 413   }
 414
 415   /**
 416    * \brief Are two Lattice equals ?
 417    * \param lattice the hkl::Lattice to compare with.
 418    * \return true if both are equals flase otherwise.
 419    */
 420   bool Lattice::operator==(const hkl::Lattice & lattice) const
 421     {
 422       return *_a == *(lattice._a)
 423              && *_b == *(lattice._b)
 424              && *_c == *(lattice._c)
 425              && *_alpha == *(lattice._alpha)
 426              && *_beta == *(lattice._beta)
 427              && *_gamma == *(lattice._gamma);
 428     }
 429
 430   /**
 431    * @brief print the Lattice into a flux
 432    * @param flux The stream to print into.
 433    * @return The modified flux.
 434    */
 435   std::ostream & Lattice::printToStream(std::ostream & flux) const
 436     {
 437       flux << *_a << std::endl
 438       << *_b << std::endl
 439       << *_c << std::endl
 440       << *_alpha << std::endl
 441       << *_beta << std::endl
 442       << *_gamma << std::endl;
 443
 444       return flux;
 445     }
 446
 447   /**
 448    * @brief compute the B matrix from the fitParameters.
 449    * @return true if the calculus is valid.
 450    */
 451   bool Lattice::_computeB() const
 452     {
 453       double a = _a->get_current().get_value();
 454       double b = _b->get_current().get_value();
 455       double c = _c->get_current().get_value();
 456       double alpha = _alpha->get_current().get_value();
 457       double beta = _beta->get_current().get_value();
 458       double gamma = _gamma->get_current().get_value();
 459
 460       if ((a != _old_a)
 461           ||(b != _old_b)
 462           ||(c != _old_c)
 463           ||(alpha != _old_alpha)
 464           ||(beta != _old_beta)
 465           ||(gamma != _old_gamma))
 466         {
 467           double cos_alpha = cos(alpha);
 468           double cos_beta = cos(beta);
 469           double cos_gamma = cos(gamma);
 470           double D = 1 - cos_alpha*cos_alpha - cos_beta*cos_beta - cos_gamma*cos_gamma + 2*cos_alpha*cos_beta*cos_gamma;
 471
 472           if (D > 0.)
 473             D = sqrt(D);
 474           else
 475             return false;
 476
 477           double sin_alpha = sin(alpha);
 478           double sin_beta = sin(beta);
 479           double sin_gamma = sin(gamma);
 480
 481           // optimization (18*, 3+)
 482           double a_star = HKL_TAU * sin_alpha / (a * D);
 483
 484           double b_star_sin_gamma_star = HKL_TAU / (b * sin_alpha);
 485           double b_star_cos_gamma_star = b_star_sin_gamma_star / D * (cos_alpha*cos_beta - cos_gamma);
 486
 487           double tmp = HKL_TAU / (c * sin_alpha);
 488           double c_star_cos_beta_star = tmp / D * (cos_gamma*cos_alpha - cos_beta);
 489           double c_star_sin_beta_star_cos_alpha_star = tmp / (sin_beta * sin_gamma) * (cos_beta*cos_gamma - cos_alpha);
 490           // end of optimization
 491
 492           _B.data[0][0] = a_star;
 493           _B.data[0][1] = b_star_cos_gamma_star;
 494           _B.data[0][2] = c_star_cos_beta_star;
 495
 496           _B.data[1][0] = 0;
 497           _B.data[1][1] = b_star_sin_gamma_star;
 498           _B.data[1][2] = c_star_sin_beta_star_cos_alpha_star;
 499
 500           _B.data[2][0] = 0;
 501           _B.data[2][1] = 0;
 502           _B.data[2][2] = HKL_TAU / c;
 503
 504           _old_a = a;
 505           _old_b = b;
 506           _old_c = c;
 507           _old_alpha = alpha;
 508           _old_beta = beta;
 509           _old_gamma = gamma;
 510         }
 511       return true;
 512     }
 513
 514   /**
 515    * @brief compute the reciprocal parameters of the Lattice.
 516    * @param[out] a_star the a_star value.
 517    * @param[out] b_star the b_star value.
 518    * @param[out] c_star the c_star value.
 519    * @param[out] alpha_star the alpha_star value.
 520    * @param[out] beta_star the beta_star value.
 521    * @param[out] gamma_star the gamma_star value.
 522    * @throw HKLException if the reciprocal calculus is not possible.
 523    */
 524   void Lattice::_compute_reciprocal(double & a_star, double & b_star, double & c_star, double & alpha_star, double & beta_star, double & gamma_star) const throw(hkl::HKLException)
 525   {
 526     double a = _a->get_current().get_value();
 527     double b = _b->get_current().get_value();
 528     double c = _c->get_current().get_value();
 529     double alpha = _alpha->get_current().get_value();
 530     double beta = _beta->get_current().get_value();
 531     double gamma = _gamma->get_current().get_value();
 532
 533
 534     // D = 1 - cos(alpha)*cos(alpha) - cos(beta)*cos(beta) - cos(gamma)*cos(gamma) + 2*cos(alpha)*cos(beta)*cos(gamma);
 535     double cos_alpha = cos(alpha);
 536     double cos_beta = cos(beta);
 537     double cos_gamma = cos(gamma);
 538     double D = 1 - cos_alpha*cos_alpha - cos_beta*cos_beta - cos_gamma*cos_gamma + 2*cos_alpha*cos_beta*cos_gamma;
 539
 540     if (D > 0.)
 541       D = sqrt(D);
 542     else
 543       HKLEXCEPTION("Incorrect lattice parameters", "Please check lattice parameters");
 544
 545     //
 546     double sin_alpha = sin(alpha);
 547     double sin_beta = sin(beta);
 548     double sin_gamma = sin(gamma);
 549
 550     double sin_beta_sin_gamma = sin_beta*sin_gamma;
 551     double sin_gamma_sin_alpha = sin_gamma*sin_alpha;
 552     double sin_alpha_sin_beta = sin_alpha*sin_beta;
 553
 554     double cos_beta1 = (cos_beta*cos_gamma - cos_alpha) / sin_beta_sin_gamma;
 555     double cos_beta2 = (cos_gamma*cos_alpha - cos_beta) / sin_gamma_sin_alpha;
 556     double cos_beta3 = (cos_alpha*cos_beta - cos_gamma) / sin_alpha_sin_beta;
 557     double sin_beta1 = D / sin_beta_sin_gamma;
 558     double sin_beta2 = D / sin_gamma_sin_alpha;
 559     double sin_beta3 = D / sin_alpha_sin_beta;
 560
 561     a_star = HKL_TAU * sin_alpha / (a * D);
 562     b_star = HKL_TAU * sin_beta / (b * D);
 563     c_star = HKL_TAU * sin_gamma / (c * D);
 564
 565     alpha_star = atan2(sin_beta1, cos_beta1);
 566     beta_star = atan2(sin_beta2, cos_beta2);
 567     gamma_star = atan2(sin_beta3, cos_beta3);
 568   }
 569
 570
 571 } // namespace hkl