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-2019, 2022 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-private.h"
  40 #include "hkl-pseudoaxis-common-hkl-private.h"  // for HklEngineHkl
  41 #include "hkl-pseudoaxis-common-q-private.h"  // for HklEngineHkl
  42 #include "hkl-pseudoaxis-private.h"     // for _HklEngine, _HklMode, etc
  43 #include "hkl-quaternion-private.h"     // for hkl_quaternion_init, etc
  44 #include "hkl-sample-private.h"         // for _HklSample
  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         kf = hkl_geometry_kf_get(fitp->geometry, fitp->detector);
  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 static const HklVector reciprocal_plan(const HklEngineHkl *engine_hkl)
 259 {
 260         const HklVector hkl = {
 261                 .data = {
 262                         engine_hkl->h->_value,
 263                         engine_hkl->k->_value,
 264                         engine_hkl->l->_value,
 265                 },
 266         };
 267         return hkl;
 268 }
 269
 270
 271 int hkl_is_reachable(HklEngine *engine, double wavelength, GError **error)
 272 {
 273         const HklEngineHkl *engine_hkl = container_of(engine, HklEngineHkl, engine);
 274         HklVector Hkl = reciprocal_plan(engine_hkl);
 275
 276         hkl_matrix_times_vector(&engine->sample->UB, &Hkl);
 277         if (hkl_vector_norm2(&Hkl) > qmax(wavelength)){
 278                 g_set_error(error,
 279                             HKL_ENGINE_ERROR,
 280                             HKL_ENGINE_ERROR_SET,
 281                             "unreachable hkl, try to change the wavelength");
 282                 return FALSE;
 283         }
 284
 285         return TRUE;
 286 }
 287
 288 /**
 289  * _RUBh_minus_Q_func: (skip)
 290  * @x:
 291  * @params:
 292  * @f:
 293  *
 294  * Only usefull if you need to create a new hkl mode.
 295  *
 296  * Returns:
 297  **/
 298 int _RUBh_minus_Q_func(const gsl_vector *x, void *params, gsl_vector *f)
 299 {
 300         CHECK_NAN(x->data, x->size);
 301
 302         return RUBh_minus_Q(x->data, params, f->data);
 303 }
 304
 305 /**
 306  * RUBh_minus_Q: (skip)
 307  * @x:
 308  * @params:
 309  * @f:
 310  *
 311  *
 312  *
 313  * Returns:
 314  **/
 315 int RUBh_minus_Q(double const x[], void *params, double f[])
 316 {
 317         HklEngine *engine = params;
 318         HklEngineHkl *engine_hkl = container_of(engine, HklEngineHkl, engine);
 319         const HklVector Hkl = reciprocal_plan(engine_hkl);
 320
 321         /* update the workspace from x; */
 322         set_geometry_axes(engine, x);
 323
 324         const struct HklHklWrite w = hkl_hkl_write(engine->geometry,
 325                                                    engine->detector,
 326                                                    engine->sample,
 327                                                    &Hkl);
 328         f[0] = w.dQ.data[0];
 329         f[1] = w.dQ.data[1];
 330         f[2] = w.dQ.data[2];
 331
 332         return GSL_SUCCESS;
 333 }
 334
 335 int hkl_mode_get_hkl_real(HklMode *self,
 336                           HklEngine *engine,
 337                           HklGeometry *geometry,
 338                           HklDetector *detector,
 339                           HklSample *sample,
 340                           GError **error)
 341 {
 342         HklEngineHkl *engine_hkl = container_of(engine, HklEngineHkl, engine);
 343
 344         /* update the geometry internals */
 345         hkl_geometry_update(geometry);
 346
 347         const struct HklHklRead r = hkl_hkl_read(geometry, detector, sample);
 348
 349         engine_hkl->h->_value = r.hkl.data[0];
 350         engine_hkl->k->_value = r.hkl.data[1];
 351         engine_hkl->l->_value = r.hkl.data[2];
 352
 353         return TRUE;
 354 }
 355
 356 int hkl_mode_set_hkl_real(HklMode *self,
 357                           HklEngine *engine,
 358                           HklGeometry *geometry,
 359                           HklDetector *detector,
 360                           HklSample *sample,
 361                           GError **error)
 362 {
 363         int last_axis;
 364
 365         hkl_error (error == NULL || *error == NULL);
 366
 367         /* check the input parameters */
 368         if(!hkl_is_reachable(engine, geometry->source.wave_length,
 369                              error)){
 370                 hkl_assert(error == NULL || *error != NULL);
 371                 return FALSE;
 372         }
 373         hkl_assert(error == NULL || *error == NULL);
 374
 375         /* compute the mode */
 376         if(!hkl_mode_auto_set_real(self, engine,
 377                                    geometry, detector, sample,
 378                                    error)){
 379                 hkl_assert(error == NULL || *error != NULL);
 380                 //fprintf(stdout, "message :%s\n", (*error)->message);
 381                 return FALSE;
 382         }
 383         hkl_assert(error == NULL || *error == NULL);
 384
 385         /* check that the mode allow to move a sample axis */
 386         /* FIXME for now the sample holder is the first one */
 387         last_axis = get_last_sample_axis_idx(geometry, sample, &self->info->axes_w);
 388         if(last_axis >= 0){
 389                 uint i;
 390                 const HklGeometryListItem *item;
 391                 uint len = engine->engines->geometries->n_items;
 392
 393                 /* For each solution already found we will generate another one */
 394                 /* using the Ewalds construction by rotating Q around the last sample */
 395                 /* axis of the mode until it intersect again the Ewald sphere. */
 396                 /* FIXME do not work if ki is colinear with the axis. */
 397
 398                 /* for this we needs : */
 399                 /* - the coordinates of the end of the Q vector (q) */
 400                 /* - the last sample axis orientation of the mode (axis_v) */
 401                 /* - the coordinates of the center of the ewalds sphere (c) */
 402                 /* - the coordinates of the center of rotation of the sample (o = 0, 0, 0) */
 403
 404                 /* then we can : */
 405                 /* - project the origin in plane of normal axis_v containing q (o') */
 406                 /* - project the center of the ewalds sphere into the same plan (c') */
 407                 /* - rotate q around this (o', c') line of 180° to find the (q2) solution */
 408                 /* - compute the (kf2) corresponding to this q2 solution */
 409                 /* at the end we just need to solve numerically the position of the detector */
 410
 411                 /* we will add solution to the geometries so save its length before */
 412                 for(i=0, item=list_top(&engine->engines->geometries->items, HklGeometryListItem, list);
 413                     i<len && NULL != item;
 414                     ++i, item=list_next(&engine->engines->geometries->items, item, list)){
 415                         int j;
 416                         HklVector ki;
 417                         HklVector kf2;
 418                         HklVector q;
 419                         HklVector axis_v;
 420                         HklQuaternion qr;
 421                         HklAxis *axis;
 422                         HklVector cp = {{0}};
 423                         HklVector op = {{0}};
 424                         double angle;
 425                         HklGeometry *geom = hkl_geometry_new_copy(item->geometry);
 426                         HklHolder *sample_holder = hkl_geometry_sample_holder_get(geom, sample);
 427
 428                         /* get the Q vector kf - ki */
 429                         ki = hkl_geometry_ki_get(geom);
 430                         q = hkl_geometry_kf_get(geom, detector);
 431                         hkl_vector_minus_vector(&q, &ki);
 432
 433                         /* compute the current orientation of the last axis */
 434                         axis = container_of(darray_item(geom->axes,
 435                                                         sample_holder->config->idx[last_axis]),
 436                                             HklAxis, parameter);
 437                         axis_v = axis->axis_v;
 438                         hkl_quaternion_init(&qr, 1, 0, 0, 0);
 439                         for(j=0; j<last_axis; ++j)
 440                                 hkl_quaternion_times_quaternion(
 441                                         &qr,
 442                                         &container_of(darray_item(geom->axes,
 443                                                                   sample_holder->config->idx[j]),
 444                                                       HklAxis, parameter)->q);
 445                         hkl_vector_rotated_quaternion(&axis_v, &qr);
 446
 447                         /* - project the center of the ewalds sphere into the same plan (c') */
 448                         hkl_vector_minus_vector(&cp, &ki);
 449                         hkl_vector_project_on_plan_with_point(&cp, &axis_v, &q);
 450                         hkl_vector_project_on_plan_with_point(&op, &axis_v, &q);
 451
 452                         /* - rotate q around this (o', c') line of 180° to find the (q2) solution */
 453                         kf2 = q;
 454                         hkl_vector_rotated_around_line(&kf2, M_PI, &cp, &op);
 455                         angle = hkl_vector_oriented_angle_points(&q, &op, &kf2, &axis_v);
 456                         /* TODO parameter list for geometry */
 457                         if(!hkl_parameter_value_set(&axis->parameter,
 458                                                     hkl_parameter_value_get(&axis->parameter, HKL_UNIT_DEFAULT) + angle,
 459                                                     HKL_UNIT_DEFAULT, error)){
 460                                 hkl_geometry_free(geom);
 461                                 return FALSE;
 462                         }
 463                         hkl_geometry_update(geom);
 464 #ifdef DEBUG
 465                         fprintf(stdout, "\n- try to add a solution by rotating Q <%f, %f, %f> around the \"%s\" axis <%f, %f, %f> of %f radian",
 466                                 q.data[0], q.data[1], q.data[2],
 467                                 ((HklParameter *)axis)->name,
 468                                 axis_v.data[0], axis_v.data[1], axis_v.data[2],
 469                                 angle);
 470                         fprintf(stdout, "\n   op: <%f, %f, %f>", op.data[0], op.data[1], op.data[2]);
 471                         fprintf(stdout, "\n   q2: <%f, %f, %f>", kf2.data[0], kf2.data[1], kf2.data[2]);
 472 #endif
 473                         hkl_vector_add_vector(&kf2, &ki);
 474
 475                         /* at the end we just need to solve numerically the position of the detector */
 476                         if(fit_detector_position(self, geom, detector, sample, &kf2))
 477                                 hkl_geometry_list_add(engine->engines->geometries,
 478                                                       geom);
 479
 480                         hkl_geometry_free(geom);
 481                 }
 482         }
 483         return TRUE;
 484 }
 485
 486 /***************************************/
 487 /* the double diffraction get set part */
 488 /***************************************/
 489
 490 static HklVector reciprocal_plan2(const HklMode *mode)
 491 {
 492         const HklVector hkl2 = {
 493                 .data = {
 494                         darray_item(mode->parameters, 0)->_value,
 495                         darray_item(mode->parameters, 1)->_value,
 496                         darray_item(mode->parameters, 2)->_value,
 497                 },
 498         };
 499         return hkl2;
 500 }
 501
 502 /**
 503  * double_diffraction: (skip)
 504  * @x:
 505  * @params:
 506  * @f:
 507  *
 508  *
 509  *
 510  * Returns:
 511  **/
 512 int _double_diffraction(double const x[], void *params, double f[])
 513 {
 514         HklEngine *engine = params;
 515         HklEngineHkl *engine_hkl = container_of(engine, HklEngineHkl, engine);
 516         const HklVector hkl1 = reciprocal_plan(engine_hkl);
 517         const HklVector hkl2 = reciprocal_plan2(engine->mode);
 518
 519         /* update the workspace from x; */
 520         set_geometry_axes(engine, x);
 521
 522         const struct HklDoubleDiffractionWrite w = hkl_double_diffraction_write(engine->geometry,
 523                                                                                 engine->detector,
 524                                                                                 engine->sample,
 525                                                                                 &hkl1, &hkl2);
 526         f[0] = w.hkl1W.dQ.data[0];
 527         f[1] = w.hkl1W.dQ.data[1];
 528         f[2] = w.hkl1W.dQ.data[2];
 529         f[3] = hkl_vector_norm2(&w.kf2) - hkl_vector_norm2(&w.hkl2W.ki); /* 2nd diffraction condition */
 530
 531         return GSL_SUCCESS;
 532 }
 533
 534 /**
 535  * double_diffraction_func: (skip)
 536  * @x:
 537  * @params:
 538  * @f:
 539  *
 540  *
 541  *
 542  * Returns:
 543  **/
 544 int _double_diffraction_func(gsl_vector const *x, void *params, gsl_vector *f)
 545 {
 546         CHECK_NAN(x->data, x->size);
 547
 548         _double_diffraction(x->data, params, f->data);
 549
 550         return  GSL_SUCCESS;
 551 }
 552
 553 /******************************************/
 554 /* the psi_constant_vertical get set part */
 555 /******************************************/
 556
 557 /**
 558  * psi_constant_vertical_func: (skip)
 559  * @x:
 560  * @params:
 561  * @f:
 562  *
 563  *
 564  *
 565  * Returns:
 566  **/
 567 int _psi_constant_vertical_func(gsl_vector const *x, void *params, gsl_vector *f)
 568 {
 569         CHECK_NAN(x->data, x->size);
 570
 571         HklEngine *engine = params;
 572         HklEngineHkl *engine_hkl = container_of(engine, HklEngineHkl, engine);
 573         const HklVector hkl = reciprocal_plan(engine_hkl);
 574         const HklVector hkl2 = reciprocal_plan2(engine->mode);
 575
 576         /* update the workspace from x; */
 577         set_geometry_axes(engine, x->data);
 578
 579         const struct HklPsiWrite w = hkl_psi_write(engine->geometry,
 580                                                    engine->detector,
 581                                                    engine->sample,
 582                                                    &hkl, &hkl2);
 583         if(w.status){
 584 #ifdef DEBUG
 585                 fprintf(stdout, "\n");
 586                 hkl_geometry_fprintf(stdout, engine->geometry);
 587                 fprintf(stdout, "\n");
 588                 fprintf(stdout, "%s n : <%f, %f, %f> hkl : <%f, %f, %f> Q : <%f, %f, %f> angle : %f\n",
 589                         __func__,
 590                         w.n.data[0], w.n.data[1], w.n.data[2],
 591                         w.hkl.data[0], w.hkl.data[1], w.hkl.data[2],
 592                         w.Qn.data[0], w.Qn.data[1], w.Qn.data[2],
 593                         w.psi * HKL_RADTODEG);
 594 #endif
 595                 f->data[0] = w.hklW.dQ.data[0];
 596                 f->data[1] = w.hklW.dQ.data[1];
 597                 f->data[2] = w.hklW.dQ.data[2];
 598                 if(hkl_vector_norm2(&w.hkl2) > HKL_EPSILON)
 599                         f->data[3] = darray_item(engine->mode->parameters, 3)->_value - w.psi;
 600                 else
 601                         f->data[3] = 1;
 602         }else{
 603                 f->data[0] = 1;
 604                 f->data[1] = 1;
 605                 f->data[2] = 1;
 606                 f->data[3] = 1;
 607         }
 608
 609         return  GSL_SUCCESS;
 610 }
 611
 612 #define HKL_MODE_PSI_CONSTANT_VERTICAL_ERROR hkl_mode_psi_constant_vertical_error_quark ()
 613
 614 static GQuark hkl_mode_psi_constant_vertical_error_quark (void)
 615 {
 616         return g_quark_from_static_string ("hkl-mode-psi-constant-vertical-error-quark");
 617 }
 618
 619 typedef enum {
 620         HKL_MODE_PSI_CONSTANT_VERTICAL_ERROR_INITIALIZED_SET, /* can not init the engine */
 621 } HklModePsiConstantVerticalError;
 622
 623 int hkl_mode_initialized_set_psi_constant_vertical_real(HklMode *self,
 624                                                         HklEngine *engine,
 625                                                         HklGeometry *geometry,
 626                                                         HklDetector *detector,
 627                                                         HklSample *sample,
 628                                                         int initialized,
 629                                                         GError **error)
 630 {
 631         if(initialized){
 632                 HklEngineHkl *engine_hkl = container_of(engine, HklEngineHkl, engine);
 633                 const HklVector hkl = reciprocal_plan(engine_hkl);
 634                 const HklVector hkl2 = reciprocal_plan2(engine->mode);
 635
 636                 const struct HklPsiWrite w = hkl_psi_write(geometry, detector, sample, &hkl, &hkl2);
 637
 638                 /* kf - ki = Q */
 639                 if (hkl_vector_is_null(&w.hklW.Q)){
 640                         g_set_error(error,
 641                                     HKL_MODE_PSI_CONSTANT_VERTICAL_ERROR,
 642                                     HKL_MODE_PSI_CONSTANT_VERTICAL_ERROR_INITIALIZED_SET,
 643                                     "can not initialize the \"%s\" mode with a null hkl (kf == ki)"
 644                                     "\nplease select a non-null hkl", self->info->name);
 645                         return FALSE;
 646                 }else{
 647                         if (hkl_vector_is_null(&w.hkl2)){
 648                                 g_set_error(error,
 649                                             HKL_MODE_PSI_CONSTANT_VERTICAL_ERROR,
 650                                             HKL_MODE_PSI_CONSTANT_VERTICAL_ERROR_INITIALIZED_SET,
 651                                             "can not initialize the \"%s\" mode"
 652                                             "\nwhen Q and the <h2, k2, l2> ref vector are colinear."
 653                                             "\nplease change one or both of them", engine->mode->info->name);
 654                                 return FALSE;
 655                         }else{
 656                                 /* compute the angle beetween hkl and n and
 657                                  * store in in the fourth parameter */
 658                                 if (!hkl_parameter_value_set(darray_item(self->parameters, 3),
 659                                                              w.psi,
 660                                                              HKL_UNIT_DEFAULT, error))
 661                                         return FALSE;
 662                         }
 663                 }
 664         }
 665
 666         self->initialized = initialized;
 667
 668         return TRUE;
 669 }
 670
 671 /*******************/
 672 /* emergence fixed */
 673 /*******************/
 674
 675 typedef struct _HklModeAutoHklEmergenceFixed HklModeAutoHklEmergenceFixed;
 676
 677 struct _HklModeAutoHklEmergenceFixed
 678 {
 679         HklMode parent;
 680         HklParameter *n_x; /* not owned */
 681         HklParameter *n_y; /* not owned */
 682         HklParameter *n_z; /* not owned */
 683         HklParameter *emergence; /* not owned */
 684 };
 685
 686 #define HKL_MODE_HKL_EMERGENCE_FIXED_ERROR hkl_mode_hkl_emergence_fixed_error_quark ()
 687
 688 static GQuark hkl_mode_hkl_emergence_fixed_error_quark (void)
 689 {
 690         return g_quark_from_static_string ("hkl-mode-hkl-emergence-fixed-error-quark");
 691 }
 692
 693 typedef enum {
 694         HKL_MODE_HKL_EMERGENCE_FIXED_ERROR_INITIALIZED_SET, /* can not init the engine */
 695         HKL_MODE_HKL_EMERGENCE_FIXED_ERROR_SET, /* can not set the engine */
 696 } HklModeAutoHklEmergenceFixedError;
 697
 698
 699 static HklVector surface(const HklModeAutoHklEmergenceFixed *mode){
 700         HklVector n = {
 701                 .data = {
 702                         mode->n_x->_value,
 703                         mode->n_y->_value,
 704                         mode->n_z->_value,
 705                 }
 706         };
 707         return n;
 708 }
 709
 710 static double expected_emergence(const HklModeAutoHklEmergenceFixed *mode){
 711         return mode->emergence->_value;
 712 }
 713
 714 static int hkl_mode_hkl_emergence_fixed_initialized_set_real(HklMode *self,
 715                                                              HklEngine *engine,
 716                                                              HklGeometry *geometry,
 717                                                              HklDetector *detector,
 718                                                              HklSample *sample,
 719                                                              int initialized,
 720                                                              GError **error)
 721 {
 722         const HklModeAutoHklEmergenceFixed *mode = container_of(self, HklModeAutoHklEmergenceFixed, parent);
 723         const HklEngineHkl *engine_hkl = container_of(engine, HklEngineHkl, engine);
 724         const HklVector hkl = reciprocal_plan(engine_hkl);
 725         const HklVector n = surface(mode);
 726
 727         /* first check the parameters */
 728         if (hkl_vector_is_null(&n)){
 729                 g_set_error(error,
 730                             HKL_MODE_HKL_EMERGENCE_FIXED_ERROR,
 731                             HKL_MODE_HKL_EMERGENCE_FIXED_ERROR_INITIALIZED_SET,
 732                             "Can not compute emergence fixed when the surface vector is null.");
 733                 return FALSE;
 734         }
 735
 736         /* compute emergence and keep it */
 737         const struct HklEmergenceFixedWrite w = hkl_emergence_fixed_write(geometry, detector, sample, &hkl, &n);
 738         mode->emergence->_value = w.emergence;
 739
 740         self->initialized = initialized;
 741
 742         return TRUE;
 743 }
 744
 745
 746 int _emergence_fixed_func(const gsl_vector *x, void *params, gsl_vector *f)
 747 {
 748         CHECK_NAN(x->data, x->size);
 749
 750         HklEngine *engine = params;
 751         const HklEngineHkl *engine_hkl = container_of(engine, HklEngineHkl, engine);
 752         const HklModeAutoHklEmergenceFixed *mode = container_of(engine->mode, HklModeAutoHklEmergenceFixed, parent);
 753         const HklVector hkl = reciprocal_plan(engine_hkl);
 754         const HklVector n = surface(mode);
 755
 756         /* update the workspace from x; */
 757         set_geometry_axes(engine, x->data);
 758
 759         const struct HklEmergenceFixedWrite w = hkl_emergence_fixed_write(engine->geometry,
 760                                                                           engine->detector,
 761                                                                           engine->sample,
 762                                                                           &hkl, &n);
 763
 764         f->data[0] = w.hklW.dQ.data[0];
 765         f->data[1] = w.hklW.dQ.data[1];
 766         f->data[2] = w.hklW.dQ.data[2];
 767         f->data[3] = expected_emergence(mode) - w.emergence;
 768
 769         return  GSL_SUCCESS;
 770 }
 771
 772 int hkl_mode_hkl_emergence_fixed_set_real(HklMode *self,
 773                                           HklEngine *engine,
 774                                           HklGeometry *geometry,
 775                                           HklDetector *detector,
 776                                           HklSample *sample,
 777                                           GError **error)
 778 {
 779         const HklModeAutoHklEmergenceFixed *mode = container_of(self, HklModeAutoHklEmergenceFixed, parent);
 780         HklVector n = surface(mode);
 781
 782         /* first check the parameters */
 783         if (hkl_vector_is_null(&n)){
 784                 g_set_error(error,
 785                             HKL_MODE_HKL_EMERGENCE_FIXED_ERROR,
 786                             HKL_MODE_HKL_EMERGENCE_FIXED_ERROR_SET,
 787                             "Can not compute hkl with emergence fixed when the surface vector is null.");
 788                 return FALSE;
 789         }
 790
 791         return hkl_mode_set_hkl_real(self, engine, geometry, detector, sample, error);
 792 }
 793
 794 HklMode *hkl_mode_hkl_emergence_fixed_new(const HklModeAutoInfo *auto_info)
 795 {
 796         static const HklModeOperations operations = {
 797                 HKL_MODE_OPERATIONS_HKL_FULL_DEFAULTS,
 798                 .capabilities = HKL_ENGINE_CAPABILITIES_READABLE | HKL_ENGINE_CAPABILITIES_WRITABLE | HKL_ENGINE_CAPABILITIES_INITIALIZABLE,
 799                 .initialized_set = hkl_mode_hkl_emergence_fixed_initialized_set_real,
 800                 .set = hkl_mode_hkl_emergence_fixed_set_real,
 801         };
 802         HklModeAutoHklEmergenceFixed *self;
 803
 804         if (darray_size(auto_info->info.axes_w) != 4){
 805                 fprintf(stderr, "This generic HklModeAutoHklEmergenceFixed need exactly 4 axes");
 806                 exit(128);
 807         }
 808
 809         self = g_new(HklModeAutoHklEmergenceFixed, 1);
 810
 811         /* the base constructor; */
 812         hkl_mode_auto_init(&self->parent,
 813                            auto_info,
 814                            &operations, FALSE);
 815
 816         self->n_x = register_mode_parameter(&self->parent, 0);
 817         self->n_y = register_mode_parameter(&self->parent, 1);
 818         self->n_z = register_mode_parameter(&self->parent, 2);
 819         self->emergence = register_mode_parameter(&self->parent, 3);
 820
 821         return &self->parent;
 822 }
 823
 824 /*************/
 825 /* HklEngine */
 826 /*************/
 827
 828 static void hkl_engine_hkl_free_real(HklEngine *base)
 829 {
 830         HklEngineHkl *self = container_of(base, HklEngineHkl, engine);
 831         hkl_engine_release(&self->engine);
 832         free(self);
 833 }
 834
 835 HklEngine *hkl_engine_hkl_new(HklEngineList *engines)
 836 {
 837         HklEngineHkl *self = g_new(HklEngineHkl, 1);
 838         static const HklParameter h = {
 839                 HKL_PARAMETER_DEFAULTS, .name = "h",
 840                 .description = "h coordinate of the diffracting plan",
 841                 .range = { .min=-1, .max=1 },
 842         };
 843         static const HklParameter k = {
 844                 HKL_PARAMETER_DEFAULTS, .name = "k",
 845                 .description = "k coordinate of the diffracting plan",
 846                 .range = { .min=-1, .max=1 },
 847         };
 848         static const HklParameter l = {
 849                 HKL_PARAMETER_DEFAULTS, .name = "l",
 850                 .description = "l coordinate of the diffracting plan",
 851                 .range={ .min=-1, .max=1 },
 852         };
 853         static const HklParameter *pseudo_axes[] = {&h, &k, &l};
 854         static HklEngineInfo info = {
 855                 HKL_ENGINE_INFO("hkl",
 856                                 pseudo_axes,
 857                                 HKL_ENGINE_DEPENDENCIES_AXES | HKL_ENGINE_DEPENDENCIES_ENERGY | HKL_ENGINE_DEPENDENCIES_SAMPLE),
 858         };
 859         static HklEngineOperations operations = {
 860                 HKL_ENGINE_OPERATIONS_DEFAULTS,
 861                 .free=hkl_engine_hkl_free_real,
 862         };
 863
 864         hkl_engine_init(&self->engine, &info, &operations, engines);
 865
 866         self->h = register_pseudo_axis(&self->engine, engines, &h);
 867         self->k = register_pseudo_axis(&self->engine, engines, &k);
 868         self->l = register_pseudo_axis(&self->engine, engines, &l);
 869
 870         return &self->engine;
 871 }