src/hkl-axis.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 <stdlib.h>
  23 #include <math.h>
  24
  25 #include <gsl/gsl_math.h>
  26
  27 #include <hkl/hkl-axis.h>
  28 #include <hkl/hkl-quaternion.h>
  29
  30 /***********/
  31 /* HklAxis */
  32 /***********/
  33 static void hkl_axis_update(HklAxis *self)
  34 {
  35         hkl_quaternion_from_angle_and_axe(&self->q, self->parent.value, &self->axis_v);
  36 }
  37
  38 /*
  39  * given a current position of angle a min and max interval find the closest
  40  * equivalent angle + n delta_angle in a given direction.
  41  * CAUSION angle MUST be in [min, max] otherwise...
  42  */
  43 static void find_angle(double current, double *angle, double *distance,
  44                        double min, double max, double delta_angle)
  45 {
  46         double new_angle = *angle;
  47         double new_distance = *distance;
  48
  49         while(new_angle >= min && new_angle <= max) {
  50                 new_distance = fabs(new_angle - current);
  51                 if (new_distance <= *distance) {
  52                         *angle = new_angle;
  53                         *distance = new_distance;
  54                 }
  55                 new_angle += delta_angle;
  56         }
  57 }
  58
  59
  60 /*
  61  * check if the angle or its equivalent is in between [min, max]
  62  */
  63 static int hkl_axis_is_value_compatible_with_range(HklAxis const *self)
  64 {
  65         double c = cos(self->parent.value);
  66         HklInterval c_r = self->parent.range;
  67
  68         hkl_interval_cos(&c_r);
  69
  70         if (c_r.min <= c && c <= c_r.max) {
  71                 double s = sin(self->parent.value);
  72                 HklInterval s_r = self->parent.range;
  73
  74                 hkl_interval_sin(&s_r);
  75
  76                 if (s_r.min <= s && s <= s_r.max)
  77                         return HKL_TRUE;
  78         }
  79         return HKL_FALSE;
  80 }
  81
  82 HklAxis *hkl_axis_new(char const *name, HklVector const *axis_v)
  83 {
  84         HklAxis *self = NULL;
  85
  86         self = calloc(1, sizeof(*self));
  87         if (!self)
  88                 die("Can not allocate memory for an HklAxis");
  89
  90         hkl_axis_init(self, name, axis_v);
  91
  92         return self;
  93 }
  94
  95 void hkl_axis_free(HklAxis *self)
  96 {
  97         if(self)
  98                 free(self);
  99 }
 100
 101 void hkl_axis_init(HklAxis *self, char const * name, HklVector const *axis_v)
 102 {
 103         static HklQuaternion q0 = {{1, 0, 0, 0}};
 104
 105         // base initializer
 106         hkl_parameter_init((HklParameter *)self, name, -M_PI, 0, M_PI,
 107                            HKL_FALSE, HKL_TRUE,
 108                            &hkl_unit_angle_rad, &hkl_unit_angle_deg);
 109
 110         self->axis_v = *axis_v;
 111         self->q = q0;
 112 }
 113
 114 double hkl_axis_get_value_unit(HklAxis const *self)
 115 {
 116         return hkl_parameter_get_value_unit(&self->parent);
 117 }
 118
 119 double hkl_axis_get_value_closest(HklAxis const *self, HklAxis const *axis)
 120 {
 121         double angle = self->parent.value;
 122
 123         if(hkl_axis_is_value_compatible_with_range(self)){
 124                 if(hkl_interval_length(&self->parent.range) >= 2*M_PI){
 125                         int k;
 126                         double current = axis->parent.value;
 127                         double distance = fabs(current - angle);
 128                         double delta = 2. * M_PI;
 129                         double min = self->parent.range.min;
 130                         double max = self->parent.range.max;
 131
 132                         // three cases
 133                         if (angle > max) {
 134                                 k = (int)(floor((max - angle) / delta));
 135                                 angle += k * delta;
 136                                 find_angle(current, &angle, &distance, min, max, -delta);
 137                         } else if (angle < min) {
 138                                 k = (int) (ceil((min - angle) / delta));
 139                                 angle += k * delta;
 140                                 find_angle(current, &angle, &distance, min, max, delta);
 141                         } else {
 142                                 find_angle(current, &angle, &distance, min, max, -delta);
 143                                 find_angle(current, &angle, &distance, min, max, delta);
 144                         }
 145                 }
 146
 147         }else
 148                 angle = GSL_NAN;
 149         return angle;
 150 }
 151
 152 double hkl_axis_get_value_closest_unit(HklAxis const *self, HklAxis const *axis)
 153 {
 154         double factor = hkl_unit_factor(self->parent.unit, self->parent.punit);
 155         return factor * hkl_axis_get_value_closest(self, axis);
 156 }
 157
 158 void hkl_axis_get_range_unit(HklAxis const *self, double *min, double *max)
 159 {
 160         hkl_parameter_get_range_unit(&self->parent, min, max);
 161 }
 162
 163 void hkl_axis_set_value(HklAxis *self, double value)
 164 {
 165         hkl_parameter_set_value(&self->parent, value);
 166         hkl_axis_update(self);
 167 }
 168
 169 void hkl_axis_set_value_unit(HklAxis *self, double value)
 170 {
 171         hkl_parameter_set_value_unit(&self->parent, value);
 172         hkl_axis_update(self);
 173 }
 174
 175 void hkl_axis_set_range(HklAxis *self, double min, double max)
 176 {
 177         hkl_parameter_set_range(&self->parent, min, max);
 178 }
 179
 180 void hkl_axis_set_range_unit(HklAxis *self, double min, double max)
 181 {
 182         hkl_parameter_set_range_unit(&self->parent, min, max);
 183 }
 184
 185 void hkl_axis_randomize(HklAxis *self)
 186 {
 187         hkl_parameter_randomize(&self->parent);
 188         hkl_axis_update(self);
 189 }
 190
 191 void hkl_axis_get_quaternion(HklAxis const *self, HklQuaternion *q)
 192 {
 193         hkl_quaternion_from_angle_and_axe(q, ((HklParameter *)self)->value,
 194                                           &self->axis_v);
 195 }
 196
 197 void hkl_axis_fprintf(FILE *f, HklAxis *self)
 198 {
 199         hkl_parameter_fprintf(f, &self->parent);
 200         hkl_vector_fprintf(f, &self->axis_v);
 201         hkl_quaternion_fprintf(f, &self->q);
 202 }