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