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-2017 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 typedef struct _HklModeEulerians HklModeEulerians;
  38 typedef struct _HklEngineEulerians HklEngineEulerians;
  39
  40 struct _HklModeEulerians
  41 {
  42         HklMode parent;
  43         HklParameter *solutions;
  44 };
  45
  46 struct _HklEngineEulerians
  47 {
  48         HklEngine engine;
  49         HklParameter *omega;
  50         HklParameter *chi;
  51         HklParameter *phi;
  52 };
  53
  54 #define HKL_MODE_EULERIANS_ERROR hkl_mode_eulerians_error_quark ()
  55
  56 static GQuark hkl_mode_eulerians_error_quark (void)
  57 {
  58         return g_quark_from_static_string ("hkl-mode-eulerians-error-quark");
  59 }
  60
  61 typedef enum {
  62         HKL_MODE_EULERIANS_ERROR_SET, /* can not set the engine */
  63 } HklModeEuleriansError;
  64
  65 static int kappa_to_eulerian(const double angles[],
  66                              double *omega, double *chi, double *phi,
  67                              double alpha, int solution)
  68 {
  69         const double komega = angles[0];
  70         const double kappa = gsl_sf_angle_restrict_symm(angles[1]);
  71         const double kphi = angles[2];
  72         const double p = atan(tan(kappa/2.) * cos(alpha));
  73
  74         if (solution){
  75                 *omega = komega + p - M_PI_2;
  76                 *chi = 2 * asin(sin(kappa/2.) * sin(alpha));
  77                 *phi = kphi + p + M_PI_2;
  78         }else{
  79                 *omega = komega + p + M_PI_2;
  80                 *chi = -2 * asin(sin(kappa/2.) * sin(alpha));
  81                 *phi = kphi + p - M_PI_2;
  82         }
  83
  84         return TRUE;
  85 }
  86
  87 static int eulerian_to_kappa(const double omega, const double chi, const double phi,
  88                              double angles[],
  89                              double alpha, double solution)
  90 {
  91         int status = TRUE;
  92         double *komega = &angles[0];
  93         double *kappa = &angles[1];
  94         double *kphi = &angles[2];
  95
  96         if (fabs(chi) <= alpha * 2){
  97                 const double p = asin(tan(chi/2.)/tan(alpha));
  98
  99                 if (solution){
 100                         *komega = omega - p + M_PI_2;
 101                         *kappa = 2 * asin(sin(chi/2.)/sin(alpha));
 102                         *kphi = phi - p - M_PI_2;
 103                 }else{
 104                         *komega = omega + p - M_PI_2;
 105                         *kappa = -2 * asin(sin(chi/2.)/sin(alpha));
 106                         *kphi = phi + p + M_PI_2;
 107                 }
 108         }else
 109                 status = FALSE;
 110
 111         return status;
 112 }
 113
 114 void kappa_2_kappap(double komega, double kappa, double kphi, double alpha,
 115                     double *komegap, double *kappap, double *kphip)
 116 {
 117         double p = atan(tan(kappa/2.) * cos(alpha));
 118         double omega = komega + p - M_PI_2;
 119         double phi = kphi + p + M_PI_2;
 120
 121         *komegap = gsl_sf_angle_restrict_symm(2*omega - komega);
 122         *kappap = -kappa;
 123         *kphip = gsl_sf_angle_restrict_symm(2*phi - kphi);
 124
 125 }
 126
 127 /***********/
 128 /* HklMode */
 129 /***********/
 130
 131 static int hkl_mode_get_eulerians_real(HklMode *base,
 132                                        HklEngine *engine,
 133                                        HklGeometry *geometry,
 134                                        HklDetector *detector,
 135                                        HklSample *sample,
 136                                        GError **error)
 137 {
 138         HklModeEulerians *self = container_of(base, HklModeEulerians, parent);
 139         HklEngineEulerians *eulerians = container_of(engine, HklEngineEulerians, engine);;
 140         const double angles[] = {
 141                 hkl_parameter_value_get(
 142                         hkl_geometry_get_axis_by_name(geometry, "komega"),
 143                         HKL_UNIT_DEFAULT),
 144                 hkl_parameter_value_get(
 145                         hkl_geometry_get_axis_by_name(geometry, "kappa"),
 146                         HKL_UNIT_DEFAULT),
 147                 hkl_parameter_value_get(
 148                         hkl_geometry_get_axis_by_name(geometry, "kphi"),
 149                         HKL_UNIT_DEFAULT),
 150         };
 151
 152         hkl_geometry_update(geometry);
 153
 154         kappa_to_eulerian(angles,
 155                           &eulerians->omega->_value,
 156                           &eulerians->chi->_value,
 157                           &eulerians->phi->_value,
 158                           50 * HKL_DEGTORAD, self->solutions->_value);
 159
 160         return TRUE;
 161 }
 162
 163 static int hkl_mode_set_eulerians_real(HklMode *base,
 164                                        HklEngine *engine,
 165                                        HklGeometry *geometry,
 166                                        HklDetector *detector,
 167                                        HklSample *sample,
 168                                        GError **error)
 169 {
 170         HklModeEulerians *self = container_of(base, HklModeEulerians, parent);
 171         HklEngineEulerians *eulerians = container_of(engine, HklEngineEulerians, engine);;
 172         double angles[3];
 173
 174         if(!eulerian_to_kappa(eulerians->omega->_value,
 175                               eulerians->chi->_value,
 176                               eulerians->phi->_value,
 177                               angles, 50 * HKL_DEGTORAD,
 178                               self->solutions->_value)){
 179                 g_set_error(error,
 180                             HKL_MODE_EULERIANS_ERROR,
 181                             HKL_MODE_EULERIANS_ERROR_SET,
 182                             "unreachable solution : 0° < chi < 50°");
 183                 return FALSE;
 184         }else
 185                 hkl_engine_add_geometry(engine, angles);
 186
 187         return TRUE;
 188 }
 189
 190 static HklMode *mode_eulerians()
 191 {
 192         static const char *axes[] = {"komega", "kappa", "kphi"};
 193         static const HklParameter parameters[] = {
 194                 { HKL_PARAMETER_DEFAULTS,.name = "solutions", ._value = 1,
 195                   .description = "(0/1) to select the first or second solution",
 196                   .range = { .max = 1 },
 197                 },
 198         };
 199         static const HklModeInfo info = {
 200                 HKL_MODE_INFO_WITH_PARAMS("eulerians", axes, axes, parameters),
 201         };
 202         static const HklModeOperations operations = {
 203                 HKL_MODE_OPERATIONS_DEFAULTS,
 204                 .get = hkl_mode_get_eulerians_real,
 205                 .set = hkl_mode_set_eulerians_real,
 206         };
 207
 208         HklModeEulerians *self = HKL_MALLOC(HklModeEulerians);
 209
 210         /* the base constructor; */
 211         hkl_mode_init(&self->parent,
 212                       &info,
 213                       &operations, TRUE);
 214
 215         self->solutions = register_mode_parameter(&self->parent, 0);
 216
 217         return &self->parent;
 218 };
 219
 220 /*************/
 221 /* HklEngine */
 222 /*************/
 223
 224 static void hkl_engine_eulerians_free_real(HklEngine *base)
 225 {
 226         HklEngineEulerians *self;
 227
 228         self = container_of(base, HklEngineEulerians, engine);
 229         hkl_engine_release(&self->engine);
 230         free(self);
 231 }
 232
 233 HklEngine *hkl_engine_eulerians_new(HklEngineList *engines)
 234 {
 235         HklEngineEulerians *self;
 236         HklMode *mode;
 237         static const HklParameter omega = {
 238                 HKL_PARAMETER_DEFAULTS_ANGLE, .name = "omega",
 239                 .description = "omega equivalent for a four circle eulerian geometry",
 240         };
 241         static const HklParameter chi = {
 242                 HKL_PARAMETER_DEFAULTS_ANGLE, .name = "chi",
 243                 .description = "chi equivalent for a four circle eulerian geometry",
 244         };
 245         static const HklParameter phi = {
 246                 HKL_PARAMETER_DEFAULTS_ANGLE, .name = "phi",
 247                 .description = "phi equivalent for a four circle eulerian geometry",
 248         };
 249         static const HklParameter *pseudo_axes[] = {&omega, &chi, &phi};
 250         static HklEngineInfo info = {
 251                 HKL_ENGINE_INFO("eulerians",
 252                                 pseudo_axes,
 253                                 HKL_ENGINE_DEPENDENCIES_AXES),
 254         };
 255         static HklEngineOperations operations = {
 256                 HKL_ENGINE_OPERATIONS_DEFAULTS,
 257                 .free=hkl_engine_eulerians_free_real,
 258         };
 259
 260         self = HKL_MALLOC(HklEngineEulerians);
 261         hkl_engine_init(&self->engine, &info, &operations, engines);
 262
 263         /* add the pseudo axes with the new API */
 264         self->omega = register_pseudo_axis(&self->engine, engines, &omega);
 265         self->chi = register_pseudo_axis(&self->engine, engines, &chi);
 266         self->phi = register_pseudo_axis(&self->engine, engines, &phi);
 267
 268         /* eulerians [default] */
 269         mode = mode_eulerians();
 270         hkl_engine_add_mode(&self->engine, mode);
 271         hkl_engine_mode_set(&self->engine, mode);
 272
 273         return &self->engine;
 274 }