hkl/hkl-engine-soleil-sixs-med.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-2015 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 <gsl/gsl_multiroots.h>
  23 #include "hkl-factory-private.h"        // for autodata_factories_, etc
  24 #include "hkl-axis-private.h"           // for HklAxis
  25 #include "hkl-pseudoaxis-common-hkl-private.h"  // for hkl_engine_hkl_new, etc
  26 #include "hkl-pseudoaxis-common-q-private.h"  // for hkl_engine_q2_new, etc
  27 #include "hkl-pseudoaxis-common-tth-private.h"  // for hkl_engine_tth2_new, etc
  28 #include "hkl-pseudoaxis-common-readonly-private.h"
  29
  30 #define PITCH "pitch"
  31 #define BETA "beta"
  32 #define MU "mu"
  33 #define OMEGA "omega"
  34 #define GAMMA "gamma"
  35 #define DELTA "delta"
  36 #define ETA_A "eta_a"
  37
  38 /* #define DEBUG */
  39
  40 /*********************/
  41 /* MED 2+2 HklEngine */
  42 /*********************/
  43
  44 static int _reflectivity_func(const gsl_vector *x, void *params, gsl_vector *f)
  45 {
  46         const double mu = x->data[0];
  47         const double gamma = x->data[2];
  48
  49         CHECK_NAN(x->data, x->size);
  50
  51         RUBh_minus_Q(x->data, params, f->data);
  52         f->data[3] = gamma - 2 * mu;
  53
  54         return  GSL_SUCCESS;
  55 }
  56
  57 static const HklFunction reflectivity_func = {
  58         .function = _reflectivity_func,
  59         .size = 4,
  60 };
  61
  62 static HklMode* mu_fixed_2_2()
  63 {
  64         static const char* axes_r[] = {BETA, MU, OMEGA, GAMMA, DELTA};
  65         static const char* axes_w[] = {OMEGA, GAMMA, DELTA};
  66         static const HklFunction *functions[] = {&RUBh_minus_Q_func};
  67         static const HklModeAutoInfo info = {
  68                 HKL_MODE_AUTO_INFO("mu_fixed", axes_r, axes_w, functions),
  69         };
  70
  71         return hkl_mode_auto_new(&info,
  72                                  &hkl_full_mode_operations,
  73                                  TRUE);
  74 }
  75
  76 static HklMode* reflectivity_2_2()
  77 {
  78         static const char* axes_r[] = {BETA, MU, OMEGA, GAMMA, DELTA};
  79         static const char* axes_w[] = {MU, OMEGA, GAMMA, DELTA};
  80         static const HklFunction *functions[] = {&reflectivity_func};
  81         static const HklModeAutoInfo info = {
  82                 HKL_MODE_AUTO_INFO("reflectivity", axes_r, axes_w, functions),
  83         };
  84
  85         return hkl_mode_auto_new(&info,
  86                                  &hkl_full_mode_operations,
  87                                  TRUE);
  88 }
  89
  90 static HklEngine *hkl_engine_soleil_sixs_med_2_2_hkl_new(HklEngineList *engines)
  91 {
  92         HklEngine *self;
  93         HklMode *default_mode;
  94
  95         self = hkl_engine_hkl_new(engines);
  96
  97         default_mode = mu_fixed_2_2();
  98         hkl_engine_add_mode(self, default_mode);
  99         hkl_engine_mode_set(self, default_mode);
 100
 101         hkl_engine_add_mode(self, reflectivity_2_2());
 102
 103         return self;
 104 }
 105
 106 /* mode incidence */
 107
 108 REGISTER_READONLY_INCIDENCE(hkl_engine_soleil_sixs_med_2_2_incidence_new,
 109                             P99_PROTECT({BETA, MU, OMEGA}));
 110
 111 /*********************/
 112 /* MED 1+2 HklEngine */
 113 /*********************/
 114
 115 static HklMode* pitch_fixed()
 116 {
 117         static const char *axes_r[] = {PITCH, MU, GAMMA, DELTA};
 118         static const char* axes_w[] = {MU, GAMMA, DELTA};
 119         static const HklFunction *functions[] = {&RUBh_minus_Q_func};
 120         static const HklModeAutoInfo info = {
 121                 HKL_MODE_AUTO_INFO(__func__, axes_r, axes_w, functions),
 122         };
 123
 124         return hkl_mode_auto_new(&info,
 125                                  &hkl_full_mode_operations,
 126                                  TRUE);
 127 }
 128
 129 static HklEngine *hkl_engine_soleil_sixs_med_1_2_hkl_new(HklEngineList *engines)
 130 {
 131         HklEngine *self;
 132         HklMode *default_mode;
 133
 134         self = hkl_engine_hkl_new(engines);
 135
 136         default_mode = pitch_fixed();
 137         hkl_engine_add_mode(self, default_mode);
 138         hkl_engine_mode_set(self, default_mode);
 139
 140         return self;
 141 }
 142
 143 /* mode incidence */
 144
 145 REGISTER_READONLY_INCIDENCE(hkl_engine_soleil_sixs_med_1_2_incidence_new,
 146                             P99_PROTECT({PITCH, MU}));
 147
 148 /*********************/
 149 /* MED 2+3 HklEngine */
 150 /*********************/
 151
 152 typedef struct _HklSlitsFit HklSlitsFit;
 153 struct _HklSlitsFit
 154 {
 155         HklGeometry *geometry;
 156         HklVector surface;
 157         unsigned int slits_id;
 158         unsigned int len;
 159         HklParameter *axis;
 160 };
 161
 162 static int slits_func(const gsl_vector *x, void *params, gsl_vector *f)
 163 {
 164         double const *x_data = gsl_vector_const_ptr(x, 0);
 165         double *f_data = gsl_vector_ptr(f, 0);
 166         HklVector n_slits = {{0, 0, 1}};
 167         HklSlitsFit *parameters = params;
 168
 169         hkl_parameter_value_set(parameters->axis, x_data[0], HKL_UNIT_DEFAULT, NULL);
 170         hkl_geometry_update(parameters->geometry);
 171
 172         /* compute the orientation of the slits */
 173         hkl_vector_rotated_quaternion(&n_slits,
 174                                       &darray_item(parameters->geometry->holders, 1)->q);
 175
 176         /* both directions must be perpendicular */
 177         f_data[0] = hkl_vector_scalar_product(&parameters->surface, &n_slits);
 178
 179         return  GSL_SUCCESS;
 180 }
 181
 182 static int fit_slits_orientation(HklSlitsFit *params)
 183 {
 184         size_t i;
 185         gsl_multiroot_fsolver_type const *T;
 186         gsl_multiroot_fsolver *s;
 187         gsl_multiroot_function f;
 188         gsl_vector *x;
 189         double *x_data;
 190         int status;
 191         int res = FALSE;
 192         int iter;
 193
 194         /* now solve the system */
 195         /* Initialize method  */
 196         T = gsl_multiroot_fsolver_hybrid;
 197         s = gsl_multiroot_fsolver_alloc (T, params->len);
 198         x = gsl_vector_alloc(params->len);
 199         x_data = gsl_vector_ptr(x, 0);
 200
 201         /* initialize x with the right values */
 202         x_data[0] = params->axis->_value;
 203
 204         f.f = slits_func;
 205         f.n = params->len;
 206         f.params = params;
 207         gsl_multiroot_fsolver_set (s, &f, x);
 208
 209         /* iterate to find the solution */
 210         iter = 0;
 211         do {
 212                 ++iter;
 213                 status = gsl_multiroot_fsolver_iterate(s);
 214                 if (status || iter % 100 == 0) {
 215                         /* Restart from another point. */
 216                         for(i=0; i<params->len; ++i)
 217                                 x_data[i] = (double)rand() / RAND_MAX * 180. / M_PI;
 218                         gsl_multiroot_fsolver_set(s, &f, x);
 219                         gsl_multiroot_fsolver_iterate(s);
 220                 }
 221                 status = gsl_multiroot_test_residual (s->f, HKL_EPSILON);
 222         } while (status == GSL_CONTINUE && iter < 1000);
 223
 224 #ifdef DEBUG
 225         fprintf(stdout, "\n  fitting the detector position using thoses axes :");
 226         for(i=0; i<params->len; ++i)
 227                 fprintf(stdout, " \"%s\"", params->axis->name);
 228         fprintf(stdout, " status : %d iter : %d", status, iter);
 229         fprintf(stdout, " x: [");
 230         for(i=0; i<params->len; ++i)
 231                 fprintf(stdout, " %.7f", s->x->data[i]);
 232         fprintf(stdout, "] f: [");
 233         for(i=0; i<params->len; ++i)
 234                 fprintf(stdout, " %.7f", s->f->data[i]);
 235         fprintf(stdout, "]\n");
 236         hkl_geometry_fprintf(stdout, params->geometry);
 237 #endif
 238         if(status != GSL_CONTINUE){
 239                 res = TRUE;
 240                 /* put the axes in the -pi, pi range. */
 241                 gsl_sf_angle_restrict_pos_e(&params->axis->_value);
 242         }
 243         /* release memory */
 244         gsl_vector_free(x);
 245         gsl_multiroot_fsolver_free(s);
 246
 247         return res;
 248 }
 249
 250 static void hkl_geometry_list_multiply_soleil_sixs_med_2_3(HklGeometryList *self,
 251                                                            HklGeometryListItem *item)
 252 {
 253         HklSlitsFit params;
 254         HklGeometry *geometry;
 255         double slits_position;
 256         HklHolder *sample_holder;
 257         HklHolder *detector_holder;
 258
 259         /* For each solution already found we will generate another one */
 260         /* we will set the right slit orientation for a given detector arm position */
 261         geometry = item->geometry;
 262         sample_holder = darray_item(geometry->holders, 0);
 263         detector_holder = darray_item(geometry->holders, 1);
 264
 265         /* get the index of the axis corresponding to the slits */
 266         /* for now the last holder is the detector one */
 267         params.slits_id = detector_holder->config->idx[detector_holder->config->len-1];
 268         params.len = 1; /* only one axis to fit */
 269         params.geometry = geometry;
 270         params.axis = darray_item(params.geometry->axes, params.slits_id);
 271
 272         /* compute the surface orientation fixed during the fit */
 273         /* use the last sample axis as sample surface normal */
 274         params.surface = container_of(darray_item(geometry->axes,
 275                                                   sample_holder->config->idx[sample_holder->config->len - 1]),
 276                                       HklAxis, parameter)->axis_v;
 277         hkl_vector_rotated_quaternion(&params.surface,
 278                                       &sample_holder->q);
 279
 280
 281         /* we just need to fit the slits orientation */
 282         /* save it's value before */
 283         slits_position = hkl_parameter_value_get(params.axis, HKL_UNIT_DEFAULT);
 284         if (fit_slits_orientation(&params) != TRUE)
 285                 hkl_parameter_value_set(params.axis, slits_position, HKL_UNIT_DEFAULT, NULL);
 286 }
 287
 288 static HklMode* mu_fixed_2_3()
 289 {
 290         static const char *axes_r[] = {BETA, MU, OMEGA, GAMMA, DELTA, ETA_A};
 291         static const char* axes_w[] = {OMEGA, GAMMA, DELTA};
 292         static const HklFunction *functions[] = {&RUBh_minus_Q_func};
 293         static const HklModeAutoInfo info = {
 294                 HKL_MODE_AUTO_INFO("mu_fixed", axes_r, axes_w, functions),
 295         };
 296
 297         return hkl_mode_auto_new(&info,
 298                                  &hkl_full_mode_operations,
 299                                  TRUE);
 300 }
 301
 302 static HklEngine *hkl_engine_soleil_sixs_med_2_3_hkl_new(HklEngineList *engines)
 303 {
 304         HklEngine *self;
 305         HklMode *default_mode;
 306
 307         self = hkl_engine_hkl_new(engines);
 308
 309         default_mode = mu_fixed_2_3();
 310         hkl_engine_add_mode(self, default_mode);
 311         hkl_engine_mode_set(self, default_mode);
 312
 313         return self;
 314 }
 315
 316 /***********************/
 317 /* SOLEIL SIXS MED 2+2 */
 318 /***********************/
 319
 320 #define HKL_GEOMETRY_TYPE_SOLEIL_SIXS_MED_2_2_DESCRIPTION               \
 321         "+ xrays source fix allong the :math:`\\vec{x}` direction (1, 0, 0)\n" \
 322         "+ 3 axes for the sample\n"                                     \
 323         "\n"                                                            \
 324         "  + **" BETA "** : rotation around the :math:`-\\vec{y}` direction (0, -1, 0)\n" \
 325         "  + **" MU "** : rotation around the :math:`\\vec{z}` direction (0, 0, 1)\n" \
 326         "  + **" OMEGA "** : rotating around the :math:`-\\vec{y}` direction (0, -1, 0)\n" \
 327         "\n"                                                            \
 328         "+ 3 axis for the detector\n"                                   \
 329         "\n"                                                            \
 330         "  + **" BETA "** : rotation around the :math:`-\\vec{y}` direction (0, -1, 0)\n" \
 331         "  + **" GAMMA "** : rotation around the :math:`\\vec{z}` direction (0, 0, 1)\n" \
 332         "  + **" DELTA "** : rotation around the :math:`-\\vec{y}` direction (0, -1, 0)\n"
 333
 334 static const char* hkl_geometry_soleil_sixs_med_2_2_axes[] = {BETA, MU, OMEGA, GAMMA, DELTA};
 335
 336 static HklGeometry *hkl_geometry_new_soleil_sixs_med_2_2(const HklFactory *factory)
 337 {
 338         HklGeometry *self = hkl_geometry_new(factory);
 339         HklHolder *h;
 340
 341         h = hkl_geometry_add_holder(self);
 342         hkl_holder_add_rotation_axis(h, BETA, 0, -1, 0);
 343         hkl_holder_add_rotation_axis(h, MU, 0, 0, 1);
 344         hkl_holder_add_rotation_axis(h, OMEGA, 0, -1, 0);
 345
 346         h = hkl_geometry_add_holder(self);
 347         hkl_holder_add_rotation_axis(h, BETA, 0, -1, 0);
 348         hkl_holder_add_rotation_axis(h, GAMMA, 0, 0, 1);
 349         hkl_holder_add_rotation_axis(h, DELTA, 0, -1, 0);
 350
 351         return self;
 352 }
 353
 354 static HklEngineList *hkl_engine_list_new_soleil_sixs_med_2_2(const HklFactory *factory)
 355 {
 356         HklEngineList *self = hkl_engine_list_new();
 357
 358         hkl_engine_soleil_sixs_med_2_2_hkl_new(self);
 359         hkl_engine_q2_new(self);
 360         hkl_engine_qper_qpar_new(self);
 361         hkl_engine_tth2_new(self);
 362         hkl_engine_soleil_sixs_med_2_2_incidence_new(self);
 363
 364         return self;
 365 }
 366
 367 REGISTER_DIFFRACTOMETER(soleil_sixs_med_2_2,"SOLEIL SIXS MED2+2", HKL_GEOMETRY_TYPE_SOLEIL_SIXS_MED_2_2_DESCRIPTION);
 368
 369 /***********************/
 370 /* SOLEIL SIXS MED 1+2 */
 371 /***********************/
 372
 373 #define HKL_GEOMETRY_TYPE_SOLEIL_SIXS_MED_1_2_DESCRIPTION               \
 374         "+ xrays source fix allong the :math:`\\vec{x}` direction (1, 0, 0)\n" \
 375         "+ 2 axes for the sample\n"                                     \
 376         "\n"                                                            \
 377         "  + **" PITCH "** : rotation around the :math:`-\\vec{y}` direction (0, -1, 0)\n" \
 378         "  + **" MU "** : rotation around the :math:`\\vec{z}` direction (0, 0, 1)\n" \
 379         "\n"                                                            \
 380         "+ 3 axis for the detector\n"                                   \
 381         "\n"                                                            \
 382         "  + **" PITCH "** : rotation around the :math:`-\\vec{y}` direction (0, -1, 0)\n" \
 383         "  + **" GAMMA "** : rotation around the :math:`\\vec{z}` direction (0, 0, 1)\n" \
 384         "  + **" DELTA "** : rotation around the :math:`-\\vec{y}` direction (0, -1, 0)\n"
 385
 386 static const char* hkl_geometry_soleil_sixs_med_1_2_axes[] = {PITCH, MU, GAMMA, DELTA};
 387
 388 static HklGeometry *hkl_geometry_new_soleil_sixs_med_1_2(const HklFactory *factory)
 389 {
 390         HklGeometry *self = hkl_geometry_new(factory);
 391         HklHolder *h;
 392
 393         h = hkl_geometry_add_holder(self);
 394         hkl_holder_add_rotation_axis(h, PITCH, 0, -1, 0);
 395         hkl_holder_add_rotation_axis(h, MU, 0, 0, 1);
 396
 397         h = hkl_geometry_add_holder(self);
 398         hkl_holder_add_rotation_axis(h, PITCH, 0, -1, 0);
 399         hkl_holder_add_rotation_axis(h, GAMMA, 0, 0, 1);
 400         hkl_holder_add_rotation_axis(h, DELTA, 0, -1, 0);
 401
 402         return self;
 403 }
 404
 405 static HklEngineList *hkl_engine_list_new_soleil_sixs_med_1_2(const HklFactory *factory)
 406 {
 407         HklEngineList *self = hkl_engine_list_new();
 408
 409         hkl_engine_soleil_sixs_med_1_2_hkl_new(self);
 410         hkl_engine_q2_new(self);
 411         hkl_engine_qper_qpar_new(self);
 412         hkl_engine_tth2_new(self);
 413         hkl_engine_soleil_sixs_med_1_2_incidence_new(self);
 414
 415         return self;
 416 }
 417
 418 REGISTER_DIFFRACTOMETER(soleil_sixs_med_1_2, "SOLEIL SIXS MED1+2", HKL_GEOMETRY_TYPE_SOLEIL_SIXS_MED_1_2_DESCRIPTION);
 419
 420
 421 /***********************/
 422 /* SOLEIL SIXS MED 2+3 */
 423 /***********************/
 424
 425 #define HKL_GEOMETRY_TYPE_SOLEIL_SIXS_MED_2_3_DESCRIPTION               \
 426         "+ xrays source fix allong the :math:`\\vec{x}` direction (1, 0, 0)\n" \
 427         "+ 3 axes for the sample\n"                                     \
 428         "\n"                                                            \
 429         "  + **" BETA "** : rotation around the :math:`-\\vec{y}` direction (0, -1, 0)\n" \
 430         "  + **" MU "** : rotation around the :math:`\\vec{z}` direction (0, 0, 1)\n" \
 431         "  + **" OMEGA "** : rotating around the :math:`-\\vec{y}` direction (0, -1, 0)\n" \
 432         "\n"                                                            \
 433         "+ 4 axis for the detector\n"                                   \
 434         "\n"                                                            \
 435         "  + **" BETA "** : rotation around the :math:`-\\vec{y}` direction (0, -1, 0)\n" \
 436         "  + **" GAMMA "** : rotation around the :math:`\\vec{z}` direction (0, 0, 1)\n" \
 437         "  + **" DELTA "** : rotation around the :math:`-\\vec{y}` direction (0, -1, 0)\n" \
 438         "  + **" ETA_A "** : rotation around the :math:`-\\vec{x}` direction (-1, 0, 0)\n"
 439
 440 static const char* hkl_geometry_soleil_sixs_med_2_3_axes[] = {BETA, MU, OMEGA, GAMMA, DELTA, ETA_A};
 441
 442 static HklGeometry *hkl_geometry_new_soleil_sixs_med_2_3(const HklFactory *factory)
 443 {
 444         HklGeometry *self = hkl_geometry_new(factory);
 445         HklHolder *h;
 446
 447         h = hkl_geometry_add_holder(self);
 448         hkl_holder_add_rotation_axis(h, BETA, 0, -1, 0);
 449         hkl_holder_add_rotation_axis(h, MU, 0, 0, 1);
 450         hkl_holder_add_rotation_axis(h, OMEGA, 0, -1, 0);
 451
 452         h = hkl_geometry_add_holder(self);
 453         hkl_holder_add_rotation_axis(h, BETA, 0, -1, 0);
 454         hkl_holder_add_rotation_axis(h, GAMMA, 0, 0, 1);
 455         hkl_holder_add_rotation_axis(h, DELTA, 0, -1, 0);
 456         hkl_holder_add_rotation_axis(h, ETA_A, -1, 0, 0);
 457
 458         return self;
 459 }
 460
 461 static HklEngineList *hkl_engine_list_new_soleil_sixs_med_2_3(const HklFactory *factory)
 462 {
 463         HklEngineList *self = hkl_engine_list_new();
 464
 465         self->geometries->multiply = hkl_geometry_list_multiply_soleil_sixs_med_2_3;
 466         hkl_engine_soleil_sixs_med_2_3_hkl_new(self);
 467         hkl_engine_q2_new(self);
 468         hkl_engine_qper_qpar_new(self);
 469         hkl_engine_tth2_new(self);
 470         hkl_engine_soleil_sixs_med_2_2_incidence_new(self);
 471
 472         return self;
 473 }
 474
 475 REGISTER_DIFFRACTOMETER(soleil_sixs_med_2_3, "SOLEIL SIXS MED2+3", HKL_GEOMETRY_TYPE_SOLEIL_SIXS_MED_2_3_DESCRIPTION);