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