1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
|
/*
* basis.c
*
* Handle basis structures
*
* (c) 2007 Thomas White <taw27@cam.ac.uk>
*
* dtr - Diffraction Tomography Reconstruction
*
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include "reflections.h"
#include "basis.h"
double basis_efom(ReflectionList *reflectionlist, Basis *basis) {
int n_indexed, n_counted;
Reflection *cur;
cur = reflectionlist->reflections;
n_indexed = 0;
n_counted = 0;
while ( cur ) {
if ( cur->type == REFLECTION_NORMAL ) {
/* Can this basis "approximately" account for this reflection? */
double det;
double a11, a12, a13, a21, a22, a23, a31, a32, a33;
double h, k, l;
/* Set up the coordinate transform from hkl to xyz */
a11 = basis->a.x; a12 = basis->a.y; a13 = basis->a.z;
a21 = basis->b.x; a22 = basis->b.y; a23 = basis->b.z;
a31 = basis->c.x; a32 = basis->c.y; a33 = basis->c.z;
/* Invert the matrix to get hkl from xyz */
det = a11*(a22*a33 - a23*a32) - a12*(a21*a33 - a23*a31) + a13*(a21*a32 - a22*a31);
h = ((a22*a33-a23*a32)*cur->x + (a23*a31-a21*a33)*cur->y + (a21*a32-a22*a31)*cur->z) / det;
k = ((a13*a32-a12*a33)*cur->x + (a11*a33-a13*a31)*cur->y + (a12*a31-a11*a32)*cur->z) / det;
l = ((a12*a23-a13*a22)*cur->x + (a13*a21-a11*a23)*cur->y + (a11*a22-a12*a21)*cur->z) / det;
/* Calculate the deviations in terms of |a|, |b| and |c| */
h = fabs(h); k = fabs(k); l = fabs(l);
h -= floor(h); k -= floor(k); l -= floor(l);
if ( h == 1.0 ) h = 0.0;
if ( k == 1.0 ) k = 0.0;
if ( l == 1.0 ) l = 0.0;
/* Define "approximately" here. Circle in basis space becomes an ellipsoid in reciprocal space */
if ( h*h + k*k + l*l <= 0.1*0.1*0.1 ) n_indexed++;
n_counted++;
}
cur = cur->next;
}
return (double)n_indexed / n_counted;
}
Basis basis_add(Basis u, Basis v) {
Basis ans;
ans.a.x = u.a.x + v.a.x; ans.a.y = u.a.y + v.a.y; ans.a.z = u.a.z + v.a.z;
ans.b.x = u.b.x + v.b.x; ans.b.y = u.b.y + v.b.y; ans.b.z = u.b.z + v.b.z;
ans.c.x = u.c.x + v.c.x; ans.c.y = u.c.y + v.c.y; ans.c.z = u.c.z + v.c.z;
return ans;
}
UnitCell basis_get_cell(Basis *basis) {
UnitCell cell;
gsl_matrix *m;
gsl_matrix *inv;
gsl_permutation *perm;
double ax, ay, az, bx, by, bz, cx, cy, cz;
int s;
m = gsl_matrix_alloc(3, 3);
gsl_matrix_set(m, 0, 0, basis->a.x); gsl_matrix_set(m, 0, 1, basis->b.x); gsl_matrix_set(m, 0, 2, basis->c.x);
gsl_matrix_set(m, 1, 0, basis->a.y); gsl_matrix_set(m, 1, 1, basis->b.y); gsl_matrix_set(m, 1, 2, basis->c.y);
gsl_matrix_set(m, 2, 0, basis->a.z); gsl_matrix_set(m, 2, 1, basis->b.z); gsl_matrix_set(m, 2, 2, basis->c.z);
gsl_matrix_transpose(m);
perm = gsl_permutation_alloc(m->size1);
inv = gsl_matrix_alloc(m->size1, m->size2);
gsl_linalg_LU_decomp(m, perm, &s);
gsl_linalg_LU_invert(m, perm, inv);
gsl_permutation_free(perm);
gsl_matrix_free(m);
ax = gsl_matrix_get(inv, 0, 0); bx = gsl_matrix_get(inv, 0, 1); cx = gsl_matrix_get(inv, 0, 2);
ay = gsl_matrix_get(inv, 1, 0); by = gsl_matrix_get(inv, 1, 1); cy = gsl_matrix_get(inv, 1, 2);
az = gsl_matrix_get(inv, 2, 0); bz = gsl_matrix_get(inv, 2, 1); cz = gsl_matrix_get(inv, 2, 2);
cell.a = sqrt(ax*ax + ay*ay + az*az);
cell.b = sqrt(bx*bx + by*by + bz*bz);
cell.c = sqrt(cx*cx + cy*cy + cz*cz);
cell.alpha = acos((bx*cx + by*cy + bz*cz)/(cell.b * cell.c));
cell.beta = acos((ax*cx + ay*cy + az*cz)/(cell.a * cell.c));
cell.gamma = acos((bx*ax + by*ay + bz*az)/(cell.b * cell.a));
gsl_matrix_free(inv);
return cell;
}
|
