hkl/hkl-pseudoaxis-common-hkl.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-2010 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  *          Maria-Teresa Nunez-Pardo-de-Verra <tnunez@mail.desy.de>
  22  */
  23 #include <string.h>
  24 #include <gsl/gsl_sf_trig.h>
  25 #include <hkl/hkl-pseudoaxis.h>
  26 #include <hkl/hkl-pseudoaxis-common.h>
  27 #include <hkl/hkl-pseudoaxis-common-hkl.h>
  28 #include <hkl/hkl-pseudoaxis-auto.h>
  29 /***************************************/
  30 /* common methode use by getter/setter */
  31 /***************************************/
  32
  33 int RUBh_minus_Q_func(const gsl_vector *x, void *params, gsl_vector *f)
  34 {
  35         double const *x_data = gsl_vector_const_ptr(x, 0);
  36         double *f_data = gsl_vector_ptr(f, 0);
  37
  38         RUBh_minus_Q(x_data, params, f_data);
  39
  40         return  GSL_SUCCESS;
  41 }
  42
  43 int RUBh_minus_Q(double const x[], void *params, double f[])
  44 {
  45         HklVector Hkl;
  46         HklVector ki, dQ;
  47         HklPseudoAxisEngine *engine;
  48         HklHolder *holder;
  49         size_t i;
  50
  51         engine = params;
  52
  53         // update the workspace from x;
  54         for(i=0; i<HKL_LIST_LEN(engine->axes); ++i)
  55                 hkl_axis_set_value(engine->axes[i], x[i]);
  56         hkl_geometry_update(engine->geometry);
  57
  58         hkl_vector_init(&Hkl,
  59                         ((HklParameter *)engine->pseudoAxes[0])->value,
  60                         ((HklParameter *)engine->pseudoAxes[1])->value,
  61                         ((HklParameter *)engine->pseudoAxes[2])->value);
  62
  63         // R * UB * h = Q
  64         // for now the 0 holder is the sample holder.
  65         holder = &engine->geometry->holders[0];
  66         hkl_matrix_times_vector(&engine->sample->UB, &Hkl);
  67         hkl_vector_rotated_quaternion(&Hkl, &holder->q);
  68
  69         // kf - ki = Q
  70         hkl_source_compute_ki(&engine->geometry->source, &ki);
  71         hkl_detector_compute_kf(engine->detector, engine->geometry, &dQ);
  72         hkl_vector_minus_vector(&dQ, &ki);
  73
  74         hkl_vector_minus_vector(&dQ, &Hkl);
  75
  76         f[0] = dQ.data[0];
  77         f[1] = dQ.data[1];
  78         f[2] = dQ.data[2];
  79
  80         return GSL_SUCCESS;
  81 }
  82
  83 int hkl_pseudo_axis_engine_mode_get_hkl_real(HklPseudoAxisEngineMode *self,
  84                                              HklPseudoAxisEngine *engine,
  85                                              HklGeometry *geometry,
  86                                              HklDetector *detector,
  87                                              HklSample *sample,
  88                                              HklError **error)
  89 {
  90         HklHolder *holder;
  91         HklMatrix RUB;
  92         HklVector hkl, ki, Q;
  93         double min, max;
  94         size_t i;
  95
  96         // update the geometry internals
  97         hkl_geometry_update(geometry);
  98
  99         // R * UB
 100         // for now the 0 holder is the sample holder.
 101         holder = &geometry->holders[0];
 102         hkl_quaternion_to_matrix(&holder->q, &RUB);
 103         hkl_matrix_times_matrix(&RUB, &sample->UB);
 104
 105         // kf - ki = Q
 106         hkl_source_compute_ki(&geometry->source, &ki);
 107         hkl_detector_compute_kf(detector, geometry, &Q);
 108         hkl_vector_minus_vector(&Q, &ki);
 109
 110         hkl_matrix_solve(&RUB, &hkl, &Q);
 111
 112         // compute the min and max
 113         min = -1;
 114         max = 1;
 115
 116         // update the pseudoAxes config part
 117         for(i=0;i<HKL_LIST_LEN(engine->pseudoAxes);++i){
 118                 HklParameter *parameter = (HklParameter *)(engine->pseudoAxes[i]);
 119                 parameter->value = hkl.data[i];
 120                 parameter->range.min = min;
 121                 parameter->range.max = max;
 122         }
 123
 124         return HKL_SUCCESS;
 125 }
 126
 127 /***************************************/
 128 /* the double diffraction get set part */
 129 /***************************************/
 130
 131 int double_diffraction_func(gsl_vector const *x, void *params, gsl_vector *f)
 132 {
 133         double const *x_data = gsl_vector_const_ptr(x, 0);
 134         double *f_data = gsl_vector_ptr(f, 0);
 135
 136         double_diffraction(x_data, params, f_data);
 137
 138         return  GSL_SUCCESS;
 139 }
 140
 141 int double_diffraction(double const x[], void *params, double f[])
 142 {
 143         HklPseudoAxisEngine *engine = params;
 144         HklVector hkl, kf2;
 145         HklVector ki;
 146         HklVector dQ;
 147         size_t i;
 148         HklHolder *holder;
 149
 150         // update the workspace from x;
 151         for(i=0; i<HKL_LIST_LEN(engine->axes); ++i)
 152                 hkl_axis_set_value(engine->axes[i], x[i]);
 153         hkl_geometry_update(engine->geometry);
 154
 155         hkl_vector_init(&hkl,
 156                         ((HklParameter *)engine->pseudoAxes[0])->value,
 157                         ((HklParameter *)engine->pseudoAxes[1])->value,
 158                         ((HklParameter *)engine->pseudoAxes[2])->value);
 159
 160         hkl_vector_init(&kf2,
 161                         engine->mode->parameters[0].value,
 162                         engine->mode->parameters[1].value,
 163                         engine->mode->parameters[2].value);
 164
 165         // R * UB * hkl = Q
 166         // for now the 0 holder is the sample holder.
 167         holder = &engine->geometry->holders[0];
 168         hkl_matrix_times_vector(&engine->sample->UB, &hkl);
 169         hkl_vector_rotated_quaternion(&hkl, &holder->q);
 170
 171         // kf - ki = Q
 172         hkl_source_compute_ki(&engine->geometry->source, &ki);
 173         hkl_detector_compute_kf(engine->detector, engine->geometry, &dQ);
 174         hkl_vector_minus_vector(&dQ, &ki);
 175         hkl_vector_minus_vector(&dQ, &hkl);
 176
 177         // R * UB * hlk2 = Q2
 178         hkl_matrix_times_vector(&engine->sample->UB, &kf2);
 179         hkl_vector_rotated_quaternion(&kf2, &holder->q);
 180         hkl_vector_add_vector(&kf2, &ki);
 181
 182         f[0] = dQ.data[0];
 183         f[1] = dQ.data[1];
 184         f[2] = dQ.data[2];
 185         f[3] = hkl_vector_norm2(&kf2) - hkl_vector_norm2(&ki);
 186
 187         return GSL_SUCCESS;
 188 }
 189
 190 /******************************************/
 191 /* the psi_constant_vertical get set part */
 192 /******************************************/
 193
 194 int psi_constant_vertical_func(gsl_vector const *x, void *params, gsl_vector *f)
 195 {
 196
 197         double const *x_data = gsl_vector_const_ptr(x, 0);
 198         double *f_data = gsl_vector_ptr(f, 0);
 199
 200         HklVector hkl;
 201         HklVector ki, kf, Q, n;
 202         HklPseudoAxisEngine *engine;
 203         size_t i;
 204
 205         RUBh_minus_Q(x_data, params, f_data);
 206         engine = params;
 207
 208         // update the workspace from x;
 209         for(i=0; i<HKL_LIST_LEN(engine->axes); ++i)
 210                 hkl_axis_set_value(engine->axes[i], x_data[i]);
 211         hkl_geometry_update(engine->geometry);
 212
 213         hkl_vector_init(&hkl, 1, 0, 0);
 214
 215         // kf - ki = Q
 216         hkl_source_compute_ki(&engine->geometry->source, &ki);
 217         hkl_detector_compute_kf(engine->detector, engine->geometry, &kf);
 218         Q = kf;
 219         hkl_vector_minus_vector(&Q, &ki);
 220
 221         hkl_vector_normalize(&Q);
 222         n = kf;
 223         hkl_vector_vectorial_product(&n, &ki);
 224         hkl_vector_vectorial_product(&n, &Q);
 225
 226
 227         hkl_vector_times_smatrix(&hkl, &engine->sample->UB);
 228         hkl_vector_rotated_quaternion(&hkl, &engine->geometry->holders[0].q);
 229
 230         // project hkl on the plan of normal Q
 231         hkl_vector_project_on_plan(&hkl, &Q);
 232
 233         f_data[3] =  engine->mode->parameters[3].value - hkl_vector_oriented_angle(&n, &hkl, &Q);
 234
 235         return  GSL_SUCCESS;
 236 }
 237
 238 int hkl_pseudo_axis_engine_mode_init_psi_constant_vertical_real(HklPseudoAxisEngineMode *self,
 239                                                                 HklPseudoAxisEngine *engine,
 240                                                                 HklGeometry *geometry,
 241                                                                 HklDetector *detector,
 242                                                                 HklSample *sample,
 243                                                                 HklError **error)
 244 {
 245         HklVector hkl;
 246         HklVector ki, kf, Q, n;
 247         int status = HKL_SUCCESS;
 248
 249         if (!self || !engine || !engine->mode || !geometry || !detector || !sample){
 250                 status = HKL_FAIL;
 251                 return status;
 252         }
 253
 254         status = hkl_pseudo_axis_engine_init_func(self, engine, geometry, detector, sample);
 255         if(status == HKL_FAIL)
 256                 return status;
 257
 258         // Compute the constant psi value (to be kept)
 259         hkl_vector_init(&hkl, 1, 0, 0);
 260
 261         // kf - ki = Q
 262         hkl_source_compute_ki(&geometry->source, &ki);
 263         hkl_detector_compute_kf(detector, geometry, &kf);
 264         Q = kf;
 265         hkl_vector_minus_vector(&Q, &ki);
 266
 267         if (hkl_vector_is_null(&Q))
 268                 status = HKL_FAIL;
 269         else{
 270                 hkl_vector_normalize(&Q); // needed for a problem of precision
 271
 272                 // compute the intersection of the plan P(kf, ki) and PQ (normal Q)
 273                 n = kf;
 274                 hkl_vector_vectorial_product(&n, &ki);
 275                 hkl_vector_vectorial_product(&n, &Q);
 276
 277                 // compute hkl in the laboratory referentiel
 278                 // the geometry was already updated in the detector compute kf
 279                 // for now the 0 holder is the sample holder
 280                 hkl.data[0] = self->parameters[0].value;
 281                 hkl.data[1] = self->parameters[1].value;
 282                 hkl.data[2] = self->parameters[2].value;
 283                 hkl_vector_times_smatrix(&hkl, &sample->UB);
 284                 hkl_vector_rotated_quaternion(&hkl, &geometry->holders[0].q);
 285
 286                 // project hkl on the plan of normal Q
 287                 hkl_vector_project_on_plan(&hkl, &Q);
 288
 289                 if (hkl_vector_is_null(&hkl))
 290                         status = HKL_FAIL;
 291                 else
 292                         // compute the angle beetween hkl and n and store in in the fourth parameter
 293                         hkl_parameter_set_value(&self->parameters[3],
 294                                                 hkl_vector_oriented_angle(&n, &hkl, &Q));
 295         }
 296
 297         return status;
 298 }
 299
 300
 301 HklPseudoAxisEngine *hkl_pseudo_axis_engine_hkl_new(void)
 302 {
 303         HklPseudoAxisEngine *self;
 304
 305         self = hkl_pseudo_axis_engine_new("hkl", 3, "h", "k", "l");
 306
 307         // h
 308         hkl_parameter_init((HklParameter *)self->pseudoAxes[0],
 309                            "h",
 310                            -1, 0., 1,
 311                            HKL_TRUE, HKL_TRUE,
 312                            NULL, NULL);
 313         // k
 314         hkl_parameter_init((HklParameter *)self->pseudoAxes[1],
 315                            "k",
 316                            -1, 0., 1,
 317                            HKL_TRUE, HKL_TRUE,
 318                            NULL, NULL);
 319         // l
 320         hkl_parameter_init((HklParameter *)self->pseudoAxes[2],
 321                            "l",
 322                            -1, 0., 1,
 323                            HKL_TRUE, HKL_TRUE,
 324                            NULL, NULL);
 325
 326         return self;
 327 }