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-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  *          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         ki = hkl_geometry_ki_get(geometry);
 105         kf = hkl_geometry_kf_get(geometry, detector);
 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 = g_new(HklEngineQ, 1);
 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         hkl_engine_init(&self->engine, &info, &operations, engines);
 165         self->q = register_pseudo_axis(&self->engine, engines, &q);
 166
 167         /* q [default] */
 168         mode = mode_q();
 169         hkl_engine_add_mode(&self->engine, mode);
 170         hkl_engine_mode_set(&self->engine, mode);
 171
 172         return &self->engine;
 173 }
 174
 175 /******/
 176 /* q2 */
 177 /******/
 178
 179 struct _HklEngineQ2
 180 {
 181         HklEngine engine;
 182         HklParameter *q;
 183         HklParameter *alpha;
 184 };
 185
 186 static void _q2(HklGeometry *geometry, HklDetector *detector,
 187                 double *q, double *alpha)
 188 {
 189         double wavelength, theta;
 190         HklVector kf, ki;
 191         static HklVector x = {
 192                 .data = {1, 0, 0},
 193         };
 194
 195         wavelength = hkl_source_get_wavelength(&geometry->source);
 196         ki = hkl_geometry_ki_get(geometry);
 197         kf = hkl_geometry_kf_get(geometry, detector);
 198         theta = hkl_vector_angle(&ki, &kf) / 2.;
 199
 200         *q = qmax(wavelength) * sin(theta);
 201
 202         /* project kf on the x plan to compute alpha */
 203         hkl_vector_project_on_plan(&kf, &x);
 204
 205         *alpha = atan2(kf.data[2], kf.data[1]);
 206 }
 207
 208 static int _q2_func(const gsl_vector *x, void *params, gsl_vector *f)
 209 {
 210         HklEngine *engine = params;
 211         const HklEngineQ2 *engine_q2 = container_of(engine, HklEngineQ2, engine);
 212         double q;
 213         double alpha;
 214
 215         CHECK_NAN(x->data, x->size);
 216
 217         /* update the workspace from x */
 218         set_geometry_axes(engine, x->data);
 219
 220         _q2(engine->geometry, engine->detector, &q, &alpha);
 221
 222         f->data[0] = engine_q2->q->_value - q;
 223         f->data[1] = engine_q2->alpha->_value - alpha;
 224
 225         return  GSL_SUCCESS;
 226 }
 227
 228 static const HklFunction q2_func = {
 229         .function = _q2_func,
 230         .size = 2,
 231 };
 232
 233 static int get_q2_real(HklMode *self,
 234                        HklEngine *engine,
 235                        HklGeometry *geometry,
 236                        HklDetector *detector,
 237                        HklSample *sample,
 238                        GError **error)
 239 {
 240         HklEngineQ2 *engine_q2 = container_of(engine, HklEngineQ2, engine);
 241
 242         _q2(geometry, detector, &engine_q2->q->_value, &engine_q2->alpha->_value);
 243
 244         return TRUE;
 245 }
 246
 247 static HklMode *mode_q2(void)
 248 {
 249         static const char* axes[] = {GAMMA, DELTA};
 250         static const HklFunction *functions[] = {&q2_func};
 251         static const HklModeAutoInfo info = {
 252                 HKL_MODE_AUTO_INFO("q2", axes, axes, functions),
 253         };
 254         static const HklModeOperations operations = {
 255                 HKL_MODE_OPERATIONS_AUTO_DEFAULTS,
 256                 .get = get_q2_real,
 257         };
 258
 259         return hkl_mode_auto_new(&info, &operations, TRUE);
 260 }
 261
 262 static const HklParameter alpha = {
 263         HKL_PARAMETER_DEFAULTS_ANGLE, .name = "alpha",
 264         .description = "angle of the projection of $\\vec{q}$ on the $yOz$ plan and $\\vec{y}$",
 265 };
 266
 267 static void hkl_engine_q2_free_real(HklEngine *base)
 268 {
 269         HklEngineQ2 *self = container_of(base, HklEngineQ2, engine);
 270         hkl_engine_release(&self->engine);
 271         free(self);
 272 }
 273
 274 HklEngine *hkl_engine_q2_new(HklEngineList *engines)
 275 {
 276         HklEngineQ2 *self = g_new(HklEngineQ2, 1);
 277         HklMode *mode;
 278         static const HklParameter *pseudo_axes[] = {&q, &alpha};
 279         static const HklEngineInfo info = {
 280                 HKL_ENGINE_INFO("q2",
 281                                 pseudo_axes,
 282                                 HKL_ENGINE_DEPENDENCIES_AXES | HKL_ENGINE_DEPENDENCIES_ENERGY),
 283         };
 284         static const HklEngineOperations operations = {
 285                 HKL_ENGINE_OPERATIONS_DEFAULTS,
 286                 .free=hkl_engine_q2_free_real,
 287         };
 288
 289         hkl_engine_init(&self->engine, &info, &operations, engines);
 290         self->q = register_pseudo_axis(&self->engine, engines, &q);
 291         self->alpha = register_pseudo_axis(&self->engine, engines, &alpha);
 292
 293         /* q2 [default] */
 294         mode = mode_q2();
 295         hkl_engine_add_mode(&self->engine, mode);
 296         hkl_engine_mode_set(&self->engine, mode);
 297
 298         return &self->engine;
 299 }
 300
 301 /************/
 302 /* QperQpar */
 303 /************/
 304
 305 typedef struct _HklModeIncidence HklModeQperQpar;
 306
 307 struct _HklEngineQperQpar
 308 {
 309         HklEngine engine;
 310         HklParameter *qper;
 311         HklParameter *qpar;
 312 };
 313
 314 static void _qper_qpar(HklEngine *engine,
 315                        const HklGeometry *geometry,
 316                        const HklDetector *detector,
 317                        const HklSample *sample,
 318                        double *qper, double *qpar)
 319 {
 320         HklModeQperQpar *mode = container_of(engine->mode, HklModeQperQpar, parent);
 321         HklHolder *sample_holder = hkl_geometry_sample_holder_get(geometry, sample);
 322         HklVector ki;
 323         HklVector q;
 324         HklVector n = {
 325                 .data = {
 326                         mode->n_x->_value,
 327                         mode->n_y->_value,
 328                         mode->n_z->_value,
 329                 },
 330         };
 331         HklVector npar;
 332         HklVector qper_v;
 333         HklVector qpar_v;
 334         double norm;
 335
 336         /* compute q = kf - ki */
 337         ki = hkl_geometry_ki_get(geometry);
 338         q = hkl_geometry_kf_get(geometry, detector);
 339         hkl_vector_minus_vector(&q, &ki);
 340
 341         /* compute the real orientation of the surface n */
 342         hkl_vector_rotated_quaternion(&n, &sample_holder->q);
 343         hkl_vector_normalize(&n);
 344
 345         /* compute the npar used to define the sign of qpar */
 346         npar = ki;
 347         hkl_vector_vectorial_product(&npar, &n);
 348
 349         /* qper */
 350         qper_v = n;
 351         norm = hkl_vector_scalar_product(&q, &n);
 352         hkl_vector_times_double(&qper_v, norm);
 353         *qper = hkl_vector_norm2(&qper_v);
 354         if (signbit(norm))
 355                 *qper *= -1;
 356
 357         /* qpar */
 358         qpar_v = q;
 359         norm = hkl_vector_scalar_product(&q, &npar);
 360         hkl_vector_minus_vector(&qpar_v, &qper_v);
 361         *qpar = hkl_vector_norm2(&qpar_v);
 362         if (signbit(norm))
 363                 *qpar *= -1;
 364 }
 365
 366 static int _qper_qpar_func(const gsl_vector *x, void *params, gsl_vector *f)
 367 {
 368         HklEngine *engine = params;
 369         const HklEngineQperQpar *engine_qper_qpar = container_of(engine, HklEngineQperQpar, engine);
 370         double qper;
 371         double qpar;
 372
 373         CHECK_NAN(x->data, x->size);
 374
 375         /* update the workspace from x */
 376         set_geometry_axes(engine, x->data);
 377
 378         _qper_qpar(engine, engine->geometry, engine->detector, engine->sample,
 379                    &qper, &qpar);
 380
 381         f->data[0] = engine_qper_qpar->qper->_value - qper;
 382         f->data[1] = engine_qper_qpar->qpar->_value - qpar;
 383
 384         return  GSL_SUCCESS;
 385 }
 386
 387 static const HklFunction qper_qpar_func = {
 388         .function = _qper_qpar_func,
 389         .size = 2,
 390 };
 391
 392 static int get_qper_qpar_real(HklMode *self,
 393                               HklEngine *engine,
 394                               HklGeometry *geometry,
 395                               HklDetector *detector,
 396                               HklSample *sample,
 397                               GError **error)
 398 {
 399         HklEngineQperQpar *engine_qper_qpar = container_of(engine, HklEngineQperQpar, engine);
 400
 401         _qper_qpar(engine, geometry, detector, sample,
 402                    &engine_qper_qpar->qper->_value,
 403                    &engine_qper_qpar->qpar->_value);
 404
 405         return TRUE;
 406 }
 407
 408 static HklMode *mode_qper_qpar(void)
 409 {
 410         static const char* axes[] = {GAMMA, DELTA};
 411         static const HklFunction *functions[] = {&qper_qpar_func};
 412         static const HklParameter parameters[] = {
 413                 SURFACE_PARAMETERS(0, 1, 0),
 414         };
 415         static const HklModeAutoInfo info = {
 416                 HKL_MODE_AUTO_INFO_WITH_PARAMS("qper_qpar", axes, axes, functions, parameters),
 417         };
 418         static const HklModeOperations operations = {
 419                 HKL_MODE_OPERATIONS_AUTO_DEFAULTS,
 420                 .get = get_qper_qpar_real,
 421         };
 422
 423         HklModeQperQpar *self = g_new(HklModeQperQpar, 1);
 424
 425         /* the base constructor; */
 426         hkl_mode_auto_init(&self->parent,
 427                            &info,
 428                            &operations, TRUE);
 429
 430         self->n_x = register_mode_parameter(&self->parent, 0);
 431         self->n_y = register_mode_parameter(&self->parent, 1);
 432         self->n_z = register_mode_parameter(&self->parent, 2);
 433
 434         return &self->parent;
 435 }
 436
 437 static void hkl_engine_qper_qpar_free_real(HklEngine *base)
 438 {
 439         HklEngineQperQpar *self = container_of(base, HklEngineQperQpar, engine);
 440         hkl_engine_release(&self->engine);
 441         free(self);
 442 }
 443
 444 HklEngine *hkl_engine_qper_qpar_new(HklEngineList *engines)
 445 {
 446         HklEngineQperQpar *self = g_new(HklEngineQperQpar, 1);
 447         HklMode *mode;
 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
 469         hkl_engine_init(&self->engine, &info, &operations, engines);
 470         self->qper = register_pseudo_axis(&self->engine, engines, &qper);
 471         self->qpar = register_pseudo_axis(&self->engine, engines, &qpar);
 472
 473         /* qper_qpar [default] */
 474         mode = mode_qper_qpar();
 475         hkl_engine_add_mode(&self->engine, mode);
 476         hkl_engine_mode_set(&self->engine, mode);
 477
 478         return &self->engine;
 479 }