src/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-2009 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 <hkl/hkl-pseudoaxis.h>
  27 #include <hkl/hkl-pseudoaxis-common.h>
  28 #include <hkl/hkl-pseudoaxis-auto.h>
  29 #include <hkl/hkl-pseudoaxis-common-psi.h>
  30
  31 static int psi(const gsl_vector *x, void *params, gsl_vector *f)
  32 {
  33
  34         HklVector dhkl0, hkl1;
  35         HklVector ki, kf, Q, n;
  36         HklMatrix RUB;
  37         HklPseudoAxisEngine *engine;
  38         HklPseudoAxisEngineModePsi *modepsi;
  39         HklPseudoAxis *psi;
  40         HklHolder *holder;
  41         size_t i;
  42         size_t len;
  43         double const *x_data = gsl_vector_const_ptr(x, 0);
  44         double *f_data = gsl_vector_ptr(f, 0);
  45
  46         engine = params;
  47         modepsi = (HklPseudoAxisEngineModePsi *)engine->mode;
  48         psi = engine->pseudoAxes[0];
  49
  50         // update the workspace from x;
  51         len = HKL_LIST_LEN(engine->axes);
  52         for(i=0; i<len; ++i)
  53                 hkl_axis_set_value(engine->axes[i], x_data[i]);
  54         hkl_geometry_update(engine->geometry);
  55
  56         // kf - ki = Q
  57         hkl_source_compute_ki(&engine->geometry->source, &ki);
  58         hkl_detector_compute_kf(engine->detector, engine->geometry, &kf);
  59         Q = kf;
  60         hkl_vector_minus_vector(&Q, &ki);
  61         if (hkl_vector_is_null(&Q)){
  62                 f_data[0] = 1;
  63                 f_data[1] = 1;
  64                 f_data[2] = 1;
  65                 f_data[3] = 1;
  66         }else{
  67                 // R * UB
  68                 // for now the 0 holder is the sample holder.
  69                 holder = &engine->geometry->holders[0];
  70                 hkl_quaternion_to_smatrix(&holder->q, &RUB);
  71                 hkl_matrix_times_smatrix(&RUB, &engine->sample->UB);
  72
  73                 // compute dhkl0
  74                 hkl_matrix_solve(&RUB, &dhkl0, &Q);
  75                 hkl_vector_minus_vector(&dhkl0, &modepsi->hkl0);
  76
  77                 // compute the intersection of the plan P(kf, ki) and PQ (normal Q)
  78                 /*
  79                  * now that dhkl0 have been computed we can use a
  80                  * normalized Q to compute n and psi
  81                  */
  82                 hkl_vector_normalize(&Q);
  83                 n = kf;
  84                 hkl_vector_vectorial_product(&n, &ki);
  85                 hkl_vector_vectorial_product(&n, &Q);
  86
  87                 // compute hkl1 in the laboratory referentiel
  88                 // for now the 0 holder is the sample holder.
  89                 hkl1.data[0] = engine->mode->parameters[0].value;
  90                 hkl1.data[1] = engine->mode->parameters[1].value;
  91                 hkl1.data[2] = engine->mode->parameters[2].value;
  92                 hkl_vector_times_smatrix(&hkl1, &engine->sample->UB);
  93                 hkl_vector_rotated_quaternion(&hkl1, &engine->geometry->holders[0].q);
  94
  95                 // project hkl1 on the plan of normal Q
  96                 hkl_vector_project_on_plan(&hkl1, &Q);
  97                 if (hkl_vector_is_null(&hkl1)){ // hkl1 colinear with Q
  98                         f_data[0] = dhkl0.data[0];
  99                         f_data[1] = dhkl0.data[1];
 100                         f_data[2] = dhkl0.data[2];
 101                         f_data[3] = 1;
 102                 }else{
 103                         f_data[0] = dhkl0.data[0];
 104                         f_data[1] = dhkl0.data[1];
 105                         f_data[2] = dhkl0.data[2];
 106                         f_data[3] = psi->parent.value - hkl_vector_oriented_angle(&n, &hkl1, &Q);
 107                 }
 108         }
 109         return GSL_SUCCESS;
 110 }
 111
 112 static int hkl_pseudo_axis_engine_mode_init_psi_real(HklPseudoAxisEngine *engine,
 113                                                         HklGeometry *geometry,
 114                                                         HklDetector const *detector,
 115                                                         HklSample const *sample)
 116 {
 117         int status = HKL_SUCCESS;
 118         HklVector ki;
 119         HklMatrix RUB;
 120         HklPseudoAxisEngineModePsi *self;
 121         HklHolder *holder;
 122
 123         status = hkl_pseudo_axis_engine_init_func(engine, geometry, detector, sample);
 124         if (status == HKL_FAIL)
 125                 return status;
 126
 127         self = (HklPseudoAxisEngineModePsi *)engine->mode;
 128
 129         // update the geometry internals
 130         hkl_geometry_update(geometry);
 131
 132         // R * UB
 133         // for now the 0 holder is the sample holder.
 134         holder = &geometry->holders[0];
 135         hkl_quaternion_to_smatrix(&holder->q, &RUB);
 136         hkl_matrix_times_smatrix(&RUB, &sample->UB);
 137
 138         // kf - ki = Q0
 139         hkl_source_compute_ki(&geometry->source, &ki);
 140         hkl_detector_compute_kf(detector, geometry, &self->Q0);
 141         hkl_vector_minus_vector(&self->Q0, &ki);
 142         if (hkl_vector_is_null(&self->Q0))
 143                 status = HKL_FAIL;
 144         else
 145                 // compute hkl0
 146                 hkl_matrix_solve(&RUB, &self->hkl0, &self->Q0);
 147
 148         return status;
 149 }
 150
 151 static int hkl_pseudo_axis_engine_mode_get_psi_real(HklPseudoAxisEngine *engine,
 152                                                        HklGeometry *geometry,
 153                                                        HklDetector const *detector,
 154                                                        HklSample const *sample)
 155 {
 156         int status = HKL_SUCCESS;
 157
 158         if (!engine || !engine->mode || !geometry || !detector || !sample){
 159                 status = HKL_FAIL;
 160                 return status;
 161         }
 162
 163         HklVector ki;
 164         HklVector kf;
 165         HklVector Q;
 166         HklVector hkl1;
 167         HklVector n;
 168         HklPseudoAxisEngineModePsi *self;
 169         HklPseudoAxisEngineMode *base;
 170
 171         self = (HklPseudoAxisEngineModePsi *)engine->mode;
 172         base = engine->mode;
 173
 174         // get kf, ki and Q
 175         hkl_source_compute_ki(&geometry->source, &ki);
 176         hkl_detector_compute_kf(detector, geometry, &kf);
 177         Q = kf;
 178         hkl_vector_minus_vector(&Q, &ki);
 179         if (hkl_vector_is_null(&Q))
 180                 status = HKL_FAIL;
 181         else{
 182                 hkl_vector_normalize(&Q); // needed for a problem of precision
 183
 184                 // compute the intersection of the plan P(kf, ki) and PQ (normal Q)
 185                 n = kf;
 186                 hkl_vector_vectorial_product(&n, &ki);
 187                 hkl_vector_vectorial_product(&n, &Q);
 188
 189                 // compute hkl1 in the laboratory referentiel
 190                 // the geometry was already updated in the detector compute kf
 191                 // for now the 0 holder is the sample holder.
 192                 hkl1.data[0] = base->parameters[0].value;
 193                 hkl1.data[1] = base->parameters[1].value;
 194                 hkl1.data[2] = base->parameters[2].value;
 195                 hkl_vector_times_smatrix(&hkl1, &sample->UB);
 196                 hkl_vector_rotated_quaternion(&hkl1, &geometry->holders[0].q);
 197
 198                 // project hkl1 on the plan of normal Q
 199                 hkl_vector_project_on_plan(&hkl1, &Q);
 200
 201                 if (hkl_vector_is_null(&hkl1))
 202                         status = HKL_FAIL;
 203                 else
 204                         // compute the angle beetween hkl1 and n
 205                         ((HklParameter *)engine->pseudoAxes[0])->value = hkl_vector_oriented_angle(&n, &hkl1, &Q);
 206         }
 207
 208         return status;
 209 }
 210
 211 static int hkl_pseudo_axis_engine_mode_set_psi_real(HklPseudoAxisEngine *engine,
 212                                                        HklGeometry *geometry,
 213                                                        HklDetector *detector,
 214                                                        HklSample *sample)
 215 {
 216         hkl_pseudo_axis_engine_prepare_internal(engine, geometry, detector,
 217                                                 sample);
 218
 219         return hkl_pseudo_axis_engine_solve_function(engine, psi);
 220 }
 221
 222 HklPseudoAxisEngineModePsi *hkl_pseudo_axis_engine_mode_psi_new(char const *name,
 223                                                                      size_t axes_names_len,
 224                                                                      char const *axes_names[])
 225 {
 226         HklPseudoAxisEngineModePsi *self;
 227         char const *parameters_names[] = {"h1", "k1", "l1"};
 228
 229         if (axes_names_len != 4)
 230                 die("This generic HklPseudoAxisEngineModePsi need exactly 4 axes");
 231
 232         self = calloc(1, sizeof(*self));
 233         if (!self)
 234                 die("Can not allocate memory for an HklPseudoAxisEngineModePsi");
 235
 236         // the base constructor;
 237         hkl_pseudo_axis_engine_mode_init(&self->parent,
 238                                             name,
 239                                             hkl_pseudo_axis_engine_mode_init_psi_real,
 240                                             hkl_pseudo_axis_engine_mode_get_psi_real,
 241                                             hkl_pseudo_axis_engine_mode_set_psi_real,
 242                                             3, parameters_names,
 243                                             axes_names_len, axes_names);
 244
 245         self->parent.parameters[0].value = 1;
 246         self->parent.parameters[0].range.min = -1;
 247         self->parent.parameters[0].range.max = 1;
 248         self->parent.parameters[0].not_to_fit = HKL_FALSE;
 249
 250         self->parent.parameters[1].value = 0;
 251         self->parent.parameters[1].range.min = -1;
 252         self->parent.parameters[1].range.max = 1;
 253         self->parent.parameters[1].not_to_fit = HKL_FALSE;
 254
 255         self->parent.parameters[2].value = 0;
 256         self->parent.parameters[2].range.min = -1;
 257         self->parent.parameters[2].range.max = 1;
 258         self->parent.parameters[2].not_to_fit = HKL_FALSE;
 259
 260         return self;
 261 }
 262
 263 HklPseudoAxisEngine *hkl_pseudo_axis_engine_psi_new(void)
 264 {
 265         HklPseudoAxisEngine *self;
 266
 267         self = hkl_pseudo_axis_engine_new("psi", 1, "psi");
 268
 269         // psi
 270         hkl_parameter_init((HklParameter *)self->pseudoAxes[0],
 271                            "psi",
 272                            -M_PI, 0., M_PI,
 273                            HKL_FALSE, HKL_TRUE,
 274                            &hkl_unit_angle_rad, &hkl_unit_angle_deg);
 275
 276         return self;
 277 }