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-2014 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 kf;
 440                         HklVector kf2;
 441                         HklVector q;
 442                         HklVector axis_v;
 443                         HklQuaternion qr;
 444                         HklAxis *axis;
 445                         HklVector cp = {0};
 446                         HklVector op = {0};
 447                         double angle;
 448                         HklGeometry *geom;
 449
 450                         geom = hkl_geometry_new_copy(item->geometry);
 451
 452                         /* get the Q vector kf - ki */
 453                         hkl_detector_compute_kf(detector, geom, &q);
 454                         hkl_source_compute_ki(&geom->source, &ki);
 455                         hkl_vector_minus_vector(&q, &ki);
 456
 457                         /* compute the current orientation of the last axis */
 458                         axis = container_of(darray_item(geom->axes,
 459                                                         darray_item(geom->holders, 0)->config->idx[last_axis]),
 460                                             HklAxis, parameter);
 461                         axis_v = axis->axis_v;
 462                         hkl_quaternion_init(&qr, 1, 0, 0, 0);
 463                         for(j=0; j<last_axis; ++j)
 464                                 hkl_quaternion_times_quaternion(
 465                                         &qr,
 466                                         &container_of(darray_item(geom->axes,
 467                                                                   darray_item(geom->holders, 0)->config->idx[j]),
 468                                                       HklAxis, parameter)->q);
 469                         hkl_vector_rotated_quaternion(&axis_v, &qr);
 470
 471                         /* - project the center of the ewalds sphere into the same plan (c') */
 472                         hkl_vector_minus_vector(&cp, &ki);
 473                         hkl_vector_project_on_plan_with_point(&cp, &axis_v, &q);
 474                         hkl_vector_project_on_plan_with_point(&op, &axis_v, &q);
 475
 476                         /* - rotate q around this (o', c') line of 180° to find the (q2) solution */
 477                         kf2 = q;
 478                         hkl_vector_rotated_around_line(&kf2, M_PI, &cp, &op);
 479                         angle = hkl_vector_oriented_angle_points(&q, &op, &kf2, &axis_v);
 480                         /* TODO parameter list for geometry */
 481                         if(!hkl_parameter_value_set(&axis->parameter,
 482                                                     hkl_parameter_value_get(&axis->parameter, HKL_UNIT_DEFAULT) + angle,
 483                                                     HKL_UNIT_DEFAULT, error))
 484                                 return FALSE;
 485                         hkl_geometry_update(geom);
 486 #ifdef DEBUG
 487                         fprintf(stdout, "\n- try to add a solution by rotating Q <%f, %f, %f> around the \"%s\" axis <%f, %f, %f> of %f radian",
 488                                 q.data[0], q.data[1], q.data[2],
 489                                 ((HklParameter *)axis)->name,
 490                                 axis_v.data[0], axis_v.data[1], axis_v.data[2],
 491                                 angle);
 492                         fprintf(stdout, "\n   op: <%f, %f, %f>", op.data[0], op.data[1], op.data[2]);
 493                         fprintf(stdout, "\n   q2: <%f, %f, %f>", kf2.data[0], kf2.data[1], kf2.data[2]);
 494 #endif
 495                         hkl_vector_add_vector(&kf2, &ki);
 496
 497                         /* at the end we just need to solve numerically the position of the detector */
 498                         if(fit_detector_position(self, geom, detector, &kf2))
 499                                 hkl_geometry_list_add(engine->engines->geometries,
 500                                                       geom);
 501
 502                         hkl_geometry_free(geom);
 503                 }
 504         }
 505         return TRUE;
 506 }
 507
 508 /***************************************/
 509 /* the double diffraction get set part */
 510 /***************************************/
 511
 512 /**
 513  * double_diffraction: (skip)
 514  * @x:
 515  * @params:
 516  * @f:
 517  *
 518  *
 519  *
 520  * Returns:
 521  **/
 522 int _double_diffraction(double const x[], void *params, double f[])
 523 {
 524         HklEngine *engine = params;
 525         HklEngineHkl *engine_hkl = container_of(engine, HklEngineHkl, engine);
 526         HklVector hkl = {
 527                 .data = {
 528                         engine_hkl->h->_value,
 529                         engine_hkl->k->_value,
 530                         engine_hkl->l->_value,
 531                 },
 532         };
 533         HklVector kf2;
 534         HklVector ki;
 535         HklVector dQ;
 536         HklHolder *sample_holder;
 537
 538         /* update the workspace from x; */
 539         set_geometry_axes(engine, x);
 540
 541         /* get the second hkl from the mode parameters */
 542         hkl_vector_init(&kf2,
 543                         darray_item(engine->mode->parameters, 0)->_value,
 544                         darray_item(engine->mode->parameters, 1)->_value,
 545                         darray_item(engine->mode->parameters, 2)->_value);
 546
 547         /* R * UB * hkl = Q */
 548         /* for now the 0 holder is the sample holder. */
 549         sample_holder = darray_item(engine->geometry->holders, 0);
 550         hkl_matrix_times_vector(&engine->sample->UB, &hkl);
 551         hkl_vector_rotated_quaternion(&hkl, &sample_holder->q);
 552
 553         /* kf - ki = Q */
 554         hkl_source_compute_ki(&engine->geometry->source, &ki);
 555         hkl_detector_compute_kf(engine->detector, engine->geometry, &dQ);
 556         hkl_vector_minus_vector(&dQ, &ki);
 557         hkl_vector_minus_vector(&dQ, &hkl);
 558
 559         /* R * UB * hlk2 = Q2 */
 560         hkl_matrix_times_vector(&engine->sample->UB, &kf2);
 561         hkl_vector_rotated_quaternion(&kf2, &sample_holder->q);
 562         hkl_vector_add_vector(&kf2, &ki);
 563
 564         f[0] = dQ.data[0];
 565         f[1] = dQ.data[1];
 566         f[2] = dQ.data[2];
 567         f[3] = hkl_vector_norm2(&kf2) - hkl_vector_norm2(&ki);
 568
 569         return GSL_SUCCESS;
 570 }
 571
 572 /**
 573  * double_diffraction_func: (skip)
 574  * @x:
 575  * @params:
 576  * @f:
 577  *
 578  *
 579  *
 580  * Returns:
 581  **/
 582 int _double_diffraction_func(gsl_vector const *x, void *params, gsl_vector *f)
 583 {
 584         CHECK_NAN(x->data, x->size);
 585
 586         _double_diffraction(x->data, params, f->data);
 587
 588         return  GSL_SUCCESS;
 589 }
 590
 591
 592 /******************************************/
 593 /* the psi_constant_vertical get set part */
 594 /******************************************/
 595
 596 /**
 597  * psi_constant_vertical_func: (skip)
 598  * @x:
 599  * @params:
 600  * @f:
 601  *
 602  *
 603  *
 604  * Returns:
 605  **/
 606 int _psi_constant_vertical_func(gsl_vector const *x, void *params, gsl_vector *f)
 607 {
 608         HklVector ki, kf, Q;
 609         HklEngine *engine = params;
 610
 611         CHECK_NAN(x->data, x->size);
 612
 613         RUBh_minus_Q(x->data, params, f->data);
 614
 615         /* update the workspace from x; */
 616         set_geometry_axes(engine, x->data);
 617
 618         /* kf - ki = Q */
 619         hkl_source_compute_ki(&engine->geometry->source, &ki);
 620         hkl_detector_compute_kf(engine->detector, engine->geometry, &kf);
 621         Q = kf;
 622         hkl_vector_minus_vector(&Q, &ki);
 623
 624         f->data[3] = darray_item(engine->mode->parameters, 3)->_value;
 625
 626         /* if |Q| > epsilon ok */
 627         if(hkl_vector_normalize(&Q)){
 628                 HklVector hkl;
 629                 HklVector n;
 630
 631                 /* compute n the intersection of the plan P(kf, ki) and PQ (normal Q) */
 632                 n = kf;
 633                 hkl_vector_vectorial_product(&n, &ki);
 634                 hkl_vector_vectorial_product(&n, &Q);
 635
 636                 /* compute the hkl ref position in the laboratory */
 637                 /* referentiel. The geometry was already updated. */
 638                 /* FIXME for now the 0 holder is the sample holder. */
 639                 hkl.data[0] = darray_item(engine->mode->parameters, 0)->_value;
 640                 hkl.data[1] = darray_item(engine->mode->parameters, 1)->_value;
 641                 hkl.data[2] = darray_item(engine->mode->parameters, 2)->_value;
 642                 hkl_matrix_times_vector(&engine->sample->UB, &hkl);
 643                 hkl_vector_rotated_quaternion(&hkl,
 644                                               &darray_item(engine->geometry->holders, 0)->q);
 645
 646                 /* project hkl on the plan of normal Q */
 647                 hkl_vector_project_on_plan(&hkl, &Q);
 648 #ifdef DEBUG
 649                 fprintf(stdout, "\n");
 650                 hkl_geometry_fprintf(stdout, engine->geometry);
 651                 fprintf(stdout, "\n");
 652                 fprintf(stdout, "%s n : <%f, %f, %f> hkl : <%f, %f, %f> Q : <%f, %f, %f> angle : %f\n",
 653                         __func__,
 654                         n.data[0], n.data[1], n.data[2],
 655                         hkl.data[0], hkl.data[1], hkl.data[2],
 656                         Q.data[0], Q.data[1], Q.data[2],
 657                         hkl_vector_oriented_angle(&n, &hkl, &Q) * HKL_RADTODEG);
 658 #endif
 659                 if(hkl_vector_norm2(&hkl) > HKL_EPSILON)
 660                         f->data[3] -=  hkl_vector_oriented_angle(&n, &hkl, &Q);
 661         }
 662
 663         return  GSL_SUCCESS;
 664 }
 665
 666 #define HKL_MODE_PSI_CONSTANT_VERTICAL_ERROR hkl_mode_psi_constant_vertical_error_quark ()
 667
 668 static GQuark hkl_mode_psi_constant_vertical_error_quark (void)
 669 {
 670         return g_quark_from_static_string ("hkl-mode-psi-constant-vertical-error-quark");
 671 }
 672
 673 typedef enum {
 674         HKL_MODE_PSI_CONSTANT_VERTICAL_ERROR_INITIALIZED_SET, /* can not init the engine */
 675 } HklModePsiConstantVerticalError;
 676
 677 int hkl_mode_initialized_set_psi_constant_vertical_real(HklMode *self,
 678                                                         HklEngine *engine,
 679                                                         HklGeometry *geometry,
 680                                                         HklDetector *detector,
 681                                                         HklSample *sample,
 682                                                         int initialized,
 683                                                         GError **error)
 684 {
 685         HklVector hkl;
 686         HklVector ki, kf, Q, n;
 687
 688         if(initialized){
 689                 /* kf - ki = Q */
 690                 hkl_source_compute_ki(&geometry->source, &ki);
 691                 hkl_detector_compute_kf(detector, geometry, &kf);
 692                 Q = kf;
 693                 hkl_vector_minus_vector(&Q, &ki);
 694
 695                 if (hkl_vector_is_null(&Q)){
 696                         g_set_error(error,
 697                                     HKL_MODE_PSI_CONSTANT_VERTICAL_ERROR,
 698                                     HKL_MODE_PSI_CONSTANT_VERTICAL_ERROR_INITIALIZED_SET,
 699                                     "can not initialize the \"%s\" mode with a null hkl (kf == ki)"
 700                                     "\nplease select a non-null hkl", self->info->name);
 701                         return FALSE;
 702                 }else{
 703                         /* needed for a problem of precision */
 704                         hkl_vector_normalize(&Q);
 705
 706                         /* compute the intersection of the plan P(kf, ki) and PQ (normal Q) */
 707                         n = kf;
 708                         hkl_vector_vectorial_product(&n, &ki);
 709                         hkl_vector_vectorial_product(&n, &Q);
 710
 711                         /* compute hkl in the laboratory referentiel */
 712                         /* the geometry was already updated in the detector compute kf */
 713                         /* for now the 0 holder is the sample holder */
 714                         hkl.data[0] = darray_item(self->parameters, 0)->_value;
 715                         hkl.data[1] = darray_item(self->parameters, 1)->_value;
 716                         hkl.data[2] = darray_item(self->parameters, 2)->_value;
 717                         hkl_matrix_times_vector(&sample->UB, &hkl);
 718                         hkl_vector_rotated_quaternion(&hkl,
 719                                                       &darray_item(geometry->holders, 0)->q);
 720
 721                         /* project hkl on the plan of normal Q */
 722                         hkl_vector_project_on_plan(&hkl, &Q);
 723
 724                         if (hkl_vector_is_null(&hkl)){
 725                                 g_set_error(error,
 726                                             HKL_MODE_PSI_CONSTANT_VERTICAL_ERROR,
 727                                             HKL_MODE_PSI_CONSTANT_VERTICAL_ERROR_INITIALIZED_SET,
 728                                             "can not initialize the \"%s\" mode"
 729                                             "\nwhen Q and the <h2, k2, l2> ref vector are colinear."
 730                                             "\nplease change one or both of them", engine->mode->info->name);
 731                                 return FALSE;
 732                         }else{
 733                                 /* compute the angle beetween hkl and n and
 734                                  * store in in the fourth parameter */
 735                                 if (!hkl_parameter_value_set(darray_item(self->parameters, 3),
 736                                                              hkl_vector_oriented_angle(&n, &hkl, &Q),
 737                                                              HKL_UNIT_DEFAULT, error))
 738                                         return FALSE;
 739                         }
 740                 }
 741         }
 742
 743         self->initialized = initialized;
 744
 745         return TRUE;
 746 }
 747
 748 /*************/
 749 /* HklEngine */
 750 /*************/
 751
 752 static void hkl_engine_hkl_free_real(HklEngine *base)
 753 {
 754         HklEngineHkl *self = container_of(base, HklEngineHkl, engine);
 755         hkl_engine_release(&self->engine);
 756         free(self);
 757 }
 758
 759 HklEngine *hkl_engine_hkl_new(void)
 760 {
 761         HklEngineHkl *self;
 762         static const HklPseudoAxis h = {
 763                 .parameter = { HKL_PARAMETER_DEFAULTS, .name = "h",
 764                                .description = "h coordinate of the diffracting plan",
 765                                .range = { .min=-1, .max=1 },
 766                 },
 767         };
 768         static const HklPseudoAxis k = {
 769                 .parameter = { HKL_PARAMETER_DEFAULTS, .name = "k",
 770                                .description = "k coordinate of the diffracting plan",
 771                                .range = { .min=-1, .max=1 },
 772                 },
 773         };
 774         static const HklPseudoAxis l = {
 775                 .parameter = { HKL_PARAMETER_DEFAULTS, .name = "l",
 776                                .description = "l coordinate of the diffracting plan",
 777                                .range={ .min=-1, .max=1 },
 778                 },
 779         };
 780         static const HklPseudoAxis *pseudo_axes[] = {&h, &k, &l};
 781         static HklEngineInfo info = {
 782                 HKL_ENGINE_INFO("hkl",
 783                                 pseudo_axes,
 784                                 HKL_ENGINE_DEPENDENCIES_AXES | HKL_ENGINE_DEPENDENCIES_ENERGY | HKL_ENGINE_DEPENDENCIES_SAMPLE),
 785         };
 786         static HklEngineOperations operations = {
 787                 HKL_ENGINE_OPERATIONS_DEFAULTS,
 788                 .free=hkl_engine_hkl_free_real,
 789         };
 790
 791         self = HKL_MALLOC(HklEngineHkl);
 792
 793         hkl_engine_init(&self->engine, &info, &operations);
 794
 795         self->h = register_pseudo_axis(&self->engine, &h.parameter);
 796         self->k = register_pseudo_axis(&self->engine, &k.parameter);
 797         self->l = register_pseudo_axis(&self->engine, &l.parameter);
 798
 799         return &self->engine;
 800 }