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-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  */
  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_func(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_matrix(&holder->q, &RUB);
  71                 hkl_matrix_times_matrix(&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_matrix(&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)){
  98                         /* hkl1 colinear with Q */
  99                         f_data[0] = dhkl0.data[0];
 100                         f_data[1] = dhkl0.data[1];
 101                         f_data[2] = dhkl0.data[2];
 102                         f_data[3] = 1;
 103                 }else{
 104                         f_data[0] = dhkl0.data[0];
 105                         f_data[1] = dhkl0.data[1];
 106                         f_data[2] = dhkl0.data[2];
 107                         f_data[3] = psi->parent.value - 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         int status = HKL_SUCCESS;
 121         HklVector ki;
 122         HklMatrix RUB;
 123         HklPseudoAxisEngineModePsi *self = (HklPseudoAxisEngineModePsi *)base;
 124         HklHolder *holder;
 125
 126         status = hkl_pseudo_axis_engine_init_func(base, engine, geometry, detector, sample);
 127         if (status == HKL_FAIL)
 128                 return status;
 129
 130         /* update the geometry internals */
 131         hkl_geometry_update(geometry);
 132
 133         /* R * UB */
 134         /* for now the 0 holder is the sample holder. */
 135         holder = &geometry->holders[0];
 136         hkl_quaternion_to_matrix(&holder->q, &RUB);
 137         hkl_matrix_times_matrix(&RUB, &sample->UB);
 138
 139         /* kf - ki = Q0 */
 140         hkl_source_compute_ki(&geometry->source, &ki);
 141         hkl_detector_compute_kf(detector, geometry, &self->Q0);
 142         hkl_vector_minus_vector(&self->Q0, &ki);
 143         if (hkl_vector_is_null(&self->Q0))
 144                 status = HKL_FAIL;
 145         else
 146                 /* compute hkl0 */
 147                 hkl_matrix_solve(&RUB, &self->hkl0, &self->Q0);
 148
 149         return status;
 150 }
 151
 152 static int hkl_pseudo_axis_engine_mode_get_psi_real(HklPseudoAxisEngineMode *base,
 153                                                     HklPseudoAxisEngine *engine,
 154                                                     HklGeometry *geometry,
 155                                                     HklDetector *detector,
 156                                                     HklSample *sample,
 157                                                     HklError **error)
 158 {
 159         int status = HKL_SUCCESS;
 160         HklVector ki;
 161         HklVector kf;
 162         HklVector Q;
 163         HklVector hkl1;
 164         HklVector n;
 165
 166         if (!base || !engine || !engine->mode || !geometry || !detector || !sample){
 167                 status = HKL_FAIL;
 168                 return status;
 169         }
 170
 171         /* get kf, ki and Q */
 172         hkl_source_compute_ki(&geometry->source, &ki);
 173         hkl_detector_compute_kf(detector, geometry, &kf);
 174         Q = kf;
 175         hkl_vector_minus_vector(&Q, &ki);
 176         if (hkl_vector_is_null(&Q))
 177                 status = HKL_FAIL;
 178         else{
 179                 /* needed for a problem of precision */
 180                 hkl_vector_normalize(&Q);
 181
 182                 /* compute the intersection of the plan P(kf, ki) and PQ (normal Q) */
 183                 n = kf;
 184                 hkl_vector_vectorial_product(&n, &ki);
 185                 hkl_vector_vectorial_product(&n, &Q);
 186
 187                 /* compute hkl1 in the laboratory referentiel */
 188                 /* the geometry was already updated in the detector compute kf */
 189                 /* for now the 0 holder is the sample holder. */
 190                 hkl1.data[0] = base->parameters[0].value;
 191                 hkl1.data[1] = base->parameters[1].value;
 192                 hkl1.data[2] = base->parameters[2].value;
 193                 hkl_vector_times_matrix(&hkl1, &sample->UB);
 194                 hkl_vector_rotated_quaternion(&hkl1, &geometry->holders[0].q);
 195
 196                 /* project hkl1 on the plan of normal Q */
 197                 hkl_vector_project_on_plan(&hkl1, &Q);
 198
 199                 if (hkl_vector_is_null(&hkl1))
 200                         status = HKL_FAIL;
 201                 else
 202                         /* compute the angle beetween hkl1 and n */
 203                         ((HklParameter *)engine->pseudoAxes[0])->value = hkl_vector_oriented_angle(&n, &hkl1, &Q);
 204         }
 205
 206         return status;
 207 }
 208
 209 HklPseudoAxisEngineModePsi *hkl_pseudo_axis_engine_mode_psi_new(char const *name,
 210                                                                 size_t axes_names_len,
 211                                                                 char const *axes_names[])
 212 {
 213         HklPseudoAxisEngineModePsi *self;
 214         HklParameter parameters[3];
 215         HklFunction functions[] = {psi_func};
 216
 217         if (axes_names_len != 4){
 218                 fprintf(stderr, "This generic HklPseudoAxisEngineModePsi need exactly 4 axes");
 219                 exit(128);
 220         }
 221
 222         self = HKL_MALLOC(HklPseudoAxisEngineModePsi);
 223
 224         /* h1  */
 225         hkl_parameter_init(&parameters[0],
 226                            "h1",
 227                            -1, 1, 1,
 228                            HKL_TRUE, HKL_FALSE,
 229                            NULL, NULL);
 230
 231         /* k1 */
 232         hkl_parameter_init(&parameters[1],
 233                            "k1",
 234                            -1, 0, 1,
 235                            HKL_TRUE, HKL_FALSE,
 236                            NULL, NULL);
 237
 238         /* l1 */
 239         hkl_parameter_init(&parameters[2],
 240                            "l1",
 241                            -1, 0, 1,
 242                            HKL_TRUE, HKL_FALSE,
 243                            NULL, NULL);
 244
 245         /* the base constructor; */
 246         hkl_pseudo_axis_engine_mode_init(&self->parent,
 247                                          name,
 248                                          hkl_pseudo_axis_engine_mode_init_psi_real,
 249                                          hkl_pseudo_axis_engine_mode_get_psi_real,
 250                                          hkl_pseudo_axis_engine_mode_set_real,
 251                                          1, functions,
 252                                          3, parameters,
 253                                          axes_names_len, axes_names);
 254
 255
 256         return self;
 257 }
 258
 259 HklPseudoAxisEngine *hkl_pseudo_axis_engine_psi_new(void)
 260 {
 261         HklPseudoAxisEngine *self;
 262
 263         self = hkl_pseudo_axis_engine_new("psi", 1, "psi");
 264
 265         /* psi */
 266         hkl_parameter_init((HklParameter *)self->pseudoAxes[0],
 267                            "psi",
 268                            -M_PI, 0., M_PI,
 269                            HKL_TRUE, HKL_TRUE,
 270                            &hkl_unit_angle_rad, &hkl_unit_angle_deg);
 271
 272         return self;
 273 }