hkl/hkl-pseudoaxis-common-hkl.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-2016 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  *          Maria-Teresa Nunez-Pardo-de-Verra <tnunez@mail.desy.de>
  22  */
  23 #include <gsl/gsl_errno.h>              // for ::GSL_SUCCESS, etc
  24 #include <gsl/gsl_multiroots.h>
  25 #include <gsl/gsl_sf_trig.h>            // for gsl_sf_angle_restrict_pos
  26 #include <gsl/gsl_vector_double.h>      // for gsl_vector, etc
  27 #include <math.h>                       // for fabs, M_PI
  28 #include <stddef.h>                     // for size_t
  29 #include <stdlib.h>                     // for free, malloc, rand, etc
  30 #include <string.h>                     // for NULL
  31 #include <sys/types.h>                  // for uint
  32 #include "hkl-axis-private.h"           // for HklAxis
  33 #include "hkl-detector-private.h"       // for hkl_detector_compute_kf
  34 #include "hkl-geometry-private.h"       // for HklHolder, _HklGeometry, etc
  35 #include "hkl-macros-private.h"         // for hkl_assert, HKL_MALLOC, etc
  36 #include "hkl-matrix-private.h"         // for hkl_matrix_times_vector, etc
  37 #include "hkl-parameter-private.h"      // for _HklParameter, etc
  38 #include "hkl-pseudoaxis-auto-private.h"  // for CHECK_NAN, etc
  39 #include "hkl-pseudoaxis-common-hkl-private.h"  // for HklEngineHkl
  40 #include "hkl-pseudoaxis-common-q-private.h"  // for HklEngineHkl
  41 #include "hkl-pseudoaxis-private.h"     // for _HklEngine, _HklMode, etc
  42 #include "hkl-quaternion-private.h"     // for hkl_quaternion_init, etc
  43 #include "hkl-sample-private.h"         // for _HklSample
  44 #include "hkl-source-private.h"         // for hkl_source_compute_ki
  45 #include "hkl-vector-private.h"         // for HklVector, etc
  46 #include "hkl.h"                        // for HklEngine, HklGeometry, etc
  47 #include "hkl/ccan/array_size/array_size.h"  // for ARRAY_SIZE
  48 #include "hkl/ccan/container_of/container_of.h"  // for container_of
  49 #include "hkl/ccan/darray/darray.h"     // for darray_item, darray_size
  50
  51 /* #define DEBUG */
  52
  53 /*******************************************/
  54 /* common methode use by hkl getter/setter */
  55 /*******************************************/
  56
  57 typedef struct _HklDetectorFit HklDetectorFit;
  58
  59 struct _HklDetectorFit
  60 {
  61         HklGeometry *geometry;
  62         HklDetector *detector;
  63         HklVector *kf0;
  64         HklParameter **axes;
  65         size_t len;
  66 };
  67
  68 /* this method is used to fit only the detector position */
  69 /* usable with only 1 or 2 axes */
  70 static int fit_detector_function(const gsl_vector *x, void *params, gsl_vector *f)
  71 {
  72         size_t i;
  73         HklDetectorFit *fitp = params;
  74         HklVector kf;
  75
  76         /* update the workspace from x; */
  77         for(i=0; i<fitp->len; ++i)
  78                 hkl_parameter_value_set(fitp->axes[i],
  79                                         x->data[i],
  80                                         HKL_UNIT_DEFAULT, NULL);
  81
  82         hkl_geometry_update(fitp->geometry);
  83
  84         hkl_detector_compute_kf(fitp->detector, fitp->geometry, &kf);
  85
  86         f->data[0] = fabs(fitp->kf0->data[0] - kf.data[0])
  87                 + fabs(fitp->kf0->data[1] - kf.data[1])
  88                 + fabs(fitp->kf0->data[2] - kf.data[2]);
  89         if (fitp->len > 1)
  90                 f->data[1] = fabs(fitp->kf0->data[1] - kf.data[1]);
  91
  92 #if 0
  93         fprintf(stdout, "\nkf0 [%f, %f, %f], kf [%f, %f, %f]",
  94                 fitp->kf0->data[0], fitp->kf0->data[1], fitp->kf0->data[2],
  95                 kf.data[0], kf.data[1], kf.data[2]);
  96         fprintf(stdout, " x : [");
  97         for(i=0; i<fitp->len; ++i)
  98                 fprintf(stdout, " %.7f", x_data[i]);
  99         fprintf(stdout, "] |  f : [");
 100         for(i=0; i<fitp->len; ++i)
 101                 fprintf(stdout, " %.7f", f_data[i]);
 102         fprintf(stdout, "]\n");
 103 #endif
 104         return GSL_SUCCESS;
 105 }
 106
 107
 108 static int fit_detector_position(HklMode *mode, HklGeometry *geometry,
 109                                  HklDetector *detector, HklVector *kf)
 110 {
 111         const char **axis_name;
 112         HklDetectorFit params;
 113         gsl_multiroot_fsolver_type const *T;
 114         gsl_multiroot_fsolver *s;
 115         gsl_multiroot_function f;
 116         gsl_vector *x;
 117         int status;
 118         int res = FALSE;
 119         int iter;
 120         HklHolder *sample_holder = darray_item(geometry->holders, 0);
 121         HklHolder *detector_holder = darray_item(geometry->holders, 1);
 122
 123         /* fit the detector part to find the position of the detector for a given kf */
 124         /* FIXME for now the sample and detector holder are respectively the first and the second one */
 125         /* we need to find the right axes to use for the fit */
 126         /* BECARFULL the sample part must not move during this fit. So exclude an axis */
 127         /* if it is also part of the sample holder. */
 128         /* For now compare the holder axes with the axes of the mode to generate the right gsl multiroot solver */
 129         params.geometry = geometry;
 130         params.detector = detector;
 131         params.kf0 = kf;
 132         params.axes = malloc(sizeof(*params.axes) * detector_holder->config->len);
 133         params.len = 0;
 134         /* for each axis of the mode */
 135         darray_foreach(axis_name, mode->info->axes_w){
 136                 size_t k;
 137                 size_t tmp;
 138
 139                 tmp = hkl_geometry_get_axis_idx_by_name(params.geometry, *axis_name);
 140                 /* check that this axis is in the detector's holder */
 141                 for(k=0; k<detector_holder->config->len; ++k)
 142                         if(tmp == detector_holder->config->idx[k]){
 143                                 size_t j;
 144                                 int ko = 0;
 145
 146                                 /* and not in the sample's holder */
 147                                 for(j=0; j<sample_holder->config->len; ++j){
 148                                         if (tmp == sample_holder->config->idx[j]){
 149                                                 ko = 1;
 150                                                 break;
 151                                         }
 152                                 }
 153                                 if(!ko)
 154                                         params.axes[params.len++] = darray_item(params.geometry->axes, tmp);
 155                         }
 156         }
 157
 158         /* if no detector axis found ???? abort */
 159         /* maybe put this at the begining of the method */
 160         if (params.len > 0){
 161                 size_t i;
 162
 163                 /* now solve the system */
 164                 /* Initialize method  */
 165                 T = gsl_multiroot_fsolver_hybrid;
 166                 s = gsl_multiroot_fsolver_alloc (T, params.len);
 167                 x = gsl_vector_alloc(params.len);
 168
 169                 /* initialize x with the right values */
 170                 for(i=0; i<params.len; ++i)
 171                         x->data[i] = hkl_parameter_value_get(params.axes[i], HKL_UNIT_DEFAULT);
 172
 173                 f.f = fit_detector_function;
 174                 f.n = params.len;
 175                 f.params = &params;
 176                 gsl_multiroot_fsolver_set (s, &f, x);
 177
 178                 /* iterate to find the solution */
 179                 iter = 0;
 180                 do {
 181                         ++iter;
 182                         status = gsl_multiroot_fsolver_iterate(s);
 183                         if (status || iter % 100 == 0) {
 184                                 /* Restart from another point. */
 185                                 for(i=0; i<params.len; ++i)
 186                                         x->data[i] = (double)rand() / RAND_MAX * 180. / M_PI;
 187                                 gsl_multiroot_fsolver_set(s, &f, x);
 188                                 gsl_multiroot_fsolver_iterate(s);
 189                         }
 190                         status = gsl_multiroot_test_residual (s->f, HKL_EPSILON);
 191                 } while (status == GSL_CONTINUE && iter < 1000);
 192
 193 #ifdef DEBUG
 194                 fprintf(stdout, "\n  fitting the detector position using thoses axes :");
 195                 for(i=0; i<params.len; ++i)
 196                         fprintf(stdout, " \"%s\"", ((HklParameter *)params.axes[i])->name);
 197                 fprintf(stdout, " status : %d iter : %d", status, iter);
 198                 fprintf(stdout, " x: [");
 199                 for(i=0; i<params.len; ++i)
 200                         fprintf(stdout, " %.7f", s->x->data[i]);
 201                 fprintf(stdout, "] f: [");
 202                 for(i=0; i<params.len; ++i)
 203                         fprintf(stdout, " %.7f", s->f->data[i]);
 204                 fprintf(stdout, "]\n");
 205                 hkl_geometry_fprintf(stdout, params.geometry);
 206 #endif
 207                 if(status != GSL_CONTINUE){
 208                         res = TRUE;
 209                         /* put the axes in the -pi, pi range. */
 210                         for(i=0; i<params.len; ++i){
 211                                 double value;
 212
 213                                 value = hkl_parameter_value_get(params.axes[i], HKL_UNIT_DEFAULT);
 214                                 /* TODO one day deal with the error for real */
 215                                 hkl_parameter_value_set(params.axes[i],
 216                                                         gsl_sf_angle_restrict_pos(value),
 217                                                         HKL_UNIT_DEFAULT, NULL);
 218                         }
 219                 }
 220                 /* release memory */
 221                 gsl_vector_free(x);
 222                 gsl_multiroot_fsolver_free(s);
 223         }
 224         free(params.axes);
 225
 226         return res;
 227 }
 228
 229 /* get the highest index of the axis in a holder */
 230 /* BEWARE, NOT the axis index in the geometry->axes */
 231 /* which is part of the axis_names of the mode */
 232 /* return -1 if there is no axes of the mode in the sample part of the geometry */
 233 static int get_last_axis_idx(HklGeometry *geometry, int holder_idx, const darray_string *axes)
 234 {
 235         int last = -1;
 236         const char **axis_name;
 237         HklHolder *holder;
 238
 239         holder = darray_item(geometry->holders, holder_idx);
 240         darray_foreach(axis_name, *axes){
 241                 size_t i;
 242                 size_t idx;
 243
 244                 /* FIXME for now the sample holder is the first one */
 245                 idx = hkl_geometry_get_axis_idx_by_name(geometry, *axis_name);
 246                 for(i=0; i<holder->config->len; ++i)
 247                         if(idx == holder->config->idx[i]){
 248                                 last = last > (int)i ? last : (int)i;
 249                                 break;
 250                         }
 251         }
 252         return last;
 253 }
 254
 255
 256 static int hkl_is_reachable(HklEngine *engine, double wavelength, GError **error)
 257 {
 258         HklEngineHkl *engine_hkl = container_of(engine, HklEngineHkl, engine);
 259         HklVector Hkl = {
 260                 .data = {
 261                         engine_hkl->h->_value,
 262                         engine_hkl->k->_value,
 263                         engine_hkl->l->_value,
 264                 },
 265         };
 266
 267         hkl_matrix_times_vector(&engine->sample->UB, &Hkl);
 268         if (hkl_vector_norm2(&Hkl) > qmax(wavelength)){
 269                 g_set_error(error,
 270                             HKL_ENGINE_ERROR,
 271                             HKL_ENGINE_ERROR_SET,
 272                             "unreachable hkl, try to change the wavelength");
 273                 return FALSE;
 274         }
 275
 276         return TRUE;
 277 }
 278
 279 /**
 280  * _RUBh_minus_Q_func: (skip)
 281  * @x:
 282  * @params:
 283  * @f:
 284  *
 285  * Only usefull if you need to create a new hkl mode.
 286  *
 287  * Returns:
 288  **/
 289 int _RUBh_minus_Q_func(const gsl_vector *x, void *params, gsl_vector *f)
 290 {
 291         CHECK_NAN(x->data, x->size);
 292
 293         return RUBh_minus_Q(x->data, params, f->data);
 294 }
 295
 296 /**
 297  * RUBh_minus_Q: (skip)
 298  * @x:
 299  * @params:
 300  * @f:
 301  *
 302  *
 303  *
 304  * Returns:
 305  **/
 306 int RUBh_minus_Q(double const x[], void *params, double f[])
 307 {
 308         HklEngine *engine = params;
 309         HklEngineHkl *engine_hkl = container_of(engine, HklEngineHkl, engine);
 310         HklVector Hkl = {
 311                 .data = {
 312                         engine_hkl->h->_value,
 313                         engine_hkl->k->_value,
 314                         engine_hkl->l->_value,
 315                 },
 316         };
 317         HklVector ki, dQ;
 318         HklHolder *sample_holder;
 319
 320         /* update the workspace from x; */
 321         set_geometry_axes(engine, x);
 322
 323         /* R * UB * h = Q */
 324         /* for now the 0 holder is the sample holder. */
 325         sample_holder = darray_item(engine->geometry->holders, 0);
 326         hkl_matrix_times_vector(&engine->sample->UB, &Hkl);
 327         hkl_vector_rotated_quaternion(&Hkl, &sample_holder->q);
 328
 329         /* kf - ki = Q */
 330         hkl_source_compute_ki(&engine->geometry->source, &ki);
 331         hkl_detector_compute_kf(engine->detector, engine->geometry, &dQ);
 332         hkl_vector_minus_vector(&dQ, &ki);
 333
 334         hkl_vector_minus_vector(&dQ, &Hkl);
 335
 336         f[0] = dQ.data[0];
 337         f[1] = dQ.data[1];
 338         f[2] = dQ.data[2];
 339
 340         return GSL_SUCCESS;
 341 }
 342
 343 int hkl_mode_get_hkl_real(HklMode *self,
 344                           HklEngine *engine,
 345                           HklGeometry *geometry,
 346                           HklDetector *detector,
 347                           HklSample *sample,
 348                           GError **error)
 349 {
 350         HklHolder *sample_holder;
 351         HklMatrix RUB;
 352         HklVector hkl, ki, Q;
 353         HklEngineHkl *engine_hkl = container_of(engine, HklEngineHkl, engine);
 354
 355         /* update the geometry internals */
 356         hkl_geometry_update(geometry);
 357
 358         /* R * UB */
 359         /* for now the 0 holder is the sample holder. */
 360         sample_holder = darray_item(geometry->holders, 0);
 361         hkl_quaternion_to_matrix(&sample_holder->q, &RUB);
 362         hkl_matrix_times_matrix(&RUB, &sample->UB);
 363
 364         /* kf - ki = Q */
 365         hkl_source_compute_ki(&geometry->source, &ki);
 366         hkl_detector_compute_kf(detector, geometry, &Q);
 367         hkl_vector_minus_vector(&Q, &ki);
 368
 369         hkl_matrix_solve(&RUB, &hkl, &Q);
 370
 371         engine_hkl->h->_value = hkl.data[0];
 372         engine_hkl->k->_value = hkl.data[1];
 373         engine_hkl->l->_value = hkl.data[2];
 374
 375         return TRUE;
 376 }
 377
 378 int hkl_mode_set_hkl_real(HklMode *self,
 379                           HklEngine *engine,
 380                           HklGeometry *geometry,
 381                           HklDetector *detector,
 382                           HklSample *sample,
 383                           GError **error)
 384 {
 385         int last_axis;
 386
 387         hkl_error (error == NULL || *error == NULL);
 388
 389         /* check the input parameters */
 390         if(!hkl_is_reachable(engine, geometry->source.wave_length,
 391                              error)){
 392                 hkl_assert(error == NULL || *error != NULL);
 393                 return FALSE;
 394         }
 395         hkl_assert(error == NULL || *error == NULL);
 396
 397         /* compute the mode */
 398         if(!hkl_mode_auto_set_real(self, engine,
 399                                    geometry, detector, sample,
 400                                    error)){
 401                 hkl_assert(error == NULL || *error != NULL);
 402                 //fprintf(stdout, "message :%s\n", (*error)->message);
 403                 return FALSE;
 404         }
 405         hkl_assert(error == NULL || *error == NULL);
 406
 407         /* check that the mode allow to move a sample axis */
 408         /* FIXME for now the sample holder is the first one */
 409         last_axis = get_last_axis_idx(geometry, 0, &self->info->axes_w);
 410         if(last_axis >= 0){
 411                 uint i;
 412                 const HklGeometryListItem *item;
 413                 uint len = engine->engines->geometries->n_items;
 414
 415                 /* For each solution already found we will generate another one */
 416                 /* using the Ewalds construction by rotating Q around the last sample */
 417                 /* axis of the mode until it intersect again the Ewald sphere. */
 418                 /* FIXME do not work if ki is colinear with the axis. */
 419
 420                 /* for this we needs : */
 421                 /* - the coordinates of the end of the Q vector (q) */
 422                 /* - the last sample axis orientation of the mode (axis_v) */
 423                 /* - the coordinates of the center of the ewalds sphere (c) */
 424                 /* - the coordinates of the center of rotation of the sample (o = 0, 0, 0) */
 425
 426                 /* then we can : */
 427                 /* - project the origin in plane of normal axis_v containing q (o') */
 428                 /* - project the center of the ewalds sphere into the same plan (c') */
 429                 /* - rotate q around this (o', c') line of 180° to find the (q2) solution */
 430                 /* - compute the (kf2) corresponding to this q2 solution */
 431                 /* at the end we just need to solve numerically the position of the detector */
 432
 433                 /* we will add solution to the geometries so save its length before */
 434                 for(i=0, item=list_top(&engine->engines->geometries->items, HklGeometryListItem, list);
 435                     i<len;
 436                     ++i, item=list_next(&engine->engines->geometries->items, item, list)){
 437                         int j;
 438                         HklVector ki;
 439                         HklVector kf2;
 440                         HklVector q;
 441                         HklVector axis_v;
 442                         HklQuaternion qr;
 443                         HklAxis *axis;
 444                         HklVector cp = {{0}};
 445                         HklVector op = {{0}};
 446                         double angle;
 447                         HklGeometry *geom;
 448
 449                         geom = hkl_geometry_new_copy(item->geometry);
 450
 451                         /* get the Q vector kf - ki */
 452                         hkl_detector_compute_kf(detector, geom, &q);
 453                         hkl_source_compute_ki(&geom->source, &ki);
 454                         hkl_vector_minus_vector(&q, &ki);
 455
 456                         /* compute the current orientation of the last axis */
 457                         axis = container_of(darray_item(geom->axes,
 458                                                         darray_item(geom->holders, 0)->config->idx[last_axis]),
 459                                             HklAxis, parameter);
 460                         axis_v = axis->axis_v;
 461                         hkl_quaternion_init(&qr, 1, 0, 0, 0);
 462                         for(j=0; j<last_axis; ++j)
 463                                 hkl_quaternion_times_quaternion(
 464                                         &qr,
 465                                         &container_of(darray_item(geom->axes,
 466                                                                   darray_item(geom->holders, 0)->config->idx[j]),
 467                                                       HklAxis, parameter)->q);
 468                         hkl_vector_rotated_quaternion(&axis_v, &qr);
 469
 470                         /* - project the center of the ewalds sphere into the same plan (c') */
 471                         hkl_vector_minus_vector(&cp, &ki);
 472                         hkl_vector_project_on_plan_with_point(&cp, &axis_v, &q);
 473                         hkl_vector_project_on_plan_with_point(&op, &axis_v, &q);
 474
 475                         /* - rotate q around this (o', c') line of 180° to find the (q2) solution */
 476                         kf2 = q;
 477                         hkl_vector_rotated_around_line(&kf2, M_PI, &cp, &op);
 478                         angle = hkl_vector_oriented_angle_points(&q, &op, &kf2, &axis_v);
 479                         /* TODO parameter list for geometry */
 480                         if(!hkl_parameter_value_set(&axis->parameter,
 481                                                     hkl_parameter_value_get(&axis->parameter, HKL_UNIT_DEFAULT) + angle,
 482                                                     HKL_UNIT_DEFAULT, error))
 483                                 return FALSE;
 484                         hkl_geometry_update(geom);
 485 #ifdef DEBUG
 486                         fprintf(stdout, "\n- try to add a solution by rotating Q <%f, %f, %f> around the \"%s\" axis <%f, %f, %f> of %f radian",
 487                                 q.data[0], q.data[1], q.data[2],
 488                                 ((HklParameter *)axis)->name,
 489                                 axis_v.data[0], axis_v.data[1], axis_v.data[2],
 490                                 angle);
 491                         fprintf(stdout, "\n   op: <%f, %f, %f>", op.data[0], op.data[1], op.data[2]);
 492                         fprintf(stdout, "\n   q2: <%f, %f, %f>", kf2.data[0], kf2.data[1], kf2.data[2]);
 493 #endif
 494                         hkl_vector_add_vector(&kf2, &ki);
 495
 496                         /* at the end we just need to solve numerically the position of the detector */
 497                         if(fit_detector_position(self, geom, detector, &kf2))
 498                                 hkl_geometry_list_add(engine->engines->geometries,
 499                                                       geom);
 500
 501                         hkl_geometry_free(geom);
 502                 }
 503         }
 504         return TRUE;
 505 }
 506
 507 /***************************************/
 508 /* the double diffraction get set part */
 509 /***************************************/
 510
 511 /**
 512  * double_diffraction: (skip)
 513  * @x:
 514  * @params:
 515  * @f:
 516  *
 517  *
 518  *
 519  * Returns:
 520  **/
 521 int _double_diffraction(double const x[], void *params, double f[])
 522 {
 523         HklEngine *engine = params;
 524         HklEngineHkl *engine_hkl = container_of(engine, HklEngineHkl, engine);
 525         HklVector hkl = {
 526                 .data = {
 527                         engine_hkl->h->_value,
 528                         engine_hkl->k->_value,
 529                         engine_hkl->l->_value,
 530                 },
 531         };
 532         HklVector kf2;
 533         HklVector ki;
 534         HklVector dQ;
 535         HklHolder *sample_holder;
 536
 537         /* update the workspace from x; */
 538         set_geometry_axes(engine, x);
 539
 540         /* get the second hkl from the mode parameters */
 541         hkl_vector_init(&kf2,
 542                         darray_item(engine->mode->parameters, 0)->_value,
 543                         darray_item(engine->mode->parameters, 1)->_value,
 544                         darray_item(engine->mode->parameters, 2)->_value);
 545
 546         /* R * UB * hkl = Q */
 547         /* for now the 0 holder is the sample holder. */
 548         sample_holder = darray_item(engine->geometry->holders, 0);
 549         hkl_matrix_times_vector(&engine->sample->UB, &hkl);
 550         hkl_vector_rotated_quaternion(&hkl, &sample_holder->q);
 551
 552         /* kf - ki = Q */
 553         hkl_source_compute_ki(&engine->geometry->source, &ki);
 554         hkl_detector_compute_kf(engine->detector, engine->geometry, &dQ);
 555         hkl_vector_minus_vector(&dQ, &ki);
 556         hkl_vector_minus_vector(&dQ, &hkl);
 557
 558         /* R * UB * hlk2 = Q2 */
 559         hkl_matrix_times_vector(&engine->sample->UB, &kf2);
 560         hkl_vector_rotated_quaternion(&kf2, &sample_holder->q);
 561         hkl_vector_add_vector(&kf2, &ki);
 562
 563         f[0] = dQ.data[0];
 564         f[1] = dQ.data[1];
 565         f[2] = dQ.data[2];
 566         f[3] = hkl_vector_norm2(&kf2) - hkl_vector_norm2(&ki);
 567
 568         return GSL_SUCCESS;
 569 }
 570
 571 /**
 572  * double_diffraction_func: (skip)
 573  * @x:
 574  * @params:
 575  * @f:
 576  *
 577  *
 578  *
 579  * Returns:
 580  **/
 581 int _double_diffraction_func(gsl_vector const *x, void *params, gsl_vector *f)
 582 {
 583         CHECK_NAN(x->data, x->size);
 584
 585         _double_diffraction(x->data, params, f->data);
 586
 587         return  GSL_SUCCESS;
 588 }
 589
 590
 591 /******************************************/
 592 /* the psi_constant_vertical get set part */
 593 /******************************************/
 594
 595 /**
 596  * psi_constant_vertical_func: (skip)
 597  * @x:
 598  * @params:
 599  * @f:
 600  *
 601  *
 602  *
 603  * Returns:
 604  **/
 605 int _psi_constant_vertical_func(gsl_vector const *x, void *params, gsl_vector *f)
 606 {
 607         HklVector ki, kf, Q;
 608         HklEngine *engine = params;
 609
 610         CHECK_NAN(x->data, x->size);
 611
 612         RUBh_minus_Q(x->data, params, f->data);
 613
 614         /* update the workspace from x; */
 615         set_geometry_axes(engine, x->data);
 616
 617         /* kf - ki = Q */
 618         hkl_source_compute_ki(&engine->geometry->source, &ki);
 619         hkl_detector_compute_kf(engine->detector, engine->geometry, &kf);
 620         Q = kf;
 621         hkl_vector_minus_vector(&Q, &ki);
 622
 623         f->data[3] = darray_item(engine->mode->parameters, 3)->_value;
 624
 625         /* if |Q| > epsilon ok */
 626         if(hkl_vector_normalize(&Q)){
 627                 HklVector hkl;
 628                 HklVector n;
 629
 630                 /* compute n the intersection of the plan P(kf, ki) and PQ (normal Q) */
 631                 n = kf;
 632                 hkl_vector_vectorial_product(&n, &ki);
 633                 hkl_vector_vectorial_product(&n, &Q);
 634
 635                 /* compute the hkl ref position in the laboratory */
 636                 /* referentiel. The geometry was already updated. */
 637                 /* FIXME for now the 0 holder is the sample holder. */
 638                 hkl.data[0] = darray_item(engine->mode->parameters, 0)->_value;
 639                 hkl.data[1] = darray_item(engine->mode->parameters, 1)->_value;
 640                 hkl.data[2] = darray_item(engine->mode->parameters, 2)->_value;
 641                 hkl_matrix_times_vector(&engine->sample->UB, &hkl);
 642                 hkl_vector_rotated_quaternion(&hkl,
 643                                               &darray_item(engine->geometry->holders, 0)->q);
 644
 645                 /* project hkl on the plan of normal Q */
 646                 hkl_vector_project_on_plan(&hkl, &Q);
 647 #ifdef DEBUG
 648                 fprintf(stdout, "\n");
 649                 hkl_geometry_fprintf(stdout, engine->geometry);
 650                 fprintf(stdout, "\n");
 651                 fprintf(stdout, "%s n : <%f, %f, %f> hkl : <%f, %f, %f> Q : <%f, %f, %f> angle : %f\n",
 652                         __func__,
 653                         n.data[0], n.data[1], n.data[2],
 654                         hkl.data[0], hkl.data[1], hkl.data[2],
 655                         Q.data[0], Q.data[1], Q.data[2],
 656                         hkl_vector_oriented_angle(&n, &hkl, &Q) * HKL_RADTODEG);
 657 #endif
 658                 if(hkl_vector_norm2(&hkl) > HKL_EPSILON)
 659                         f->data[3] -=  hkl_vector_oriented_angle(&n, &hkl, &Q);
 660         }
 661
 662         return  GSL_SUCCESS;
 663 }
 664
 665 #define HKL_MODE_PSI_CONSTANT_VERTICAL_ERROR hkl_mode_psi_constant_vertical_error_quark ()
 666
 667 static GQuark hkl_mode_psi_constant_vertical_error_quark (void)
 668 {
 669         return g_quark_from_static_string ("hkl-mode-psi-constant-vertical-error-quark");
 670 }
 671
 672 typedef enum {
 673         HKL_MODE_PSI_CONSTANT_VERTICAL_ERROR_INITIALIZED_SET, /* can not init the engine */
 674 } HklModePsiConstantVerticalError;
 675
 676 int hkl_mode_initialized_set_psi_constant_vertical_real(HklMode *self,
 677                                                         HklEngine *engine,
 678                                                         HklGeometry *geometry,
 679                                                         HklDetector *detector,
 680                                                         HklSample *sample,
 681                                                         int initialized,
 682                                                         GError **error)
 683 {
 684         HklVector hkl;
 685         HklVector ki, kf, Q, n;
 686
 687         if(initialized){
 688                 /* kf - ki = Q */
 689                 hkl_source_compute_ki(&geometry->source, &ki);
 690                 hkl_detector_compute_kf(detector, geometry, &kf);
 691                 Q = kf;
 692                 hkl_vector_minus_vector(&Q, &ki);
 693
 694                 if (hkl_vector_is_null(&Q)){
 695                         g_set_error(error,
 696                                     HKL_MODE_PSI_CONSTANT_VERTICAL_ERROR,
 697                                     HKL_MODE_PSI_CONSTANT_VERTICAL_ERROR_INITIALIZED_SET,
 698                                     "can not initialize the \"%s\" mode with a null hkl (kf == ki)"
 699                                     "\nplease select a non-null hkl", self->info->name);
 700                         return FALSE;
 701                 }else{
 702                         /* needed for a problem of precision */
 703                         hkl_vector_normalize(&Q);
 704
 705                         /* compute the intersection of the plan P(kf, ki) and PQ (normal Q) */
 706                         n = kf;
 707                         hkl_vector_vectorial_product(&n, &ki);
 708                         hkl_vector_vectorial_product(&n, &Q);
 709
 710                         /* compute hkl in the laboratory referentiel */
 711                         /* the geometry was already updated in the detector compute kf */
 712                         /* for now the 0 holder is the sample holder */
 713                         hkl.data[0] = darray_item(self->parameters, 0)->_value;
 714                         hkl.data[1] = darray_item(self->parameters, 1)->_value;
 715                         hkl.data[2] = darray_item(self->parameters, 2)->_value;
 716                         hkl_matrix_times_vector(&sample->UB, &hkl);
 717                         hkl_vector_rotated_quaternion(&hkl,
 718                                                       &darray_item(geometry->holders, 0)->q);
 719
 720                         /* project hkl on the plan of normal Q */
 721                         hkl_vector_project_on_plan(&hkl, &Q);
 722
 723                         if (hkl_vector_is_null(&hkl)){
 724                                 g_set_error(error,
 725                                             HKL_MODE_PSI_CONSTANT_VERTICAL_ERROR,
 726                                             HKL_MODE_PSI_CONSTANT_VERTICAL_ERROR_INITIALIZED_SET,
 727                                             "can not initialize the \"%s\" mode"
 728                                             "\nwhen Q and the <h2, k2, l2> ref vector are colinear."
 729                                             "\nplease change one or both of them", engine->mode->info->name);
 730                                 return FALSE;
 731                         }else{
 732                                 /* compute the angle beetween hkl and n and
 733                                  * store in in the fourth parameter */
 734                                 if (!hkl_parameter_value_set(darray_item(self->parameters, 3),
 735                                                              hkl_vector_oriented_angle(&n, &hkl, &Q),
 736                                                              HKL_UNIT_DEFAULT, error))
 737                                         return FALSE;
 738                         }
 739                 }
 740         }
 741
 742         self->initialized = initialized;
 743
 744         return TRUE;
 745 }
 746
 747 /*******************/
 748 /* emergence fixed */
 749 /*******************/
 750
 751 typedef struct _HklModeAutoHklEmergenceFixed HklModeAutoHklEmergenceFixed;
 752
 753 struct _HklModeAutoHklEmergenceFixed
 754 {
 755         HklMode parent;
 756         HklParameter *n_x; /* not owned */
 757         HklParameter *n_y; /* not owned */
 758         HklParameter *n_z; /* not owned */
 759         HklParameter *emergence; /* not owned */
 760 };
 761
 762 #define HKL_MODE_HKL_EMERGENCE_FIXED_ERROR hkl_mode_hkl_emergence_fixed_error_quark ()
 763
 764 static GQuark hkl_mode_hkl_emergence_fixed_error_quark (void)
 765 {
 766         return g_quark_from_static_string ("hkl-mode-hkl-emergence-fixed-error-quark");
 767 }
 768
 769 typedef enum {
 770         HKL_MODE_HKL_EMERGENCE_FIXED_ERROR_INITIALIZED_SET, /* can not init the engine */
 771         HKL_MODE_HKL_EMERGENCE_FIXED_ERROR_SET, /* can not set the engine */
 772 } HklModeAutoHklEmergenceFixedError;
 773
 774
 775 static HklVector surface(const HklModeAutoHklEmergenceFixed *mode){
 776         HklVector n = {
 777                 .data = {
 778                         mode->n_x->_value,
 779                         mode->n_y->_value,
 780                         mode->n_z->_value,
 781                 }
 782         };
 783         return n;
 784 }
 785
 786 static double expected_emergence(const HklModeAutoHklEmergenceFixed *mode){
 787         return mode->emergence->_value;
 788 }
 789
 790 static int hkl_mode_hkl_emergence_fixed_initialized_set_real(HklMode *self,
 791                                                              HklEngine *engine,
 792                                                              HklGeometry *geometry,
 793                                                              HklDetector *detector,
 794                                                              HklSample *sample,
 795                                                              int initialized,
 796                                                              GError **error)
 797 {
 798         const HklModeAutoHklEmergenceFixed *mode = container_of(self, HklModeAutoHklEmergenceFixed, parent);
 799         HklVector kf;
 800         HklVector n = surface(mode);
 801
 802         /* first check the parameters */
 803         if (hkl_vector_is_null(&n)){
 804                 g_set_error(error,
 805                             HKL_MODE_HKL_EMERGENCE_FIXED_ERROR,
 806                             HKL_MODE_HKL_EMERGENCE_FIXED_ERROR_INITIALIZED_SET,
 807                             "Can not compute emergence fixed when the surface vector is null.");
 808                 return FALSE;
 809         }
 810
 811         /* compute the orientation of the surface */
 812         hkl_vector_rotated_quaternion(&n, &darray_item(geometry->holders, 0)->q);
 813
 814         hkl_detector_compute_kf(detector, geometry, &kf);
 815
 816         /* compute emergence and keep it */
 817         mode->emergence->_value = _emergence(&n, &kf);
 818
 819         self->initialized = initialized;
 820
 821         return TRUE;
 822 }
 823
 824
 825 int _emergence_fixed_func(const gsl_vector *x, void *params, gsl_vector *f)
 826 {
 827         HklEngine *engine = params;
 828         HklModeAutoHklEmergenceFixed *mode = container_of(engine->mode,
 829                                                           HklModeAutoHklEmergenceFixed,
 830                                                           parent);
 831         HklGeometry *geometry = engine->geometry;
 832         const HklDetector *detector = engine->detector;
 833         HklVector n = surface(mode);
 834         HklVector kf;
 835
 836         CHECK_NAN(x->data, x->size);
 837
 838         RUBh_minus_Q(x->data, params, f->data);
 839
 840         /* compute the orientation of the surface */
 841         hkl_vector_rotated_quaternion(&n, &darray_item(geometry->holders, 0)->q);
 842         hkl_detector_compute_kf(detector, geometry, &kf);
 843
 844         f->data[3] = expected_emergence(mode) - _emergence(&n, &kf);
 845
 846         return  GSL_SUCCESS;
 847 }
 848
 849 int hkl_mode_hkl_emergence_fixed_set_real(HklMode *self,
 850                                           HklEngine *engine,
 851                                           HklGeometry *geometry,
 852                                           HklDetector *detector,
 853                                           HklSample *sample,
 854                                           GError **error)
 855 {
 856         const HklModeAutoHklEmergenceFixed *mode = container_of(self, HklModeAutoHklEmergenceFixed, parent);
 857         HklVector n = surface(mode);
 858
 859         /* first check the parameters */
 860         if (hkl_vector_is_null(&n)){
 861                 g_set_error(error,
 862                             HKL_MODE_HKL_EMERGENCE_FIXED_ERROR,
 863                             HKL_MODE_HKL_EMERGENCE_FIXED_ERROR_SET,
 864                             "Can not compute hkl with emergence fixed when the surface vector is null.");
 865                 return FALSE;
 866         }
 867
 868         return hkl_mode_set_hkl_real(self, engine, geometry, detector, sample, error);
 869 }
 870
 871 HklMode *hkl_mode_hkl_emergence_fixed_new(const HklModeAutoInfo *auto_info)
 872 {
 873         static const HklModeOperations operations = {
 874                 HKL_MODE_OPERATIONS_HKL_FULL_DEFAULTS,
 875                 .capabilities = HKL_ENGINE_CAPABILITIES_READABLE | HKL_ENGINE_CAPABILITIES_WRITABLE | HKL_ENGINE_CAPABILITIES_INITIALIZABLE,
 876                 .initialized_set = hkl_mode_hkl_emergence_fixed_initialized_set_real,
 877                 .set = hkl_mode_hkl_emergence_fixed_set_real,
 878         };
 879         HklModeAutoHklEmergenceFixed *self;
 880
 881         if (darray_size(auto_info->info.axes_w) != 4){
 882                 fprintf(stderr, "This generic HklModeAutoHklEmergenceFixed need exactly 4 axes");
 883                 exit(128);
 884         }
 885
 886         self = HKL_MALLOC(HklModeAutoHklEmergenceFixed);
 887
 888         /* the base constructor; */
 889         hkl_mode_auto_init(&self->parent,
 890                            auto_info,
 891                            &operations, FALSE);
 892
 893         self->n_x = register_mode_parameter(&self->parent, 0);
 894         self->n_y = register_mode_parameter(&self->parent, 1);
 895         self->n_z = register_mode_parameter(&self->parent, 2);
 896         self->emergence = register_mode_parameter(&self->parent, 3);
 897
 898         return &self->parent;
 899 }
 900
 901 /*************/
 902 /* HklEngine */
 903 /*************/
 904
 905 static void hkl_engine_hkl_free_real(HklEngine *base)
 906 {
 907         HklEngineHkl *self = container_of(base, HklEngineHkl, engine);
 908         hkl_engine_release(&self->engine);
 909         free(self);
 910 }
 911
 912 HklEngine *hkl_engine_hkl_new(HklEngineList *engines)
 913 {
 914         HklEngineHkl *self;
 915         static const HklParameter h = {
 916                 HKL_PARAMETER_DEFAULTS, .name = "h",
 917                 .description = "h coordinate of the diffracting plan",
 918                 .range = { .min=-1, .max=1 },
 919         };
 920         static const HklParameter k = {
 921                 HKL_PARAMETER_DEFAULTS, .name = "k",
 922                 .description = "k coordinate of the diffracting plan",
 923                 .range = { .min=-1, .max=1 },
 924         };
 925         static const HklParameter l = {
 926                 HKL_PARAMETER_DEFAULTS, .name = "l",
 927                 .description = "l coordinate of the diffracting plan",
 928                 .range={ .min=-1, .max=1 },
 929         };
 930         static const HklParameter *pseudo_axes[] = {&h, &k, &l};
 931         static HklEngineInfo info = {
 932                 HKL_ENGINE_INFO("hkl",
 933                                 pseudo_axes,
 934                                 HKL_ENGINE_DEPENDENCIES_AXES | HKL_ENGINE_DEPENDENCIES_ENERGY | HKL_ENGINE_DEPENDENCIES_SAMPLE),
 935         };
 936         static HklEngineOperations operations = {
 937                 HKL_ENGINE_OPERATIONS_DEFAULTS,
 938                 .free=hkl_engine_hkl_free_real,
 939         };
 940
 941         self = HKL_MALLOC(HklEngineHkl);
 942
 943         hkl_engine_init(&self->engine, &info, &operations, engines);
 944
 945         self->h = register_pseudo_axis(&self->engine, engines, &h);
 946         self->k = register_pseudo_axis(&self->engine, engines, &k);
 947         self->l = register_pseudo_axis(&self->engine, engines, &l);
 948
 949         return &self->engine;
 950 }