hkl/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-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 <stdlib.h>
  23 #include <math.h>
  24
  25 #include <gsl/gsl_math.h>
  26 #include <gsl/gsl_sf_trig.h>
  27
  28 #include <hkl/hkl-axis.h>
  29 #include <hkl/hkl-quaternion.h>
  30
  31 /***********/
  32 /* HklAxis */
  33 /***********/
  34 static void hkl_axis_update(HklAxis *self)
  35 {
  36         hkl_quaternion_init_from_angle_and_axe(&self->q,
  37                                                ((HklParameter *)self)->value, &self->axis_v);
  38 }
  39
  40 /*
  41  * given a current position of angle a min and max interval find the closest
  42  * equivalent angle + n delta_angle in a given direction.
  43  * CAUSION angle MUST be in [min, max] otherwise...
  44  */
  45 static void find_angle(double current, double *angle, double *distance,
  46                        double min, double max, double delta_angle)
  47 {
  48         double new_angle = *angle;
  49         double new_distance;
  50
  51         while(new_angle >= min && new_angle <= max) {
  52                 new_distance = fabs(new_angle - current);
  53                 if (new_distance <= *distance) {
  54                         *angle = new_angle;
  55                         *distance = new_distance;
  56                 }
  57                 new_angle += delta_angle;
  58         }
  59 }
  60
  61
  62 /*
  63  * check if the angle or its equivalent is in between [min, max]
  64  */
  65 int hkl_axis_is_value_compatible_with_range(HklAxis const *self)
  66 {
  67         double value;
  68         int res = HKL_FALSE;
  69         HklInterval range;
  70
  71         value = ((HklParameter *)self)->value;
  72         range = ((HklParameter *)self)->range;
  73
  74         if(hkl_interval_length(&range) > 2*M_PI)
  75                 res = HKL_TRUE;
  76         else{
  77                 hkl_interval_angle_restrict_symm(&range);
  78                 value = gsl_sf_angle_restrict_symm(value);
  79
  80                 if(range.min <= range.max){
  81                         if(range.min <= value && range.max >= value)
  82                                 res = HKL_TRUE;
  83                 }else{
  84                         if(value <= range.max || value >= range.min)
  85                                 res = HKL_TRUE;
  86                 }
  87         }
  88         return res;
  89 }
  90
  91 HklAxis *hkl_axis_new(char const *name, HklVector const *axis_v)
  92 {
  93         HklAxis *self = NULL;
  94
  95         self = HKL_MALLOC(HklAxis);
  96
  97         hkl_axis_init(self, name, axis_v);
  98
  99         return self;
 100 }
 101
 102 void hkl_axis_free(HklAxis *self)
 103 {
 104         if(self)
 105                 free(self);
 106 }
 107
 108 void hkl_axis_init(HklAxis *self, char const * name, HklVector const *axis_v)
 109 {
 110         static HklQuaternion q0 = {{1, 0, 0, 0}};
 111
 112         /* base initializer */
 113         hkl_parameter_init((HklParameter *)self, name, -M_PI, 0, M_PI,
 114                            HKL_TRUE, HKL_TRUE,
 115                            &hkl_unit_angle_rad, &hkl_unit_angle_deg);
 116
 117         self->axis_v = *axis_v;
 118         self->q = q0;
 119 }
 120
 121 char const *hkl_axis_get_name(HklAxis const *self)
 122 {
 123         return ((HklParameter *)self)->name;
 124 }
 125
 126 int hkl_axis_get_changed(HklAxis const *self)
 127 {
 128         return ((HklParameter *)self)->changed;
 129 }
 130
 131 void hkl_axis_set_changed(HklAxis *self, int changed)
 132 {
 133         ((HklParameter *)self)->changed = changed;
 134 }
 135
 136 double hkl_axis_get_value(HklAxis const *self)
 137 {
 138         return ((HklParameter *)self)->value;
 139 }
 140
 141 double hkl_axis_get_value_unit(HklAxis const *self)
 142 {
 143         return hkl_parameter_get_value_unit((HklParameter *)self);
 144 }
 145
 146 double hkl_axis_get_value_closest(HklAxis const *self, HklAxis const *axis)
 147 {
 148         double angle = ((HklParameter *)self)->value;
 149
 150         if(hkl_axis_is_value_compatible_with_range(self)){
 151                 if(hkl_interval_length(&((HklParameter *)self)->range) >= 2*M_PI){
 152                         int k;
 153                         double current = ((HklParameter *)axis)->value;
 154                         double distance = fabs(current - angle);
 155                         double delta = 2. * M_PI;
 156                         double min = ((HklParameter *)self)->range.min;
 157                         double max = ((HklParameter *)self)->range.max;
 158
 159                         /* three cases */
 160                         if (angle > max) {
 161                                 k = (int)(floor((max - angle) / delta));
 162                                 angle += k * delta;
 163                                 find_angle(current, &angle, &distance, min, max, -delta);
 164                         } else if (angle < min) {
 165                                 k = (int) (ceil((min - angle) / delta));
 166                                 angle += k * delta;
 167                                 find_angle(current, &angle, &distance, min, max, delta);
 168                         } else {
 169                                 find_angle(current, &angle, &distance, min, max, -delta);
 170                                 find_angle(current, &angle, &distance, min, max, delta);
 171                         }
 172                 }
 173
 174         }else
 175                 angle = GSL_NAN;
 176         return angle;
 177 }
 178
 179 double hkl_axis_get_value_closest_unit(HklAxis const *self, HklAxis const *axis)
 180 {
 181         double factor = hkl_unit_factor(((HklParameter *)self)->unit, ((HklParameter *)self)->punit);
 182         return factor * hkl_axis_get_value_closest(self, axis);
 183 }
 184
 185 double hkl_axis_get_max(HklAxis const *self)
 186 {
 187         return hkl_parameter_get_max((HklParameter *)self);
 188 }
 189
 190 void hkl_axis_get_range_unit(HklAxis const *self, double *min, double *max)
 191 {
 192         hkl_parameter_get_range_unit((HklParameter *)self, min, max);
 193 }
 194
 195 void hkl_axis_set_value(HklAxis *self, double value)
 196 {
 197         hkl_parameter_set_value((HklParameter *)self, value);
 198         hkl_axis_update(self);
 199 }
 200
 201 void hkl_axis_set_value_smallest_in_range(HklAxis *self)
 202 {
 203         double value, min;
 204
 205         value = ((HklParameter *)self)->value;
 206         min = ((HklParameter *)self)->range.min;
 207
 208         if(value < min)
 209                 hkl_axis_set_value(self, value + 2*M_PI*ceil((min - value)/(2*M_PI)));
 210         else
 211                 hkl_axis_set_value(self, value - 2*M_PI*floor((value - min)/(2*M_PI)));
 212 }
 213
 214 void hkl_axis_set_value_unit(HklAxis *self, double value)
 215 {
 216         hkl_parameter_set_value_unit((HklParameter *)self, value);
 217         hkl_axis_update(self);
 218 }
 219
 220 void hkl_axis_set_range(HklAxis *self, double min, double max)
 221 {
 222         hkl_parameter_set_range((HklParameter *)self, min, max);
 223 }
 224
 225 void hkl_axis_set_range_unit(HklAxis *self, double min, double max)
 226 {
 227         hkl_parameter_set_range_unit((HklParameter *)self, min, max);
 228 }
 229
 230 void hkl_axis_randomize(HklAxis *self)
 231 {
 232         hkl_parameter_randomize((HklParameter *)self);
 233         hkl_axis_update(self);
 234 }
 235
 236 void hkl_axis_get_quaternion(HklAxis const *self, HklQuaternion *q)
 237 {
 238         hkl_quaternion_init_from_angle_and_axe(q,
 239                                                ((HklParameter *)self)->value, &self->axis_v);
 240 }
 241
 242 int hkl_axis_is_valid(const HklAxis *self)
 243 {
 244         return hkl_parameter_is_valid((HklParameter *)self);
 245 }
 246
 247 void hkl_axis_fprintf(FILE *f, HklAxis *self)
 248 {
 249         hkl_parameter_fprintf(f, (HklParameter *)self);
 250         hkl_vector_fprintf(f, &self->axis_v);
 251         hkl_quaternion_fprintf(f, &self->q);
 252 }