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-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  *          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-private.h"     // for _HklEngine, etc
  36 #include "hkl-source-private.h"         // for hkl_source_compute_ki, etc
  37 #include "hkl-vector-private.h"         // for HklVector, hkl_vector_angle, etc
  38 #include "hkl.h"                        // for HklEngine, HklParameter, etc
  39 #include "hkl/ccan/array_size/array_size.h"  // for ARRAY_SIZE
  40 #include "hkl/ccan/container_of/container_of.h"  // for container_of
  41 #include "hkl/ccan/darray/darray.h"     // for darray_item
  42
  43
  44 double qmax(double wavelength)
  45 {
  46         return 2 * HKL_TAU / wavelength;
  47 }
  48
  49 /*****/
  50 /* q */
  51 /*****/
  52
  53 struct _HklEngineQ
  54 {
  55         HklEngine engine;
  56         HklParameter *q;
  57 };
  58
  59 static int _q_func(const gsl_vector *x, void *params, gsl_vector *f)
  60 {
  61         double q;
  62         HklEngine *engine = params;
  63         const HklEngineQ *engine_q = container_of(engine, HklEngineQ, engine);
  64         double tth;
  65
  66         CHECK_NAN(x->data, x->size);
  67
  68         /* update the workspace from x */
  69         set_geometry_axes(engine, x->data);
  70
  71         tth = gsl_sf_angle_restrict_symm(x->data[0]);
  72         q = qmax(hkl_source_get_wavelength(&engine->geometry->source)) * sin(tth/2.);
  73
  74         f->data[0] = engine_q->q->_value - q;
  75
  76         return  GSL_SUCCESS;
  77 }
  78
  79 static const HklFunction q_func = {
  80         .function = _q_func,
  81         .size = 1,
  82 };
  83
  84 static int get_q_real(HklMode *self,
  85                       HklEngine *base,
  86                       HklGeometry *geometry,
  87                       HklDetector *detector,
  88                       HklSample *sample,
  89                       GError **error)
  90 {
  91         double wavelength;
  92         double theta;
  93         HklVector ki, kf;
  94         HklEngineQ *engine = container_of(base, HklEngineQ, engine);
  95
  96         wavelength = hkl_source_get_wavelength(&geometry->source);
  97         hkl_source_compute_ki(&geometry->source, &ki);
  98         hkl_detector_compute_kf(detector, geometry, &kf);
  99         theta = hkl_vector_angle(&ki, &kf) / 2.;
 100
 101         /* we decide of the sign of theta depending on the orientation
 102          * of kf in the direct-space */
 103         if(kf.data[1] < 0 || kf.data[2] < 0)
 104                 theta = -theta;
 105
 106         /* update q */
 107         engine->q->_value = qmax(wavelength) * sin(theta);
 108
 109         return TRUE;
 110 }
 111
 112 /* not declared in the constructor as it is used also in the q2 pseudo
 113  * axis engine */
 114 static const HklPseudoAxis q = {
 115         .parameter = {HKL_PARAMETER_DEFAULTS, .name="q", .range={.max=1}},
 116 };
 117
 118 static HklMode *mode_q(void)
 119 {
 120         static const char *axes[] = {"tth"};
 121         static const HklFunction *functions[] = {&q_func};
 122         static HklModeAutoInfo info = {
 123                 HKL_MODE_AUTO_INFO("q", axes, axes, functions),
 124         };
 125         static const HklModeOperations operations = {
 126                 HKL_MODE_OPERATIONS_AUTO_DEFAULTS,
 127                 .get = get_q_real,
 128         };
 129
 130         return hkl_mode_auto_new(&info, &operations, TRUE);
 131 }
 132
 133 static void hkl_engine_q_free_real(HklEngine *base)
 134 {
 135         HklEngineQ *self=container_of(base, HklEngineQ, engine);
 136         hkl_engine_release(&self->engine);
 137         free(self);
 138 }
 139
 140 HklEngine *hkl_engine_q_new(void)
 141 {
 142         HklEngineQ *self;
 143         HklMode *mode;
 144         static const HklPseudoAxis *pseudo_axes[] = {&q};
 145         static const HklEngineInfo info = {
 146                 .name = "q",
 147                 .pseudo_axes = DARRAY(pseudo_axes),
 148         };
 149         static const HklEngineOperations operations = {
 150                 HKL_ENGINE_OPERATIONS_DEFAULTS,
 151                 .free=hkl_engine_q_free_real,
 152         };
 153
 154         self = HKL_MALLOC(HklEngineQ);
 155
 156         hkl_engine_init(&self->engine, &info, &operations);
 157         self->q = register_pseudo_axis(&self->engine, &q.parameter);
 158
 159         /* q [default] */
 160         mode = mode_q();
 161         hkl_engine_add_mode(&self->engine, mode);
 162         hkl_engine_mode_set(&self->engine, mode);
 163
 164         return &self->engine;
 165 }
 166
 167 /******/
 168 /* q2 */
 169 /******/
 170
 171 struct _HklEngineQ2
 172 {
 173         HklEngine engine;
 174         HklParameter *q;
 175         HklParameter *alpha;
 176 };
 177
 178 static void _q2(HklGeometry *geometry, HklDetector *detector,
 179                 double *q, double *alpha)
 180 {
 181         double wavelength, theta;
 182         HklVector kf, ki;
 183         static HklVector x = {
 184                 .data = {1, 0, 0},
 185         };
 186
 187         wavelength = hkl_source_get_wavelength(&geometry->source);
 188         hkl_source_compute_ki(&geometry->source, &ki);
 189         hkl_detector_compute_kf(detector, geometry, &kf);
 190         theta = hkl_vector_angle(&ki, &kf) / 2.;
 191
 192         *q = qmax(wavelength) * sin(theta);
 193
 194         /* project kf on the x plan to compute alpha */
 195         hkl_vector_project_on_plan(&kf, &x);
 196
 197         *alpha = atan2(kf.data[2], kf.data[1]);
 198
 199 }
 200
 201 static int _q2_func(const gsl_vector *x, void *params, gsl_vector *f)
 202 {
 203         HklEngine *engine = params;
 204         const HklEngineQ2 *engine_q2 = container_of(engine, HklEngineQ2, engine);
 205         double q;
 206         double alpha;
 207
 208         CHECK_NAN(x->data, x->size);
 209
 210         /* update the workspace from x */
 211         set_geometry_axes(engine, x->data);
 212
 213         _q2(engine->geometry, engine->detector, &q, &alpha);
 214
 215         f->data[0] = engine_q2->q->_value - q;
 216         f->data[1] = engine_q2->alpha->_value - alpha;
 217
 218         return  GSL_SUCCESS;
 219 }
 220
 221 static const HklFunction q2_func = {
 222         .function = _q2_func,
 223         .size = 2,
 224 };
 225
 226 static int get_q2_real(HklMode *self,
 227                        HklEngine *engine,
 228                        HklGeometry *geometry,
 229                        HklDetector *detector,
 230                        HklSample *sample,
 231                        GError **error)
 232 {
 233         HklEngineQ2 *engine_q2 = container_of(engine, HklEngineQ2, engine);
 234
 235         _q2(geometry, detector, &engine_q2->q->_value, &engine_q2->alpha->_value);
 236
 237         return TRUE;
 238 }
 239
 240 static HklMode *mode_q2(void)
 241 {
 242         static const char* axes[] = {"gamma", "delta"};
 243         static const HklFunction *functions[] = {&q2_func};
 244         static const HklModeAutoInfo info = {
 245                 HKL_MODE_AUTO_INFO("q2", axes, axes, functions),
 246         };
 247         static const HklModeOperations operations = {
 248                 HKL_MODE_OPERATIONS_AUTO_DEFAULTS,
 249                 .get = get_q2_real,
 250         };
 251
 252         return hkl_mode_auto_new(&info, &operations, TRUE);
 253 }
 254
 255 static const HklPseudoAxis alpha = {
 256         .parameter = {HKL_PARAMETER_DEFAULTS_ANGLE, .name="alpha"},
 257 };
 258
 259 static void hkl_engine_q2_free_real(HklEngine *base)
 260 {
 261         HklEngineQ2 *self = container_of(base, HklEngineQ2, engine);
 262         hkl_engine_release(&self->engine);
 263         free(self);
 264 }
 265
 266 HklEngine *hkl_engine_q2_new(void)
 267 {
 268         HklEngineQ2 *self;
 269         HklMode *mode;
 270         static const HklPseudoAxis *pseudo_axes[] = {&q, &alpha};
 271         static const HklEngineInfo info = {
 272                 .name = "q2",
 273                 .pseudo_axes = DARRAY(pseudo_axes),
 274         };
 275         static const HklEngineOperations operations = {
 276                 HKL_ENGINE_OPERATIONS_DEFAULTS,
 277                 .free=hkl_engine_q2_free_real,
 278         };
 279
 280         self = HKL_MALLOC(HklEngineQ2);
 281
 282         hkl_engine_init(&self->engine, &info, &operations);
 283         self->q = register_pseudo_axis(&self->engine, &q.parameter);
 284         self->alpha = register_pseudo_axis(&self->engine, &alpha.parameter);
 285
 286         /* q2 [default] */
 287         mode = mode_q2();
 288         hkl_engine_add_mode(&self->engine, mode);
 289         hkl_engine_mode_set(&self->engine, mode);
 290
 291         return &self->engine;
 292 }
 293
 294 /************/
 295 /* QperQpar */
 296 /************/
 297
 298 struct _HklEngineQperQpar
 299 {
 300         HklEngine engine;
 301         HklParameter *qper;
 302         HklParameter *qpar;
 303 };
 304
 305 static void _qper_qpar(HklEngine *engine,
 306                        HklGeometry *geometry, HklDetector *detector,
 307                        double *qper, double *qpar)
 308 {
 309         HklVector ki;
 310         HklVector q;
 311         HklVector n = {
 312                 .data = {
 313                         darray_item(engine->mode->parameters, 0)->_value,
 314                         darray_item(engine->mode->parameters, 1)->_value,
 315                         darray_item(engine->mode->parameters, 2)->_value,
 316                 },
 317         };
 318         HklVector npar;
 319         HklVector qper_v;
 320         HklVector qpar_v;
 321         double norm;
 322
 323         /* compute q = kf - ki */
 324         hkl_source_compute_ki(&geometry->source, &ki);
 325         hkl_detector_compute_kf(detector, geometry, &q);
 326         hkl_vector_minus_vector(&q, &ki);
 327
 328         /* compute the real orientation of the surface n */
 329         hkl_vector_rotated_quaternion(&n, &darray_item(geometry->holders, 0)->q);
 330         hkl_vector_normalize(&n);
 331
 332         /* compute the npar used to define the sign of qpar */
 333         npar = ki;
 334         hkl_vector_vectorial_product(&npar, &n);
 335
 336         /* qper */
 337         qper_v = n;
 338         norm = hkl_vector_scalar_product(&q, &n);
 339         hkl_vector_times_double(&qper_v, norm);
 340         *qper = hkl_vector_norm2(&qper_v);
 341         if (signbit(norm))
 342                 *qper *= -1;
 343
 344         /* qpar */
 345         qpar_v = q;
 346         norm = hkl_vector_scalar_product(&q, &npar);
 347         hkl_vector_minus_vector(&qpar_v, &qper_v);
 348         *qpar = hkl_vector_norm2(&qpar_v);
 349         if (signbit(norm))
 350                 *qpar *= -1;
 351 }
 352
 353 static int _qper_qpar_func(const gsl_vector *x, void *params, gsl_vector *f)
 354 {
 355         HklEngine *engine = params;
 356         const HklEngineQperQpar *engine_qper_qpar = container_of(engine, HklEngineQperQpar, engine);
 357         double qper;
 358         double qpar;
 359
 360         CHECK_NAN(x->data, x->size);
 361
 362         /* update the workspace from x */
 363         set_geometry_axes(engine, x->data);
 364
 365         _qper_qpar(engine, engine->geometry, engine->detector,
 366                    &qper, &qpar);
 367
 368         f->data[0] = engine_qper_qpar->qper->_value - qper;
 369         f->data[1] = engine_qper_qpar->qpar->_value - qpar;
 370
 371         return  GSL_SUCCESS;
 372 }
 373
 374 static const HklFunction qper_qpar_func = {
 375         .function = _qper_qpar_func,
 376         .size = 2,
 377 };
 378
 379 static int get_qper_qpar_real(HklMode *self,
 380                               HklEngine *engine,
 381                               HklGeometry *geometry,
 382                               HklDetector *detector,
 383                               HklSample *sample,
 384                               GError **error)
 385 {
 386         HklEngineQperQpar *engine_qper_qpar = container_of(engine, HklEngineQperQpar, engine);
 387
 388         _qper_qpar(engine, geometry, detector,
 389                    &engine_qper_qpar->qper->_value,
 390                    &engine_qper_qpar->qpar->_value);
 391
 392         return TRUE;
 393 }
 394
 395 static HklMode *mode_qper_qpar(void)
 396 {
 397         static const char* axes[] = {"gamma", "delta"};
 398         static const HklFunction *functions[] = {&qper_qpar_func};
 399         static const HklParameter parameters[] = {
 400                 {HKL_PARAMETER_DEFAULTS, .name = "x", .range = {.min=-1, .max=1}, ._value = 0,},
 401                 {HKL_PARAMETER_DEFAULTS, .name = "y", .range = {.min=-1, .max=1}, ._value = 1,},
 402                 {HKL_PARAMETER_DEFAULTS, .name = "z", .range = {.min=-1, .max=1}, ._value = 0,},
 403         };
 404         static const HklModeAutoInfo info = {
 405                 HKL_MODE_AUTO_INFO_WITH_PARAMS("qper_qpar", axes, axes, functions, parameters),
 406         };
 407         static const HklModeOperations operations = {
 408                 HKL_MODE_OPERATIONS_AUTO_DEFAULTS,
 409                 .get = get_qper_qpar_real,
 410         };
 411
 412         return hkl_mode_auto_new(&info, &operations, TRUE);
 413 }
 414
 415 static void hkl_engine_qper_qpar_free_real(HklEngine *base)
 416 {
 417         HklEngineQperQpar *self = container_of(base, HklEngineQperQpar, engine);
 418         hkl_engine_release(&self->engine);
 419         free(self);
 420 }
 421
 422 HklEngine *hkl_engine_qper_qpar_new(void)
 423 {
 424         static const HklPseudoAxis qper = {
 425                 .parameter = {HKL_PARAMETER_DEFAULTS, .name="qper", .range={.min=-1, .max=1}},
 426         };
 427         static const HklPseudoAxis qpar = {
 428                 .parameter = {HKL_PARAMETER_DEFAULTS, .name="qpar", .range={.min=-1, .max=1}},
 429         };
 430         static const HklPseudoAxis *pseudo_axes[] = {&qper, &qpar};
 431         static const HklEngineInfo info = {
 432                 .name = "qper_qpar",
 433                 .pseudo_axes = DARRAY(pseudo_axes),
 434         };
 435         static const HklEngineOperations operations = {
 436                 HKL_ENGINE_OPERATIONS_DEFAULTS,
 437                 .free = hkl_engine_qper_qpar_free_real,
 438         };
 439         HklEngineQperQpar *self;
 440         HklMode *mode;
 441
 442         self = HKL_MALLOC(HklEngineQperQpar);
 443
 444         hkl_engine_init(&self->engine, &info, &operations);
 445         self->qper = register_pseudo_axis(&self->engine, &qper.parameter);
 446         self->qpar = register_pseudo_axis(&self->engine, &qpar.parameter);
 447
 448         /* qper_qpar [default] */
 449         mode = mode_qper_qpar();
 450         hkl_engine_add_mode(&self->engine, mode);
 451         hkl_engine_mode_set(&self->engine, mode);
 452
 453         return &self->engine;
 454 }