1 /* This file is part of the hkl library.
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.
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.
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/>.
16 * Copyright (C) 2003-2019, 2022 Synchrotron SOLEIL
17 * L'Orme des Merisiers Saint-Aubin
18 * BP 48 91192 GIF-sur-YVETTE CEDEX
20 * Authors: Picca Frédéric-Emmanuel <picca@synchrotron-soleil.fr>
22 #include <gsl/gsl_errno.h> // for ::GSL_SUCCESS
23 #include <gsl/gsl_vector_double.h> // for gsl_vector
24 #include <gsl/gsl_sys.h> // for gsl_isnan
25 #include <stdio.h> // for fprintf, stderr
26 #include <stdlib.h> // for NULL, exit, free
27 #include <sys/types.h> // for uint
28 #include "hkl-detector-private.h" // for hkl_detector_compute_kf
29 #include "hkl-geometry-private.h" // for HklHolder, _HklGeometry, etc
30 #include "hkl-macros-private.h" // for HKL_MALLOC, hkl_assert, etc
31 #include "hkl-matrix-private.h" // for hkl_matrix_solve, etc
32 #include "hkl-parameter-private.h"
33 #include "hkl-pseudoaxis-auto-private.h" // for HklModeAutoInfo, etc
34 #include "hkl-pseudoaxis-common-psi-private.h" // for HklEnginePsi, etc
35 #include "hkl-pseudoaxis-private.h" // for _HklEngine, HklEngineInfo, etc
36 #include "hkl-quaternion-private.h" // for hkl_quaternion_to_matrix
37 #include "hkl-sample-private.h" // for _HklSample
38 #include "hkl-source-private.h" // for hkl_source_compute_ki
39 #include "hkl-vector-private.h" // for HklVector, etc
40 #include "hkl.h" // for HklEngine, HklGeometry, etc
41 #include "hkl/ccan/array_size/array_size.h" // for ARRAY_SIZE
42 #include "hkl/ccan/container_of/container_of.h" // for container_of
43 #include "hkl/ccan/darray/darray.h" // for darray_item
45 #define HKL_MODE_PSI_ERROR hkl_mode_psi_error_quark ()
47 static GQuark
hkl_mode_psi_error_quark (void)
49 return g_quark_from_static_string ("hkl-mode-psi-error-quark");
53 HKL_MODE_PSI_ERROR_INIT
, /* can not init the engine */
54 HKL_MODE_PSI_ERROR_GET
, /* can not get the engine */
57 /***********************/
58 /* numerical functions */
59 /***********************/
61 int _psi_func(const gsl_vector
*x
, void *params
, gsl_vector
*f
)
64 HklVector dhkl0
, hkl1
;
65 HklVector ki
, kf
, Q
, n
;
67 HklEngine
*engine
= params
;
68 HklEnginePsi
*psi_engine
= container_of(engine
, HklEnginePsi
, engine
);
69 HklModePsi
*modepsi
= container_of(engine
->mode
, HklModePsi
, parent
);
70 HklHolder
*sample_holder
;
72 CHECK_NAN(x
->data
, x
->size
);
74 /* update the workspace from x; */
75 set_geometry_axes(engine
, x
->data
);
78 ki
= hkl_geometry_ki_get(engine
->geometry
);
79 kf
= hkl_geometry_kf_get(engine
->geometry
, engine
->detector
);
81 hkl_vector_minus_vector(&Q
, &ki
);
82 if (hkl_vector_is_null(&Q
)){
89 /* for now the 0 holder is the sample holder. */
90 sample_holder
= hkl_geometry_sample_holder_get(engine
->geometry
,
92 hkl_quaternion_to_matrix(&sample_holder
->q
, &RUB
);
93 hkl_matrix_times_matrix(&RUB
, &engine
->sample
->UB
);
96 hkl_matrix_solve(&RUB
, &dhkl0
, &Q
);
97 hkl_vector_minus_vector(&dhkl0
, &modepsi
->hkl0
);
99 /* compute the intersection of the plan P(kf, ki) and PQ (normal Q) */
101 * now that dhkl0 have been computed we can use a
102 * normalized Q to compute n and psi
104 hkl_vector_normalize(&Q
);
106 hkl_vector_vectorial_product(&n
, &ki
);
107 hkl_vector_vectorial_product(&n
, &Q
);
109 /* compute hkl1 in the laboratory referentiel */
110 /* for now the 0 holder is the sample holder. */
111 for(unsigned int i
=0; i
<3; ++i
)
112 hkl1
.data
[i
] = darray_item(engine
->mode
->parameters
, i
)->_value
;
114 hkl_matrix_times_vector(&engine
->sample
->UB
, &hkl1
);
115 hkl_vector_rotated_quaternion(&hkl1
, &sample_holder
->q
);
117 /* project hkl1 on the plan of normal Q */
118 hkl_vector_project_on_plan(&hkl1
, &Q
);
119 if (hkl_vector_is_null(&hkl1
)){
120 /* hkl1 colinear with Q */
121 f
->data
[0] = dhkl0
.data
[0];
122 f
->data
[1] = dhkl0
.data
[1];
123 f
->data
[2] = dhkl0
.data
[2];
126 f
->data
[0] = dhkl0
.data
[0];
127 f
->data
[1] = dhkl0
.data
[1];
128 f
->data
[2] = dhkl0
.data
[2];
129 f
->data
[3] = psi_engine
->psi
->_value
- hkl_vector_oriented_angle(&n
, &hkl1
, &Q
);
135 static int hkl_mode_initialized_set_psi_real(HklMode
*self
,
137 HklGeometry
*geometry
,
138 HklDetector
*detector
,
145 HklModePsi
*psi_mode
= container_of(self
, HklModePsi
, parent
);
146 HklHolder
*sample_holder
= hkl_geometry_sample_holder_get(geometry
, sample
);
148 hkl_error (error
== NULL
|| *error
== NULL
);
151 /* update the geometry internals */
152 hkl_geometry_update(geometry
);
155 /* for now the 0 holder is the sample holder. */
156 hkl_quaternion_to_matrix(&sample_holder
->q
, &RUB
);
157 hkl_matrix_times_matrix(&RUB
, &sample
->UB
);
160 ki
= hkl_geometry_ki_get(geometry
);
161 psi_mode
->Q0
= hkl_geometry_kf_get(geometry
, detector
);
162 hkl_vector_minus_vector(&psi_mode
->Q0
, &ki
);
163 if (hkl_vector_is_null(&psi_mode
->Q0
)){
166 HKL_MODE_PSI_ERROR_INIT
,
167 "can not initialize the \"%s\" engine when hkl is null",
172 hkl_matrix_solve(&RUB
, &psi_mode
->hkl0
, &psi_mode
->Q0
);
175 self
->initialized
= initialized
;
180 static int hkl_mode_get_psi_real(HklMode
*base
,
182 HklGeometry
*geometry
,
183 HklDetector
*detector
,
193 if (!base
|| !engine
|| !engine
->mode
|| !geometry
|| !detector
|| !sample
){
196 HKL_MODE_PSI_ERROR_GET
,
201 /* get kf, ki and Q */
202 ki
= hkl_geometry_ki_get(geometry
);
203 kf
= hkl_geometry_kf_get(geometry
, detector
);
205 hkl_vector_minus_vector(&Q
, &ki
);
206 if (hkl_vector_is_null(&Q
)){
209 HKL_MODE_PSI_ERROR_GET
,
210 "can not compute psi when hkl is null (kf == ki)");
213 HklHolder
*sample_holder
= hkl_geometry_sample_holder_get(geometry
, sample
);
215 /* needed for a problem of precision */
216 hkl_vector_normalize(&Q
);
218 /* compute the intersection of the plan P(kf, ki) and PQ (normal Q) */
220 hkl_vector_vectorial_product(&n
, &ki
);
221 hkl_vector_vectorial_product(&n
, &Q
);
223 /* compute hkl1 in the laboratory referentiel */
224 /* the geometry was already updated in the detector compute kf */
225 /* for now the 0 holder is the sample holder. */
226 for(unsigned int i
=0; i
<3; ++i
)
227 hkl1
.data
[i
] = darray_item(base
->parameters
, i
)->_value
;
229 hkl_matrix_times_vector(&sample
->UB
, &hkl1
);
230 hkl_vector_rotated_quaternion(&hkl1
, &sample_holder
->q
);
232 /* project hkl1 on the plan of normal Q */
233 hkl_vector_project_on_plan(&hkl1
, &Q
);
235 if (hkl_vector_is_null(&hkl1
)){
238 HKL_MODE_PSI_ERROR_GET
,
239 "can not compute psi when Q and the ref vector are colinear");
242 HklEnginePsi
*psi_engine
= container_of(engine
, HklEnginePsi
, engine
);
244 /* compute the angle beetween hkl1 and n */
245 psi_engine
->psi
->_value
= hkl_vector_oriented_angle(&n
, &hkl1
, &Q
);
252 HklMode
*hkl_mode_psi_new(const HklModeAutoInfo
*auto_info
)
254 static const HklModeOperations operations
= {
255 HKL_MODE_OPERATIONS_AUTO_DEFAULTS
,
256 .capabilities
= HKL_ENGINE_CAPABILITIES_READABLE
| HKL_ENGINE_CAPABILITIES_WRITABLE
| HKL_ENGINE_CAPABILITIES_INITIALIZABLE
,
257 .initialized_set
= hkl_mode_initialized_set_psi_real
,
258 .get
= hkl_mode_get_psi_real
,
262 if (darray_size(auto_info
->info
.axes_w
) != 4){
263 fprintf(stderr
, "This generic HklModePsi need exactly 4 axes");
267 self
= g_new(HklModePsi
, 1);
269 /* the base constructor; */
270 hkl_mode_auto_init(&self
->parent
,
274 return &self
->parent
;
281 static void hkl_engine_psi_free_real(HklEngine
*base
)
283 HklEnginePsi
*self
=container_of(base
, HklEnginePsi
, engine
);
284 hkl_engine_release(&self
->engine
);
288 HklEngine
*hkl_engine_psi_new(HklEngineList
*engines
)
290 HklEnginePsi
*self
= g_new(HklEnginePsi
, 1);
291 static const HklParameter psi
= {
292 HKL_PARAMETER_DEFAULTS_ANGLE
, .name
= "psi",
293 .description
= "angle between the reference vector and the diffraction plan",
295 static const HklParameter
*pseudo_axes
[] = {&psi
};
296 static const HklEngineInfo info
= {
297 HKL_ENGINE_INFO("psi",
299 HKL_ENGINE_DEPENDENCIES_AXES
| HKL_ENGINE_DEPENDENCIES_ENERGY
| HKL_ENGINE_DEPENDENCIES_SAMPLE
),
301 static const HklEngineOperations operations
= {
302 HKL_ENGINE_OPERATIONS_DEFAULTS
,
303 .free
=hkl_engine_psi_free_real
,
306 hkl_engine_init(&self
->engine
, &info
, &operations
, engines
);
308 self
->psi
= register_pseudo_axis(&self
->engine
, engines
, &psi
);
310 return &self
->engine
;