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",
 116                        .description = "the norm of $\\vec{q}$",
 117                        .range = { .max=1 },
 118         },
 119 };
 120
 121 static HklMode *mode_q(void)
 122 {
 123         static const char *axes[] = {"tth"};
 124         static const HklFunction *functions[] = {&q_func};
 125         static HklModeAutoInfo info = {
 126                 HKL_MODE_AUTO_INFO("q", axes, axes, functions),
 127         };
 128         static const HklModeOperations operations = {
 129                 HKL_MODE_OPERATIONS_AUTO_DEFAULTS,
 130                 .get = get_q_real,
 131         };
 132
 133         return hkl_mode_auto_new(&info, &operations, TRUE);
 134 }
 135
 136 static void hkl_engine_q_free_real(HklEngine *base)
 137 {
 138         HklEngineQ *self=container_of(base, HklEngineQ, engine);
 139         hkl_engine_release(&self->engine);
 140         free(self);
 141 }
 142
 143 HklEngine *hkl_engine_q_new(void)
 144 {
 145         HklEngineQ *self;
 146         HklMode *mode;
 147         static const HklPseudoAxis *pseudo_axes[] = {&q};
 148         static const HklEngineInfo info = {
 149                 .name = "q",
 150                 .pseudo_axes = DARRAY(pseudo_axes),
 151         };
 152         static const HklEngineOperations operations = {
 153                 HKL_ENGINE_OPERATIONS_DEFAULTS,
 154                 .free=hkl_engine_q_free_real,
 155         };
 156
 157         self = HKL_MALLOC(HklEngineQ);
 158
 159         hkl_engine_init(&self->engine, &info, &operations);
 160         self->q = register_pseudo_axis(&self->engine, &q.parameter);
 161
 162         /* q [default] */
 163         mode = mode_q();
 164         hkl_engine_add_mode(&self->engine, mode);
 165         hkl_engine_mode_set(&self->engine, mode);
 166
 167         return &self->engine;
 168 }
 169
 170 /******/
 171 /* q2 */
 172 /******/
 173
 174 struct _HklEngineQ2
 175 {
 176         HklEngine engine;
 177         HklParameter *q;
 178         HklParameter *alpha;
 179 };
 180
 181 static void _q2(HklGeometry *geometry, HklDetector *detector,
 182                 double *q, double *alpha)
 183 {
 184         double wavelength, theta;
 185         HklVector kf, ki;
 186         static HklVector x = {
 187                 .data = {1, 0, 0},
 188         };
 189
 190         wavelength = hkl_source_get_wavelength(&geometry->source);
 191         hkl_source_compute_ki(&geometry->source, &ki);
 192         hkl_detector_compute_kf(detector, geometry, &kf);
 193         theta = hkl_vector_angle(&ki, &kf) / 2.;
 194
 195         *q = qmax(wavelength) * sin(theta);
 196
 197         /* project kf on the x plan to compute alpha */
 198         hkl_vector_project_on_plan(&kf, &x);
 199
 200         *alpha = atan2(kf.data[2], kf.data[1]);
 201
 202 }
 203
 204 static int _q2_func(const gsl_vector *x, void *params, gsl_vector *f)
 205 {
 206         HklEngine *engine = params;
 207         const HklEngineQ2 *engine_q2 = container_of(engine, HklEngineQ2, engine);
 208         double q;
 209         double alpha;
 210
 211         CHECK_NAN(x->data, x->size);
 212
 213         /* update the workspace from x */
 214         set_geometry_axes(engine, x->data);
 215
 216         _q2(engine->geometry, engine->detector, &q, &alpha);
 217
 218         f->data[0] = engine_q2->q->_value - q;
 219         f->data[1] = engine_q2->alpha->_value - alpha;
 220
 221         return  GSL_SUCCESS;
 222 }
 223
 224 static const HklFunction q2_func = {
 225         .function = _q2_func,
 226         .size = 2,
 227 };
 228
 229 static int get_q2_real(HklMode *self,
 230                        HklEngine *engine,
 231                        HklGeometry *geometry,
 232                        HklDetector *detector,
 233                        HklSample *sample,
 234                        GError **error)
 235 {
 236         HklEngineQ2 *engine_q2 = container_of(engine, HklEngineQ2, engine);
 237
 238         _q2(geometry, detector, &engine_q2->q->_value, &engine_q2->alpha->_value);
 239
 240         return TRUE;
 241 }
 242
 243 static HklMode *mode_q2(void)
 244 {
 245         static const char* axes[] = {"gamma", "delta"};
 246         static const HklFunction *functions[] = {&q2_func};
 247         static const HklModeAutoInfo info = {
 248                 HKL_MODE_AUTO_INFO("q2", axes, axes, functions),
 249         };
 250         static const HklModeOperations operations = {
 251                 HKL_MODE_OPERATIONS_AUTO_DEFAULTS,
 252                 .get = get_q2_real,
 253         };
 254
 255         return hkl_mode_auto_new(&info, &operations, TRUE);
 256 }
 257
 258 static const HklPseudoAxis alpha = {
 259         .parameter = { HKL_PARAMETER_DEFAULTS_ANGLE, .name = "alpha",
 260                        .description = "angle of the projection of $\\vec{q}$ on the $yOz$ plan and $\\vec{y}$",
 261         },
 262 };
 263
 264 static void hkl_engine_q2_free_real(HklEngine *base)
 265 {
 266         HklEngineQ2 *self = container_of(base, HklEngineQ2, engine);
 267         hkl_engine_release(&self->engine);
 268         free(self);
 269 }
 270
 271 HklEngine *hkl_engine_q2_new(void)
 272 {
 273         HklEngineQ2 *self;
 274         HklMode *mode;
 275         static const HklPseudoAxis *pseudo_axes[] = {&q, &alpha};
 276         static const HklEngineInfo info = {
 277                 .name = "q2",
 278                 .pseudo_axes = DARRAY(pseudo_axes),
 279         };
 280         static const HklEngineOperations operations = {
 281                 HKL_ENGINE_OPERATIONS_DEFAULTS,
 282                 .free=hkl_engine_q2_free_real,
 283         };
 284
 285         self = HKL_MALLOC(HklEngineQ2);
 286
 287         hkl_engine_init(&self->engine, &info, &operations);
 288         self->q = register_pseudo_axis(&self->engine, &q.parameter);
 289         self->alpha = register_pseudo_axis(&self->engine, &alpha.parameter);
 290
 291         /* q2 [default] */
 292         mode = mode_q2();
 293         hkl_engine_add_mode(&self->engine, mode);
 294         hkl_engine_mode_set(&self->engine, mode);
 295
 296         return &self->engine;
 297 }
 298
 299 /************/
 300 /* QperQpar */
 301 /************/
 302
 303 struct _HklEngineQperQpar
 304 {
 305         HklEngine engine;
 306         HklParameter *qper;
 307         HklParameter *qpar;
 308 };
 309
 310 static void _qper_qpar(HklEngine *engine,
 311                        HklGeometry *geometry, HklDetector *detector,
 312                        double *qper, double *qpar)
 313 {
 314         HklVector ki;
 315         HklVector q;
 316         HklVector n = {
 317                 .data = {
 318                         darray_item(engine->mode->parameters, 0)->_value,
 319                         darray_item(engine->mode->parameters, 1)->_value,
 320                         darray_item(engine->mode->parameters, 2)->_value,
 321                 },
 322         };
 323         HklVector npar;
 324         HklVector qper_v;
 325         HklVector qpar_v;
 326         double norm;
 327
 328         /* compute q = kf - ki */
 329         hkl_source_compute_ki(&geometry->source, &ki);
 330         hkl_detector_compute_kf(detector, geometry, &q);
 331         hkl_vector_minus_vector(&q, &ki);
 332
 333         /* compute the real orientation of the surface n */
 334         hkl_vector_rotated_quaternion(&n, &darray_item(geometry->holders, 0)->q);
 335         hkl_vector_normalize(&n);
 336
 337         /* compute the npar used to define the sign of qpar */
 338         npar = ki;
 339         hkl_vector_vectorial_product(&npar, &n);
 340
 341         /* qper */
 342         qper_v = n;
 343         norm = hkl_vector_scalar_product(&q, &n);
 344         hkl_vector_times_double(&qper_v, norm);
 345         *qper = hkl_vector_norm2(&qper_v);
 346         if (signbit(norm))
 347                 *qper *= -1;
 348
 349         /* qpar */
 350         qpar_v = q;
 351         norm = hkl_vector_scalar_product(&q, &npar);
 352         hkl_vector_minus_vector(&qpar_v, &qper_v);
 353         *qpar = hkl_vector_norm2(&qpar_v);
 354         if (signbit(norm))
 355                 *qpar *= -1;
 356 }
 357
 358 static int _qper_qpar_func(const gsl_vector *x, void *params, gsl_vector *f)
 359 {
 360         HklEngine *engine = params;
 361         const HklEngineQperQpar *engine_qper_qpar = container_of(engine, HklEngineQperQpar, engine);
 362         double qper;
 363         double qpar;
 364
 365         CHECK_NAN(x->data, x->size);
 366
 367         /* update the workspace from x */
 368         set_geometry_axes(engine, x->data);
 369
 370         _qper_qpar(engine, engine->geometry, engine->detector,
 371                    &qper, &qpar);
 372
 373         f->data[0] = engine_qper_qpar->qper->_value - qper;
 374         f->data[1] = engine_qper_qpar->qpar->_value - qpar;
 375
 376         return  GSL_SUCCESS;
 377 }
 378
 379 static const HklFunction qper_qpar_func = {
 380         .function = _qper_qpar_func,
 381         .size = 2,
 382 };
 383
 384 static int get_qper_qpar_real(HklMode *self,
 385                               HklEngine *engine,
 386                               HklGeometry *geometry,
 387                               HklDetector *detector,
 388                               HklSample *sample,
 389                               GError **error)
 390 {
 391         HklEngineQperQpar *engine_qper_qpar = container_of(engine, HklEngineQperQpar, engine);
 392
 393         _qper_qpar(engine, geometry, detector,
 394                    &engine_qper_qpar->qper->_value,
 395                    &engine_qper_qpar->qpar->_value);
 396
 397         return TRUE;
 398 }
 399
 400 static HklMode *mode_qper_qpar(void)
 401 {
 402         static const char* axes[] = {"gamma", "delta"};
 403         static const HklFunction *functions[] = {&qper_qpar_func};
 404         static const HklParameter parameters[] = {
 405                 {
 406                         HKL_PARAMETER_DEFAULTS, .name = "x", ._value = 0,
 407                         .description = "the first coordinate of the surface vector",
 408                         .range = { .min=-1, .max=1 },
 409                 },
 410                 {
 411                         HKL_PARAMETER_DEFAULTS, .name = "y", ._value = 1,
 412                         .description = "the second coordinate of the surface vector",
 413                         .range = { .min=-1, .max=1 },
 414                 },
 415                 {
 416                         HKL_PARAMETER_DEFAULTS, .name = "z", ._value = 0,
 417                         .description = "the third coordinate of the surface vector",
 418                         .range = { .min=-1, .max=1 },
 419                 },
 420         };
 421         static const HklModeAutoInfo info = {
 422                 HKL_MODE_AUTO_INFO_WITH_PARAMS("qper_qpar", axes, axes, functions, parameters),
 423         };
 424         static const HklModeOperations operations = {
 425                 HKL_MODE_OPERATIONS_AUTO_DEFAULTS,
 426                 .get = get_qper_qpar_real,
 427         };
 428
 429         return hkl_mode_auto_new(&info, &operations, TRUE);
 430 }
 431
 432 static void hkl_engine_qper_qpar_free_real(HklEngine *base)
 433 {
 434         HklEngineQperQpar *self = container_of(base, HklEngineQperQpar, engine);
 435         hkl_engine_release(&self->engine);
 436         free(self);
 437 }
 438
 439 HklEngine *hkl_engine_qper_qpar_new(void)
 440 {
 441         static const HklPseudoAxis qper = {
 442                 .parameter = { HKL_PARAMETER_DEFAULTS, .name = "qper",
 443                                .description = "perpendicular component of $\\vec{q}$ along the normal of the sample surface",
 444                                .range = { .min=-1, .max=1 },
 445                 },
 446         };
 447         static const HklPseudoAxis qpar = {
 448                 .parameter = { HKL_PARAMETER_DEFAULTS, .name = "qpar",
 449                                .description = "parallel component of $\\vec{q}$",
 450                                .range = { .min=-1, .max=1 },
 451                 },
 452         };
 453         static const HklPseudoAxis *pseudo_axes[] = {&qper, &qpar};
 454         static const HklEngineInfo info = {
 455                 .name = "qper_qpar",
 456                 .pseudo_axes = DARRAY(pseudo_axes),
 457         };
 458         static const HklEngineOperations operations = {
 459                 HKL_ENGINE_OPERATIONS_DEFAULTS,
 460                 .free = hkl_engine_qper_qpar_free_real,
 461         };
 462         HklEngineQperQpar *self;
 463         HklMode *mode;
 464
 465         self = HKL_MALLOC(HklEngineQperQpar);
 466
 467         hkl_engine_init(&self->engine, &info, &operations);
 468         self->qper = register_pseudo_axis(&self->engine, &qper.parameter);
 469         self->qpar = register_pseudo_axis(&self->engine, &qpar.parameter);
 470
 471         /* qper_qpar [default] */
 472         mode = mode_qper_qpar();
 473         hkl_engine_add_mode(&self->engine, mode);
 474         hkl_engine_mode_set(&self->engine, mode);
 475
 476         return &self->engine;
 477 }