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-2012 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_math.h>
  23 #include <gsl/gsl_vector.h>
  24 #include <gsl/gsl_sf_trig.h>
  25
  26 #include <ccan/array_size/array_size.h>
  27
  28 #include "hkl-pseudoaxis-auto-private.h"
  29 #include "hkl-pseudoaxis-common-psi-private.h"
  30 #include "hkl-pseudoaxis-common-private.h"
  31
  32
  33 /***********************/
  34 /* numerical functions */
  35 /***********************/
  36
  37 int _psi_func(const gsl_vector *x, void *params, gsl_vector *f)
  38 {
  39
  40         HklVector dhkl0, hkl1;
  41         HklVector ki, kf, Q, n;
  42         HklMatrix RUB;
  43         HklPseudoAxisEngine *engine = params;
  44         HklPseudoAxisEngineModePsi *modepsi = container_of(engine->mode, HklPseudoAxisEngineModePsi, parent);
  45         HklHolder *holder;
  46
  47         CHECK_NAN(x->data, x->size);
  48
  49         /* update the workspace from x; */
  50         set_geometry_axes(engine, x->data);
  51
  52         /* kf - ki = Q */
  53         hkl_source_compute_ki(&engine->geometry->source, &ki);
  54         hkl_detector_compute_kf(engine->detector, engine->geometry, &kf);
  55         Q = kf;
  56         hkl_vector_minus_vector(&Q, &ki);
  57         if (hkl_vector_is_null(&Q)){
  58                 f->data[0] = 1;
  59                 f->data[1] = 1;
  60                 f->data[2] = 1;
  61                 f->data[3] = 1;
  62         }else{
  63                 /* R * UB */
  64                 /* for now the 0 holder is the sample holder. */
  65                 holder = &engine->geometry->holders[0];
  66                 hkl_quaternion_to_matrix(&holder->q, &RUB);
  67                 hkl_matrix_times_matrix(&RUB, &engine->sample->UB);
  68
  69                 /* compute dhkl0 */
  70                 hkl_matrix_solve(&RUB, &dhkl0, &Q);
  71                 hkl_vector_minus_vector(&dhkl0, &modepsi->hkl0);
  72
  73                 /* compute the intersection of the plan P(kf, ki) and PQ (normal Q) */
  74                 /*
  75                  * now that dhkl0 have been computed we can use a
  76                  * normalized Q to compute n and psi
  77                  */
  78                 hkl_vector_normalize(&Q);
  79                 n = kf;
  80                 hkl_vector_vectorial_product(&n, &ki);
  81                 hkl_vector_vectorial_product(&n, &Q);
  82
  83                 /* compute hkl1 in the laboratory referentiel */
  84                 /* for now the 0 holder is the sample holder. */
  85                 hkl1.data[0] = engine->mode->parameters[0].value;
  86                 hkl1.data[1] = engine->mode->parameters[1].value;
  87                 hkl1.data[2] = engine->mode->parameters[2].value;
  88                 hkl_matrix_times_vector(&engine->sample->UB, &hkl1);
  89                 hkl_vector_rotated_quaternion(&hkl1, &engine->geometry->holders[0].q);
  90
  91                 /* project hkl1 on the plan of normal Q */
  92                 hkl_vector_project_on_plan(&hkl1, &Q);
  93                 if (hkl_vector_is_null(&hkl1)){
  94                         /* hkl1 colinear with Q */
  95                         f->data[0] = dhkl0.data[0];
  96                         f->data[1] = dhkl0.data[1];
  97                         f->data[2] = dhkl0.data[2];
  98                         f->data[3] = 1;
  99                 }else{
 100                         double psi;
 101                         uint len = 1;
 102
 103                         f->data[0] = dhkl0.data[0];
 104                         f->data[1] = dhkl0.data[1];
 105                         f->data[2] = dhkl0.data[2];
 106                         hkl_pseudo_axis_engine_get_values(engine, &psi, &len);
 107                         f->data[3] = psi - hkl_vector_oriented_angle(&n, &hkl1, &Q);
 108                 }
 109         }
 110         return GSL_SUCCESS;
 111 }
 112
 113 static int hkl_pseudo_axis_engine_mode_init_psi_real(HklPseudoAxisEngineMode *base,
 114                                                      HklPseudoAxisEngine *engine,
 115                                                      HklGeometry *geometry,
 116                                                      HklDetector *detector,
 117                                                      HklSample *sample,
 118                                                      HklError **error)
 119 {
 120         HklVector ki;
 121         HklMatrix RUB;
 122         HklPseudoAxisEngineModePsi *self = container_of(base, HklPseudoAxisEngineModePsi, parent);
 123         HklHolder *holder;
 124
 125         hkl_return_val_if_fail (error == NULL || *error == NULL, HKL_FALSE);
 126
 127         if (!hkl_pseudo_axis_engine_init_func(base, engine, geometry, detector, sample)){
 128                 hkl_error_set(error, "internal error");
 129                 return HKL_FALSE;
 130         }
 131         hkl_assert(error == NULL || *error == NULL);
 132
 133         /* update the geometry internals */
 134         hkl_geometry_update(geometry);
 135
 136         /* R * UB */
 137         /* for now the 0 holder is the sample holder. */
 138         holder = &geometry->holders[0];
 139         hkl_quaternion_to_matrix(&holder->q, &RUB);
 140         hkl_matrix_times_matrix(&RUB, &sample->UB);
 141
 142         /* kf - ki = Q0 */
 143         hkl_source_compute_ki(&geometry->source, &ki);
 144         hkl_detector_compute_kf(detector, geometry, &self->Q0);
 145         hkl_vector_minus_vector(&self->Q0, &ki);
 146         if (hkl_vector_is_null(&self->Q0)){
 147                 hkl_error_set(error, "can not initialize the \"%s\" engine when hkl is null",
 148                               engine->info->name);
 149                 return HKL_FALSE;
 150         }else
 151                 /* compute hkl0 */
 152                 hkl_matrix_solve(&RUB, &self->hkl0, &self->Q0);
 153
 154         return HKL_TRUE;
 155 }
 156
 157 static int hkl_pseudo_axis_engine_mode_get_psi_real(HklPseudoAxisEngineMode *base,
 158                                                     HklPseudoAxisEngine *engine,
 159                                                     HklGeometry *geometry,
 160                                                     HklDetector *detector,
 161                                                     HklSample *sample,
 162                                                     HklError **error)
 163 {
 164         HklVector ki;
 165         HklVector kf;
 166         HklVector Q;
 167         HklVector hkl1;
 168         HklVector n;
 169
 170         if (!base || !engine || !engine->mode || !geometry || !detector || !sample){
 171                 hkl_error_set(error, "internal error");
 172                 return HKL_FALSE;
 173         }
 174
 175         /* get kf, ki and Q */
 176         hkl_source_compute_ki(&geometry->source, &ki);
 177         hkl_detector_compute_kf(detector, geometry, &kf);
 178         Q = kf;
 179         hkl_vector_minus_vector(&Q, &ki);
 180         if (hkl_vector_is_null(&Q)){
 181                 hkl_error_set(error, "can not compute psi when hkl is null (kf == ki)");
 182                 return HKL_FALSE;
 183         }else{
 184                 /* needed for a problem of precision */
 185                 hkl_vector_normalize(&Q);
 186
 187                 /* compute the intersection of the plan P(kf, ki) and PQ (normal Q) */
 188                 n = kf;
 189                 hkl_vector_vectorial_product(&n, &ki);
 190                 hkl_vector_vectorial_product(&n, &Q);
 191
 192                 /* compute hkl1 in the laboratory referentiel */
 193                 /* the geometry was already updated in the detector compute kf */
 194                 /* for now the 0 holder is the sample holder. */
 195                 hkl1.data[0] = base->parameters[0].value;
 196                 hkl1.data[1] = base->parameters[1].value;
 197                 hkl1.data[2] = base->parameters[2].value;
 198                 hkl_matrix_times_vector(&sample->UB, &hkl1);
 199                 hkl_vector_rotated_quaternion(&hkl1, &geometry->holders[0].q);
 200
 201                 /* project hkl1 on the plan of normal Q */
 202                 hkl_vector_project_on_plan(&hkl1, &Q);
 203
 204                 if (hkl_vector_is_null(&hkl1)){
 205                         hkl_error_set(error, "can not compute psi when Q and the ref vector are colinear");
 206                         return HKL_FALSE;
 207                 }else{
 208                         double psi;
 209
 210                         /* compute the angle beetween hkl1 and n */
 211                         psi = hkl_vector_oriented_angle(&n, &hkl1, &Q);
 212                         hkl_pseudo_axis_engine_set_values(engine, &psi, 1);
 213                 }
 214         }
 215
 216         return HKL_TRUE;
 217 }
 218
 219 static const HklPseudoAxisEngineModeOperations psi_mode_operations = {
 220         .init = hkl_pseudo_axis_engine_mode_init_psi_real,
 221         .get = hkl_pseudo_axis_engine_mode_get_psi_real,
 222         .set = hkl_pseudo_axis_engine_mode_set_real
 223 };
 224
 225 HklPseudoAxisEngineMode *hkl_pseudo_axis_engine_mode_psi_new(const HklPseudoAxisEngineModeAutoInfo *info)
 226 {
 227         HklPseudoAxisEngineModePsi *self;
 228
 229         if (info->mode.n_axes != 4){
 230                 fprintf(stderr, "This generic HklPseudoAxisEngineModePsi need exactly 4 axes");
 231                 exit(128);
 232         }
 233
 234         self = HKL_MALLOC(HklPseudoAxisEngineModePsi);
 235
 236         /* the base constructor; */
 237         hkl_pseudo_axis_engine_mode_auto_init(&self->parent,
 238                                               info,
 239                                               &psi_mode_operations);
 240
 241         return &self->parent;
 242 }
 243
 244 HklPseudoAxisEngine *hkl_pseudo_axis_engine_psi_new(void)
 245 {
 246         HklPseudoAxisEngine *self;
 247         static const HklPseudoAxis psi = {
 248                 .parent = { HKL_PARAMETER_DEFAULTS_ANGLE, .name = "psi"}
 249         };
 250         static const HklPseudoAxis *pseudo_axes[] = {&psi};
 251         static const HklPseudoAxisEngineInfo info = {
 252                 .name = "psi",
 253                 .pseudo_axes = pseudo_axes,
 254                 .n_pseudo_axes = ARRAY_SIZE(pseudo_axes),
 255         };
 256
 257         self = hkl_pseudo_axis_engine_new(&info);
 258
 259         return self;
 260 }