hkl/hkl-pseudoaxis-common-psi.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-2017 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  */
  22 #include <gsl/gsl_errno.h>              // for ::GSL_SUCCESS
  23 #include <gsl/gsl_vector_double.h>      // for gsl_vector
  24 #include <gsl/gsl_sys.h>                // for gsl_isnan
  25 #include <stdio.h>                      // for fprintf, stderr
  26 #include <stdlib.h>                     // for NULL, exit, free
  27 #include <sys/types.h>                  // for uint
  28 #include "hkl-detector-private.h"       // for hkl_detector_compute_kf
  29 #include "hkl-geometry-private.h"       // for HklHolder, _HklGeometry, etc
  30 #include "hkl-macros-private.h"         // for HKL_MALLOC, hkl_assert, etc
  31 #include "hkl-matrix-private.h"         // for hkl_matrix_solve, etc
  32 #include "hkl-parameter-private.h"
  33 #include "hkl-pseudoaxis-auto-private.h"  // for HklModeAutoInfo, etc
  34 #include "hkl-pseudoaxis-common-psi-private.h"  // for HklEnginePsi, etc
  35 #include "hkl-pseudoaxis-private.h"     // for _HklEngine, HklEngineInfo, etc
  36 #include "hkl-quaternion-private.h"     // for hkl_quaternion_to_matrix
  37 #include "hkl-sample-private.h"         // for _HklSample
  38 #include "hkl-source-private.h"         // for hkl_source_compute_ki
  39 #include "hkl-vector-private.h"         // for HklVector, etc
  40 #include "hkl.h"                        // for HklEngine, HklGeometry, etc
  41 #include "hkl/ccan/array_size/array_size.h"  // for ARRAY_SIZE
  42 #include "hkl/ccan/container_of/container_of.h"  // for container_of
  43 #include "hkl/ccan/darray/darray.h"     // for darray_item
  44
  45 #define HKL_MODE_PSI_ERROR hkl_mode_psi_error_quark ()
  46
  47 static GQuark hkl_mode_psi_error_quark (void)
  48 {
  49         return g_quark_from_static_string ("hkl-mode-psi-error-quark");
  50 }
  51
  52 typedef enum {
  53         HKL_MODE_PSI_ERROR_INIT, /* can not init the engine */
  54         HKL_MODE_PSI_ERROR_GET, /* can not get the engine */
  55 } HklModePsiError;
  56
  57 /***********************/
  58 /* numerical functions */
  59 /***********************/
  60
  61 int _psi_func(const gsl_vector *x, void *params, gsl_vector *f)
  62 {
  63
  64         HklVector dhkl0, hkl1;
  65         HklVector ki, kf, Q, n;
  66         HklMatrix RUB;
  67         HklEngine *engine = params;
  68         HklEnginePsi *psi_engine = container_of(engine, HklEnginePsi, engine);
  69         HklModePsi *modepsi = container_of(engine->mode, HklModePsi, parent);
  70         HklHolder *sample_holder;
  71
  72         CHECK_NAN(x->data, x->size);
  73
  74         /* update the workspace from x; */
  75         set_geometry_axes(engine, x->data);
  76
  77         /* kf - ki = Q */
  78         hkl_source_compute_ki(&engine->geometry->source, &ki);
  79         hkl_detector_compute_kf(engine->detector, engine->geometry, &kf);
  80         Q = kf;
  81         hkl_vector_minus_vector(&Q, &ki);
  82         if (hkl_vector_is_null(&Q)){
  83                 f->data[0] = 1;
  84                 f->data[1] = 1;
  85                 f->data[2] = 1;
  86                 f->data[3] = 1;
  87         }else{
  88                 /* R * UB */
  89                 /* for now the 0 holder is the sample holder. */
  90                 sample_holder = darray_item(engine->geometry->holders, 0);
  91                 hkl_quaternion_to_matrix(&sample_holder->q, &RUB);
  92                 hkl_matrix_times_matrix(&RUB, &engine->sample->UB);
  93
  94                 /* compute dhkl0 */
  95                 hkl_matrix_solve(&RUB, &dhkl0, &Q);
  96                 hkl_vector_minus_vector(&dhkl0, &modepsi->hkl0);
  97
  98                 /* compute the intersection of the plan P(kf, ki) and PQ (normal Q) */
  99                 /*
 100                  * now that dhkl0 have been computed we can use a
 101                  * normalized Q to compute n and psi
 102                  */
 103                 hkl_vector_normalize(&Q);
 104                 n = kf;
 105                 hkl_vector_vectorial_product(&n, &ki);
 106                 hkl_vector_vectorial_product(&n, &Q);
 107
 108                 /* compute hkl1 in the laboratory referentiel */
 109                 /* for now the 0 holder is the sample holder. */
 110                 for(unsigned int i=0; i<3; ++i)
 111                         hkl1.data[i] = darray_item(engine->mode->parameters, i)->_value;
 112
 113                 hkl_matrix_times_vector(&engine->sample->UB, &hkl1);
 114                 hkl_vector_rotated_quaternion(&hkl1, &sample_holder->q);
 115
 116                 /* project hkl1 on the plan of normal Q */
 117                 hkl_vector_project_on_plan(&hkl1, &Q);
 118                 if (hkl_vector_is_null(&hkl1)){
 119                         /* hkl1 colinear with Q */
 120                         f->data[0] = dhkl0.data[0];
 121                         f->data[1] = dhkl0.data[1];
 122                         f->data[2] = dhkl0.data[2];
 123                         f->data[3] = 1;
 124                 }else{
 125                         f->data[0] = dhkl0.data[0];
 126                         f->data[1] = dhkl0.data[1];
 127                         f->data[2] = dhkl0.data[2];
 128                         f->data[3] = psi_engine->psi->_value - hkl_vector_oriented_angle(&n, &hkl1, &Q);
 129                 }
 130         }
 131         return GSL_SUCCESS;
 132 }
 133
 134 static int hkl_mode_initialized_set_psi_real(HklMode *self,
 135                                              HklEngine *engine,
 136                                              HklGeometry *geometry,
 137                                              HklDetector *detector,
 138                                              HklSample *sample,
 139                                              int initialized,
 140                                              GError **error)
 141 {
 142         HklVector ki;
 143         HklMatrix RUB;
 144         HklModePsi *psi_mode = container_of(self, HklModePsi, parent);
 145         HklHolder *sample_holder;
 146
 147         hkl_error (error == NULL || *error == NULL);
 148
 149         if(initialized){
 150                 /* update the geometry internals */
 151                 hkl_geometry_update(geometry);
 152
 153                 /* R * UB */
 154                 /* for now the 0 holder is the sample holder. */
 155                 sample_holder = darray_item(geometry->holders, 0);
 156                 hkl_quaternion_to_matrix(&sample_holder->q, &RUB);
 157                 hkl_matrix_times_matrix(&RUB, &sample->UB);
 158
 159                 /* kf - ki = Q0 */
 160                 hkl_source_compute_ki(&geometry->source, &ki);
 161                 hkl_detector_compute_kf(detector, geometry, &psi_mode->Q0);
 162                 hkl_vector_minus_vector(&psi_mode->Q0, &ki);
 163                 if (hkl_vector_is_null(&psi_mode->Q0)){
 164                         g_set_error(error,
 165                                     HKL_MODE_PSI_ERROR,
 166                                     HKL_MODE_PSI_ERROR_INIT,
 167                                     "can not initialize the \"%s\" engine when hkl is null",
 168                                     engine->info->name);
 169                         return FALSE;
 170                 }else
 171                         /* compute hkl0 */
 172                         hkl_matrix_solve(&RUB, &psi_mode->hkl0, &psi_mode->Q0);
 173         }
 174
 175         self->initialized = initialized;
 176
 177         return TRUE;
 178 }
 179
 180 static int hkl_mode_get_psi_real(HklMode *base,
 181                                  HklEngine *engine,
 182                                  HklGeometry *geometry,
 183                                  HklDetector *detector,
 184                                  HklSample *sample,
 185                                  GError **error)
 186 {
 187         HklVector ki;
 188         HklVector kf;
 189         HklVector Q;
 190         HklVector hkl1;
 191         HklVector n;
 192
 193         if (!base || !engine || !engine->mode || !geometry || !detector || !sample){
 194                 g_set_error(error,
 195                             HKL_MODE_PSI_ERROR,
 196                             HKL_MODE_PSI_ERROR_GET,
 197                             "internal error");
 198                 return FALSE;
 199         }
 200
 201         /* get kf, ki and Q */
 202         hkl_source_compute_ki(&geometry->source, &ki);
 203         hkl_detector_compute_kf(detector, geometry, &kf);
 204         Q = kf;
 205         hkl_vector_minus_vector(&Q, &ki);
 206         if (hkl_vector_is_null(&Q)){
 207                 g_set_error(error,
 208                             HKL_MODE_PSI_ERROR,
 209                             HKL_MODE_PSI_ERROR_GET,
 210                             "can not compute psi when hkl is null (kf == ki)");
 211                 return FALSE;
 212         }else{
 213                 /* needed for a problem of precision */
 214                 hkl_vector_normalize(&Q);
 215
 216                 /* compute the intersection of the plan P(kf, ki) and PQ (normal Q) */
 217                 n = kf;
 218                 hkl_vector_vectorial_product(&n, &ki);
 219                 hkl_vector_vectorial_product(&n, &Q);
 220
 221                 /* compute hkl1 in the laboratory referentiel */
 222                 /* the geometry was already updated in the detector compute kf */
 223                 /* for now the 0 holder is the sample holder. */
 224                 for(unsigned int i=0; i<3; ++i)
 225                         hkl1.data[i] = darray_item(base->parameters, i)->_value;
 226
 227                 hkl_matrix_times_vector(&sample->UB, &hkl1);
 228                 hkl_vector_rotated_quaternion(&hkl1, &darray_item(geometry->holders, 0)->q);
 229
 230                 /* project hkl1 on the plan of normal Q */
 231                 hkl_vector_project_on_plan(&hkl1, &Q);
 232
 233                 if (hkl_vector_is_null(&hkl1)){
 234                         g_set_error(error,
 235                                     HKL_MODE_PSI_ERROR,
 236                                     HKL_MODE_PSI_ERROR_GET,
 237                                     "can not compute psi when Q and the ref vector are colinear");
 238                         return FALSE;
 239                 }else{
 240                         HklEnginePsi *psi_engine = container_of(engine, HklEnginePsi, engine);
 241
 242                         /* compute the angle beetween hkl1 and n */
 243                         psi_engine->psi->_value = hkl_vector_oriented_angle(&n, &hkl1, &Q);
 244                 }
 245         }
 246
 247         return TRUE;
 248 }
 249
 250 HklMode *hkl_mode_psi_new(const HklModeAutoInfo *auto_info)
 251 {
 252         static const HklModeOperations operations = {
 253                 HKL_MODE_OPERATIONS_AUTO_DEFAULTS,
 254                 .capabilities = HKL_ENGINE_CAPABILITIES_READABLE | HKL_ENGINE_CAPABILITIES_WRITABLE | HKL_ENGINE_CAPABILITIES_INITIALIZABLE,
 255                 .initialized_set = hkl_mode_initialized_set_psi_real,
 256                 .get = hkl_mode_get_psi_real,
 257         };
 258         HklModePsi *self;
 259
 260         if (darray_size(auto_info->info.axes_w) != 4){
 261                 fprintf(stderr, "This generic HklModePsi need exactly 4 axes");
 262                 exit(128);
 263         }
 264
 265         self = HKL_MALLOC(HklModePsi);
 266
 267         /* the base constructor; */
 268         hkl_mode_auto_init(&self->parent,
 269                            auto_info,
 270                            &operations, FALSE);
 271
 272         return &self->parent;
 273 }
 274
 275 /*************/
 276 /* HklEngine */
 277 /*************/
 278
 279 static void hkl_engine_psi_free_real(HklEngine *base)
 280 {
 281         HklEnginePsi *self=container_of(base, HklEnginePsi, engine);
 282         hkl_engine_release(&self->engine);
 283         free(self);
 284 }
 285
 286 HklEngine *hkl_engine_psi_new(HklEngineList *engines)
 287 {
 288         HklEnginePsi *self;
 289         static const HklParameter psi = {
 290                 HKL_PARAMETER_DEFAULTS_ANGLE, .name = "psi",
 291                 .description = "angle between the reference vector and the diffraction plan",
 292         };
 293         static const HklParameter *pseudo_axes[] = {&psi};
 294         static const HklEngineInfo info = {
 295                 HKL_ENGINE_INFO("psi",
 296                                 pseudo_axes,
 297                                 HKL_ENGINE_DEPENDENCIES_AXES | HKL_ENGINE_DEPENDENCIES_ENERGY | HKL_ENGINE_DEPENDENCIES_SAMPLE),
 298         };
 299         static const HklEngineOperations operations = {
 300                 HKL_ENGINE_OPERATIONS_DEFAULTS,
 301                 .free=hkl_engine_psi_free_real,
 302         };
 303
 304         self = HKL_MALLOC(HklEnginePsi);
 305
 306         hkl_engine_init(&self->engine, &info, &operations, engines);
 307
 308         self->psi = register_pseudo_axis(&self->engine, engines, &psi);
 309
 310         return &self->engine;
 311 }