src/hkl-pseudoaxis-common.c

   1 #include <string.h>
   2 #include <gsl/gsl_sf_trig.h>
   3 #include <hkl/hkl-pseudoaxis-common.h>
   4 #include <hkl/hkl-pseudoaxis-auto.h>
   5
   6 /***************************************/
   7 /* common methode use by getter/setter */
   8 /***************************************/
   9
  10 static int RUBh_minus_Q_func(const gsl_vector *x, void *params, gsl_vector *f)
  11 {
  12         double const *x_data = gsl_vector_const_ptr(x, 0);
  13         double *f_data = gsl_vector_ptr(f, 0);
  14
  15         RUBh_minus_Q(x_data, params, f_data);
  16
  17         return  GSL_SUCCESS;
  18 }
  19
  20 int RUBh_minus_Q(double const x[], void *params, double f[])
  21 {
  22         HklVector Hkl;
  23         HklVector ki, dQ;
  24         HklPseudoAxisEngine *engine;
  25         HklPseudoAxis *H, *K, *L;
  26         HklHolder *holder;
  27         size_t i;
  28
  29         engine = params;
  30         H = &engine->pseudoAxes[0];
  31         K = &engine->pseudoAxes[1];
  32         L = &engine->pseudoAxes[2];
  33
  34         // update the workspace from x;
  35         for(i=0; i<engine->axes_len; ++i) {
  36                 HklAxis *axis = engine->axes[i];
  37                 axis->config.value = x[i];
  38                 axis->config.dirty = 1;
  39         }
  40         hkl_geometry_update(engine->geometry);
  41
  42         hkl_vector_init(&Hkl, H->config.value, K->config.value,
  43                         L->config.value);
  44
  45         // R * UB * h = Q
  46         // for now the 0 holder is the sample holder.
  47         holder = &engine->geometry->holders[0];
  48         hkl_matrix_times_vector(&engine->sample->UB, &Hkl);
  49         hkl_vector_rotated_quaternion(&Hkl, &holder->q);
  50
  51         // kf - ki = Q
  52         hkl_source_compute_ki(&engine->geometry->source, &ki);
  53         hkl_detector_compute_kf(engine->detector, engine->geometry, &dQ);
  54         hkl_vector_minus_vector(&dQ, &ki);
  55
  56         hkl_vector_minus_vector(&dQ, &Hkl);
  57
  58         f[0] = dQ.data[0];
  59         f[1] = dQ.data[1];
  60         f[2] = dQ.data[2];
  61
  62         return GSL_SUCCESS;
  63 }
  64
  65 int hkl_pseudo_axis_engine_init_func(HklPseudoAxisEngine *self,
  66                                      HklGeometry *geometry,
  67                                      HklDetector *detector,
  68                                      HklSample *sample)
  69 {
  70         HklPseudoAxisEngineGetSet *getset;
  71
  72         if (!self || !self->getset || !geometry || !detector || !sample)
  73                 return HKL_FAIL;
  74
  75         getset = self->getset;
  76
  77         // update the geometry internals
  78         hkl_geometry_update(geometry);
  79
  80         if(getset->geometry_init)
  81                 hkl_geometry_free(getset->geometry_init);
  82         getset->geometry_init = hkl_geometry_new_copy(geometry);
  83
  84         getset->detector_init = *detector;
  85
  86         if(getset->sample_init)
  87                 hkl_sample_free(getset->sample_init);
  88         getset->sample_init = hkl_sample_new_copy(sample);
  89
  90         return HKL_SUCCESS;
  91 }
  92
  93 int hkl_pseudo_axis_engine_getter_func_hkl(HklPseudoAxisEngine *self,
  94                                            HklGeometry *geometry,
  95                                            HklDetector const *detector,
  96                                            HklSample const *sample)
  97 {
  98         HklHolder *holder;
  99         HklMatrix RUB;
 100         HklVector hkl, ki, Q;
 101         double min, max;
 102
 103         // update the geometry internals
 104         hkl_geometry_update(geometry);
 105
 106         // R * UB
 107         // for now the 0 holder is the sample holder.
 108         holder = &geometry->holders[0];
 109         hkl_quaternion_to_smatrix(&holder->q, &RUB);
 110         hkl_matrix_times_smatrix(&RUB, &sample->UB);
 111
 112         // kf - ki = Q
 113         hkl_source_compute_ki(&geometry->source, &ki);
 114         hkl_detector_compute_kf(detector, geometry, &Q);
 115         hkl_vector_minus_vector(&Q, &ki);
 116
 117         hkl_matrix_solve(&RUB, &hkl, &Q);
 118
 119         // compute the min and max
 120         min = -1;
 121         max = 1;
 122
 123         // update the pseudoAxes config part
 124         hkl_axis_config_init(&self->pseudoAxes[0].config,
 125                              min, max, hkl.data[0], HKL_FALSE);
 126
 127         hkl_axis_config_init(&self->pseudoAxes[1].config,
 128                              min, max, hkl.data[1], HKL_FALSE);
 129
 130         hkl_axis_config_init(&self->pseudoAxes[2].config,
 131                              min, max, hkl.data[2], HKL_FALSE);
 132
 133         return HKL_SUCCESS;
 134 }
 135
 136 int hkl_pseudo_axis_engine_setter_func_hkl(HklPseudoAxisEngine *self,
 137                                            HklGeometry *geometry,
 138                                            HklDetector *detector,
 139                                            HklSample *sample)
 140 {
 141         hkl_pseudoAxeEngine_prepare_internal(self, geometry, detector,
 142                                              sample);
 143
 144         return hkl_pseudoAxeEngine_solve_function(self, RUBh_minus_Q_func);
 145 }
 146