hkl/hkl-pseudoaxis-common-q.c

   1 /* This file is part of the hkl library.
   2  *
   3  * The hkl library is free software: you can redistribute it and/or modify
   4  * it under the terms of the GNU General Public License as published by
   5  * the Free Software Foundation, either version 3 of the License, or
   6  * (at your option) any later version.
   7  *
   8  * The hkl library is distributed in the hope that it will be useful,
   9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11  * GNU General Public License for more details.
  12  *
  13  * You should have received a copy of the GNU General Public License
  14  * along with the hkl library.  If not, see <http://www.gnu.org/licenses/>.
  15  *
  16  * Copyright (C) 2003-2016 Synchrotron SOLEIL
  17  *                         L'Orme des Merisiers Saint-Aubin
  18  *                         BP 48 91192 GIF-sur-YVETTE CEDEX
  19  *
  20  * Authors: Picca Frédéric-Emmanuel <picca@synchrotron-soleil.fr>
  21  *          Jens Krüger <Jens.Krueger@frm2.tum.de>
  22  */
  23 #include <gsl/gsl_errno.h>              // for ::GSL_SUCCESS
  24 #include <gsl/gsl_sf_trig.h>            // for gsl_sf_angle_restrict_symm
  25 #include <gsl/gsl_sys.h>                // for gsl_isnan
  26 #include <gsl/gsl_vector_double.h>      // for gsl_vector
  27 #include <math.h>                       // for sin, atan2, signbit
  28 #include <stdlib.h>                     // for free
  29 #include "hkl-detector-private.h"       // for hkl_detector_compute_kf
  30 #include "hkl-geometry-private.h"       // for _HklGeometry, HklHolder
  31 #include "hkl-macros-private.h"         // for HKL_MALLOC
  32 #include "hkl-parameter-private.h"      // for _HklParameter, etc
  33 #include "hkl-pseudoaxis-auto-private.h"  // for HklFunction, etc
  34 #include "hkl-pseudoaxis-common-q-private.h"  // for HklEngineQ2, etc
  35 #include "hkl-pseudoaxis-common-readonly-private.h"
  36 #include "hkl-pseudoaxis-private.h"     // for _HklEngine, etc
  37 #include "hkl-source-private.h"         // for hkl_source_compute_ki, etc
  38 #include "hkl-vector-private.h"         // for HklVector, hkl_vector_angle, etc
  39 #include "hkl.h"                        // for HklEngine, HklParameter, etc
  40 #include "hkl/ccan/array_size/array_size.h"  // for ARRAY_SIZE
  41 #include "hkl/ccan/container_of/container_of.h"  // for container_of
  42 #include "hkl/ccan/darray/darray.h"     // for darray_item
  43
  44 #define GAMMA "gamma"
  45 #define DELTA "delta"
  46
  47 typedef struct _HklEngineQ HklEngineQ;
  48 typedef struct _HklEngineQ2 HklEngineQ2;
  49 typedef struct _HklEngineQperQpar HklEngineQperQpar;
  50
  51 double qmax(double wavelength)
  52 {
  53         return 2 * HKL_TAU / wavelength;
  54 }
  55
  56 /*****/
  57 /* q */
  58 /*****/
  59
  60 struct _HklEngineQ
  61 {
  62         HklEngine engine;
  63         HklParameter *q;
  64 };
  65
  66 static int _q_func(const gsl_vector *x, void *params, gsl_vector *f)
  67 {
  68         double q;
  69         HklEngine *engine = params;
  70         const HklEngineQ *engine_q = container_of(engine, HklEngineQ, engine);
  71         double tth;
  72
  73         CHECK_NAN(x->data, x->size);
  74
  75         /* update the workspace from x */
  76         set_geometry_axes(engine, x->data);
  77
  78         tth = gsl_sf_angle_restrict_symm(x->data[0]);
  79         q = qmax(hkl_source_get_wavelength(&engine->geometry->source)) * sin(tth/2.);
  80
  81         f->data[0] = engine_q->q->_value - q;
  82
  83         return  GSL_SUCCESS;
  84 }
  85
  86 static const HklFunction q_func = {
  87         .function = _q_func,
  88         .size = 1,
  89 };
  90
  91 static int get_q_real(HklMode *self,
  92                       HklEngine *base,
  93                       HklGeometry *geometry,
  94                       HklDetector *detector,
  95                       HklSample *sample,
  96                       GError **error)
  97 {
  98         double wavelength;
  99         double theta;
 100         HklVector ki, kf;
 101         HklEngineQ *engine = container_of(base, HklEngineQ, engine);
 102
 103         wavelength = hkl_source_get_wavelength(&geometry->source);
 104         hkl_source_compute_ki(&geometry->source, &ki);
 105         hkl_detector_compute_kf(detector, geometry, &kf);
 106         theta = hkl_vector_angle(&ki, &kf) / 2.;
 107
 108         /* we decide of the sign of theta depending on the orientation
 109          * of kf in the direct-space */
 110         if(kf.data[1] < 0 || kf.data[2] < 0)
 111                 theta = -theta;
 112
 113         /* update q */
 114         engine->q->_value = qmax(wavelength) * sin(theta);
 115
 116         return TRUE;
 117 }
 118
 119 /* not declared in the constructor as it is used also in the q2 pseudo
 120  * axis engine */
 121 static const HklParameter q = {
 122         HKL_PARAMETER_DEFAULTS, .name="q",
 123         .description = "the norm of $\\vec{q}$",
 124         .range = { .max=1 },
 125 };
 126
 127 static HklMode *mode_q(void)
 128 {
 129         static const char *axes[] = {"tth"};
 130         static const HklFunction *functions[] = {&q_func};
 131         static HklModeAutoInfo info = {
 132                 HKL_MODE_AUTO_INFO("q", axes, axes, functions),
 133         };
 134         static const HklModeOperations operations = {
 135                 HKL_MODE_OPERATIONS_AUTO_DEFAULTS,
 136                 .get = get_q_real,
 137         };
 138
 139         return hkl_mode_auto_new(&info, &operations, TRUE);
 140 }
 141
 142 static void hkl_engine_q_free_real(HklEngine *base)
 143 {
 144         HklEngineQ *self=container_of(base, HklEngineQ, engine);
 145         hkl_engine_release(&self->engine);
 146         free(self);
 147 }
 148
 149 HklEngine *hkl_engine_q_new(HklEngineList *engines)
 150 {
 151         HklEngineQ *self;
 152         HklMode *mode;
 153         static const HklParameter *pseudo_axes[] = {&q};
 154         static const HklEngineInfo info = {
 155                 HKL_ENGINE_INFO("q",
 156                                 pseudo_axes,
 157                                 HKL_ENGINE_DEPENDENCIES_AXES | HKL_ENGINE_DEPENDENCIES_ENERGY),
 158         };
 159         static const HklEngineOperations operations = {
 160                 HKL_ENGINE_OPERATIONS_DEFAULTS,
 161                 .free=hkl_engine_q_free_real,
 162         };
 163
 164         self = HKL_MALLOC(HklEngineQ);
 165
 166         hkl_engine_init(&self->engine, &info, &operations, engines);
 167         self->q = register_pseudo_axis(&self->engine, engines, &q);
 168
 169         /* q [default] */
 170         mode = mode_q();
 171         hkl_engine_add_mode(&self->engine, mode);
 172         hkl_engine_mode_set(&self->engine, mode);
 173
 174         return &self->engine;
 175 }
 176
 177 /******/
 178 /* q2 */
 179 /******/
 180
 181 struct _HklEngineQ2
 182 {
 183         HklEngine engine;
 184         HklParameter *q;
 185         HklParameter *alpha;
 186 };
 187
 188 static void _q2(HklGeometry *geometry, HklDetector *detector,
 189                 double *q, double *alpha)
 190 {
 191         double wavelength, theta;
 192         HklVector kf, ki;
 193         static HklVector x = {
 194                 .data = {1, 0, 0},
 195         };
 196
 197         wavelength = hkl_source_get_wavelength(&geometry->source);
 198         hkl_source_compute_ki(&geometry->source, &ki);
 199         hkl_detector_compute_kf(detector, geometry, &kf);
 200         theta = hkl_vector_angle(&ki, &kf) / 2.;
 201
 202         *q = qmax(wavelength) * sin(theta);
 203
 204         /* project kf on the x plan to compute alpha */
 205         hkl_vector_project_on_plan(&kf, &x);
 206
 207         *alpha = atan2(kf.data[2], kf.data[1]);
 208
 209 }
 210
 211 static int _q2_func(const gsl_vector *x, void *params, gsl_vector *f)
 212 {
 213         HklEngine *engine = params;
 214         const HklEngineQ2 *engine_q2 = container_of(engine, HklEngineQ2, engine);
 215         double q;
 216         double alpha;
 217
 218         CHECK_NAN(x->data, x->size);
 219
 220         /* update the workspace from x */
 221         set_geometry_axes(engine, x->data);
 222
 223         _q2(engine->geometry, engine->detector, &q, &alpha);
 224
 225         f->data[0] = engine_q2->q->_value - q;
 226         f->data[1] = engine_q2->alpha->_value - alpha;
 227
 228         return  GSL_SUCCESS;
 229 }
 230
 231 static const HklFunction q2_func = {
 232         .function = _q2_func,
 233         .size = 2,
 234 };
 235
 236 static int get_q2_real(HklMode *self,
 237                        HklEngine *engine,
 238                        HklGeometry *geometry,
 239                        HklDetector *detector,
 240                        HklSample *sample,
 241                        GError **error)
 242 {
 243         HklEngineQ2 *engine_q2 = container_of(engine, HklEngineQ2, engine);
 244
 245         _q2(geometry, detector, &engine_q2->q->_value, &engine_q2->alpha->_value);
 246
 247         return TRUE;
 248 }
 249
 250 static HklMode *mode_q2(void)
 251 {
 252         static const char* axes[] = {GAMMA, DELTA};
 253         static const HklFunction *functions[] = {&q2_func};
 254         static const HklModeAutoInfo info = {
 255                 HKL_MODE_AUTO_INFO("q2", axes, axes, functions),
 256         };
 257         static const HklModeOperations operations = {
 258                 HKL_MODE_OPERATIONS_AUTO_DEFAULTS,
 259                 .get = get_q2_real,
 260         };
 261
 262         return hkl_mode_auto_new(&info, &operations, TRUE);
 263 }
 264
 265 static const HklParameter alpha = {
 266         HKL_PARAMETER_DEFAULTS_ANGLE, .name = "alpha",
 267         .description = "angle of the projection of $\\vec{q}$ on the $yOz$ plan and $\\vec{y}$",
 268 };
 269
 270 static void hkl_engine_q2_free_real(HklEngine *base)
 271 {
 272         HklEngineQ2 *self = container_of(base, HklEngineQ2, engine);
 273         hkl_engine_release(&self->engine);
 274         free(self);
 275 }
 276
 277 HklEngine *hkl_engine_q2_new(HklEngineList *engines)
 278 {
 279         HklEngineQ2 *self;
 280         HklMode *mode;
 281         static const HklParameter *pseudo_axes[] = {&q, &alpha};
 282         static const HklEngineInfo info = {
 283                 HKL_ENGINE_INFO("q2",
 284                                 pseudo_axes,
 285                                 HKL_ENGINE_DEPENDENCIES_AXES | HKL_ENGINE_DEPENDENCIES_ENERGY),
 286         };
 287         static const HklEngineOperations operations = {
 288                 HKL_ENGINE_OPERATIONS_DEFAULTS,
 289                 .free=hkl_engine_q2_free_real,
 290         };
 291
 292         self = HKL_MALLOC(HklEngineQ2);
 293
 294         hkl_engine_init(&self->engine, &info, &operations, engines);
 295         self->q = register_pseudo_axis(&self->engine, engines, &q);
 296         self->alpha = register_pseudo_axis(&self->engine, engines, &alpha);
 297
 298         /* q2 [default] */
 299         mode = mode_q2();
 300         hkl_engine_add_mode(&self->engine, mode);
 301         hkl_engine_mode_set(&self->engine, mode);
 302
 303         return &self->engine;
 304 }
 305
 306 /************/
 307 /* QperQpar */
 308 /************/
 309
 310 typedef struct _HklModeIncidence HklModeQperQpar;
 311
 312 struct _HklEngineQperQpar
 313 {
 314         HklEngine engine;
 315         HklParameter *qper;
 316         HklParameter *qpar;
 317 };
 318
 319 static void _qper_qpar(HklEngine *engine,
 320                        HklGeometry *geometry, HklDetector *detector,
 321                        double *qper, double *qpar)
 322 {
 323         HklModeQperQpar *mode = container_of(engine->mode, HklModeQperQpar, parent);
 324         HklVector ki;
 325         HklVector q;
 326         HklVector n = {
 327                 .data = {
 328                         mode->n_x->_value,
 329                         mode->n_y->_value,
 330                         mode->n_z->_value,
 331                 },
 332         };
 333         HklVector npar;
 334         HklVector qper_v;
 335         HklVector qpar_v;
 336         double norm;
 337
 338         /* compute q = kf - ki */
 339         hkl_source_compute_ki(&geometry->source, &ki);
 340         hkl_detector_compute_kf(detector, geometry, &q);
 341         hkl_vector_minus_vector(&q, &ki);
 342
 343         /* compute the real orientation of the surface n */
 344         hkl_vector_rotated_quaternion(&n, &darray_item(geometry->holders, 0)->q);
 345         hkl_vector_normalize(&n);
 346
 347         /* compute the npar used to define the sign of qpar */
 348         npar = ki;
 349         hkl_vector_vectorial_product(&npar, &n);
 350
 351         /* qper */
 352         qper_v = n;
 353         norm = hkl_vector_scalar_product(&q, &n);
 354         hkl_vector_times_double(&qper_v, norm);
 355         *qper = hkl_vector_norm2(&qper_v);
 356         if (signbit(norm))
 357                 *qper *= -1;
 358
 359         /* qpar */
 360         qpar_v = q;
 361         norm = hkl_vector_scalar_product(&q, &npar);
 362         hkl_vector_minus_vector(&qpar_v, &qper_v);
 363         *qpar = hkl_vector_norm2(&qpar_v);
 364         if (signbit(norm))
 365                 *qpar *= -1;
 366 }
 367
 368 static int _qper_qpar_func(const gsl_vector *x, void *params, gsl_vector *f)
 369 {
 370         HklEngine *engine = params;
 371         const HklEngineQperQpar *engine_qper_qpar = container_of(engine, HklEngineQperQpar, engine);
 372         double qper;
 373         double qpar;
 374
 375         CHECK_NAN(x->data, x->size);
 376
 377         /* update the workspace from x */
 378         set_geometry_axes(engine, x->data);
 379
 380         _qper_qpar(engine, engine->geometry, engine->detector,
 381                    &qper, &qpar);
 382
 383         f->data[0] = engine_qper_qpar->qper->_value - qper;
 384         f->data[1] = engine_qper_qpar->qpar->_value - qpar;
 385
 386         return  GSL_SUCCESS;
 387 }
 388
 389 static const HklFunction qper_qpar_func = {
 390         .function = _qper_qpar_func,
 391         .size = 2,
 392 };
 393
 394 static int get_qper_qpar_real(HklMode *self,
 395                               HklEngine *engine,
 396                               HklGeometry *geometry,
 397                               HklDetector *detector,
 398                               HklSample *sample,
 399                               GError **error)
 400 {
 401         HklEngineQperQpar *engine_qper_qpar = container_of(engine, HklEngineQperQpar, engine);
 402
 403         _qper_qpar(engine, geometry, detector,
 404                    &engine_qper_qpar->qper->_value,
 405                    &engine_qper_qpar->qpar->_value);
 406
 407         return TRUE;
 408 }
 409
 410 static HklMode *mode_qper_qpar(void)
 411 {
 412         static const char* axes[] = {GAMMA, DELTA};
 413         static const HklFunction *functions[] = {&qper_qpar_func};
 414         static const HklParameter parameters[] = {
 415                 SURFACE_PARAMETERS(0, 1, 0),
 416         };
 417         static const HklModeAutoInfo info = {
 418                 HKL_MODE_AUTO_INFO_WITH_PARAMS("qper_qpar", axes, axes, functions, parameters),
 419         };
 420         static const HklModeOperations operations = {
 421                 HKL_MODE_OPERATIONS_AUTO_DEFAULTS,
 422                 .get = get_qper_qpar_real,
 423         };
 424
 425         HklModeQperQpar *self = HKL_MALLOC(HklModeQperQpar);
 426
 427         /* the base constructor; */
 428         hkl_mode_auto_init(&self->parent,
 429                            &info,
 430                            &operations, TRUE);
 431
 432         self->n_x = register_mode_parameter(&self->parent, 0);
 433         self->n_y = register_mode_parameter(&self->parent, 1);
 434         self->n_z = register_mode_parameter(&self->parent, 2);
 435
 436         return &self->parent;
 437 }
 438
 439 static void hkl_engine_qper_qpar_free_real(HklEngine *base)
 440 {
 441         HklEngineQperQpar *self = container_of(base, HklEngineQperQpar, engine);
 442         hkl_engine_release(&self->engine);
 443         free(self);
 444 }
 445
 446 HklEngine *hkl_engine_qper_qpar_new(HklEngineList *engines)
 447 {
 448         static const HklParameter qper = {
 449                 HKL_PARAMETER_DEFAULTS, .name = "qper",
 450                 .description = "perpendicular component of $\\vec{q}$ along the normal of the sample surface",
 451                 .range = { .min=-1, .max=1 },
 452         };
 453         static const HklParameter qpar = {
 454                 HKL_PARAMETER_DEFAULTS, .name = "qpar",
 455                 .description = "parallel component of $\\vec{q}$",
 456                 .range = { .min=-1, .max=1 },
 457         };
 458         static const HklParameter *pseudo_axes[] = {&qper, &qpar};
 459         static const HklEngineInfo info = {
 460                 HKL_ENGINE_INFO("qper_qpar",
 461                                 pseudo_axes,
 462                                 HKL_ENGINE_DEPENDENCIES_AXES | HKL_ENGINE_DEPENDENCIES_ENERGY),
 463         };
 464         static const HklEngineOperations operations = {
 465                 HKL_ENGINE_OPERATIONS_DEFAULTS,
 466                 .free = hkl_engine_qper_qpar_free_real,
 467         };
 468         HklEngineQperQpar *self;
 469         HklMode *mode;
 470
 471         self = HKL_MALLOC(HklEngineQperQpar);
 472
 473         hkl_engine_init(&self->engine, &info, &operations, engines);
 474         self->qper = register_pseudo_axis(&self->engine, engines, &qper);
 475         self->qpar = register_pseudo_axis(&self->engine, engines, &qpar);
 476
 477         /* qper_qpar [default] */
 478         mode = mode_qper_qpar();
 479         hkl_engine_add_mode(&self->engine, mode);
 480         hkl_engine_mode_set(&self->engine, mode);
 481
 482         return &self->engine;
 483 }