hkl/hkl-pseudoaxis-common-eulerians.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-2014 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  *          Jens Krüger <Jens.Krueger@frm2.tum.de>
  22  */
  23 #include <gsl/gsl_sf_trig.h>            // for gsl_sf_angle_restrict_symm
  24 #include <math.h>                       // for sin, asin, M_PI_2, tan, etc
  25 #include <stdlib.h>                     // for free
  26 #include <sys/types.h>                  // for uint
  27 #include "hkl-geometry-private.h"
  28 #include "hkl-macros-private.h"         // for HKL_MALLOC
  29 #include "hkl-parameter-private.h"      // for _HklParameter, etc
  30 #include "hkl-pseudoaxis-common-eulerians-private.h"
  31 #include "hkl-pseudoaxis-private.h"     // for _HklPseudoAxis, etc
  32 #include "hkl.h"                        // for HklParameter, HklPseudoAxis, etc
  33 #include "hkl/ccan/array_size/array_size.h"  // for ARRAY_SIZE
  34 #include "hkl/ccan/container_of/container_of.h"  // for container_of
  35 #include "hkl/ccan/darray/darray.h"     // for darray_item
  36
  37 #define HKL_MODE_EULERIANS_ERROR hkl_mode_eulerians_error_quark ()
  38
  39 static GQuark hkl_mode_eulerians_error_quark (void)
  40 {
  41         return g_quark_from_static_string ("hkl-mode-eulerians-error-quark");
  42 }
  43
  44 typedef enum {
  45         HKL_MODE_EULERIANS_ERROR_SET, /* can not set the engine */
  46 } HklModeEuleriansError;
  47
  48 static int kappa_to_eulerian(const double angles[],
  49                              double *omega, double *chi, double *phi,
  50                              double alpha, int solution)
  51 {
  52         const double komega = angles[0];
  53         const double kappa = gsl_sf_angle_restrict_symm(angles[1]);
  54         const double kphi = angles[2];
  55         const double p = atan(tan(kappa/2.) * cos(alpha));
  56
  57         if (solution){
  58                 *omega = komega + p - M_PI_2;
  59                 *chi = 2 * asin(sin(kappa/2.) * sin(alpha));
  60                 *phi = kphi + p + M_PI_2;
  61         }else{
  62                 *omega = komega + p + M_PI_2;
  63                 *chi = -2 * asin(sin(kappa/2.) * sin(alpha));
  64                 *phi = kphi + p - M_PI_2;
  65         }
  66
  67         return TRUE;
  68 }
  69
  70 static int eulerian_to_kappa(const double omega, const double chi, const double phi,
  71                              double angles[],
  72                              double alpha, double solution)
  73 {
  74         int status = TRUE;
  75         double *komega = &angles[0];
  76         double *kappa = &angles[1];
  77         double *kphi = &angles[2];
  78
  79         if (fabs(chi) <= alpha * 2){
  80                 const double p = asin(tan(chi/2.)/tan(alpha));
  81
  82                 if (solution){
  83                         *komega = omega - p + M_PI_2;
  84                         *kappa = 2 * asin(sin(chi/2.)/sin(alpha));
  85                         *kphi = phi - p - M_PI_2;
  86                 }else{
  87                         *komega = omega + p - M_PI_2;
  88                         *kappa = -2 * asin(sin(chi/2.)/sin(alpha));
  89                         *kphi = phi + p + M_PI_2;
  90                 }
  91         }else
  92                 status = FALSE;
  93
  94         return status;
  95 }
  96
  97 void kappa_2_kappap(double komega, double kappa, double kphi, double alpha,
  98                     double *komegap, double *kappap, double *kphip)
  99 {
 100         double p;
 101         double omega;
 102         double phi;
 103
 104         p = atan(tan(kappa/2.) * cos(alpha));
 105         omega = komega + p - M_PI_2;
 106         phi = kphi + p + M_PI_2;
 107
 108         *komegap = gsl_sf_angle_restrict_symm(2*omega - komega);
 109         *kappap = -kappa;
 110         *kphip = gsl_sf_angle_restrict_symm(2*phi - kphi);
 111
 112 }
 113
 114 /***********/
 115 /* HklMode */
 116 /***********/
 117
 118 static int hkl_mode_get_eulerians_real(HklMode *self,
 119                                        HklEngine *engine,
 120                                        HklGeometry *geometry,
 121                                        HklDetector *detector,
 122                                        HklSample *sample,
 123                                        GError **error)
 124 {
 125         HklEngineEulerians *eulerians;
 126         const double angles[] = {
 127                 hkl_parameter_value_get(
 128                         hkl_geometry_get_axis_by_name(geometry, "komega"),
 129                         HKL_UNIT_DEFAULT),
 130                 hkl_parameter_value_get(
 131                         hkl_geometry_get_axis_by_name(geometry, "kappa"),
 132                         HKL_UNIT_DEFAULT),
 133                 hkl_parameter_value_get(
 134                         hkl_geometry_get_axis_by_name(geometry, "kphi"),
 135                         HKL_UNIT_DEFAULT),
 136         };
 137         double values[3];
 138         double solution;
 139         HklParameter *parameter;
 140
 141         hkl_geometry_update(geometry);
 142
 143         solution = darray_item(self->parameters, 0)->_value;
 144
 145         eulerians = container_of(engine, HklEngineEulerians, engine);
 146         kappa_to_eulerian(angles,
 147                           &eulerians->omega->_value,
 148                           &eulerians->chi->_value,
 149                           &eulerians->phi->_value,
 150                           50 * HKL_DEGTORAD, solution);
 151
 152         return TRUE;
 153 }
 154
 155 static int hkl_mode_set_eulerians_real(HklMode *self,
 156                                        HklEngine *engine,
 157                                        HklGeometry *geometry,
 158                                        HklDetector *detector,
 159                                        HklSample *sample,
 160                                        GError **error)
 161 {
 162         double solution;
 163         HklEngineEulerians *engine_eulerians;
 164         double angles[3];
 165
 166         solution = darray_item(self->parameters, 0)->_value;
 167         engine_eulerians = container_of(engine, HklEngineEulerians, engine);
 168         if(!eulerian_to_kappa(engine_eulerians->omega->_value,
 169                               engine_eulerians->chi->_value,
 170                               engine_eulerians->phi->_value,
 171                               angles, 50 * HKL_DEGTORAD, solution)){
 172                 g_set_error(error,
 173                             HKL_MODE_EULERIANS_ERROR,
 174                             HKL_MODE_EULERIANS_ERROR_SET,
 175                             "unreachable solution : 0° < chi < 50°");
 176                 return FALSE;
 177         }else
 178                 hkl_engine_add_geometry(engine, angles);
 179
 180         return TRUE;
 181 }
 182
 183
 184 static HklMode *mode_eulerians()
 185 {
 186         HklMode *mode;
 187         static const char *axes[] = {"komega", "kappa", "kphi"};
 188         static const HklParameter parameters[] = {
 189                 {HKL_PARAMETER_DEFAULTS, .name="solutions", .range = {.max = 1}, ._value = 1,},
 190         };
 191         static const HklModeInfo info = {
 192                 HKL_MODE_INFO_WITH_PARAMS("eulerians", axes, axes, parameters),
 193         };
 194         static const HklModeOperations operations = {
 195                 HKL_MODE_OPERATIONS_DEFAULTS,
 196                 .get = hkl_mode_get_eulerians_real,
 197                 .set = hkl_mode_set_eulerians_real,
 198         };
 199
 200         return hkl_mode_new(&info, &operations, TRUE);
 201 };
 202
 203 /*************/
 204 /* HklEngine */
 205 /*************/
 206
 207 static void hkl_engine_eulerians_free_real(HklEngine *base)
 208 {
 209         HklEngineEulerians *self;
 210
 211         self = container_of(base, HklEngineEulerians, engine);
 212         hkl_engine_release(&self->engine);
 213         free(self);
 214 }
 215
 216 HklEngine *hkl_engine_eulerians_new(void)
 217 {
 218         HklEngineEulerians *self;
 219         HklMode *mode;
 220         static const HklPseudoAxis omega = {
 221                 .parameter = { HKL_PARAMETER_DEFAULTS_ANGLE, .name = "omega"}
 222         };
 223         static const HklPseudoAxis chi = {
 224                 .parameter = { HKL_PARAMETER_DEFAULTS_ANGLE, .name = "chi"}
 225         };
 226         static const HklPseudoAxis phi = {
 227                 .parameter = {HKL_PARAMETER_DEFAULTS_ANGLE, .name = "phi"}
 228         };
 229         static const HklPseudoAxis *pseudo_axes[] = {&omega, &chi, &phi};
 230         static HklEngineInfo info = {
 231                 .name = "eulerians",
 232                 .pseudo_axes = DARRAY(pseudo_axes),
 233         };
 234         static HklEngineOperations operations = {
 235                 HKL_ENGINE_OPERATIONS_DEFAULTS,
 236                 .free=hkl_engine_eulerians_free_real,
 237         };
 238
 239         self = HKL_MALLOC(HklEngineEulerians);
 240         hkl_engine_init(&self->engine, &info, &operations);
 241
 242         /* add the pseudo axes with the new API */
 243         self->omega = register_pseudo_axis(&self->engine, &omega.parameter);
 244         self->chi = register_pseudo_axis(&self->engine, &chi.parameter);
 245         self->phi = register_pseudo_axis(&self->engine, &phi.parameter);
 246
 247         /* eulerians [default] */
 248         mode = mode_eulerians();
 249         hkl_engine_add_mode(&self->engine, mode);
 250         hkl_engine_mode_set(&self->engine, mode);
 251
 252         return &self->engine;
 253 }