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
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
|
/*
* 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"
#include "utils.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;
}
static void basis_print(Basis *cell) {
printf("%12.8f %12.8f %12.8f\n", cell->a.x/1e9, cell->a.y/1e9, cell->a.z/1e9);
printf("%12.8f %12.8f %12.8f\n", cell->b.x/1e9, cell->b.y/1e9, cell->b.z/1e9);
printf("%12.8f %12.8f %12.8f\n", cell->c.x/1e9, cell->c.y/1e9, cell->c.z/1e9);
}
static void cell_print(UnitCell *cell) {
printf("%12.8f %12.8f %12.8f nm\n", cell->a*1e9, cell->b*1e9, cell->c*1e9);
printf("%12.8f %12.8f %12.8f deg\n", rad2deg(cell->alpha), rad2deg(cell->beta), rad2deg(cell->gamma));
}
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;
// basis->a.x = 0.5; basis->a.y = 0.0; basis->a.z = 0.0;
// basis->b.x = 0.0; basis->b.y = 0.2; basis->b.z = 0.0;
// basis->c.x = 0.0; basis->c.y = 0.0; basis->c.z = 1.0;
printf("Reciprocal-space cell (nm^-1):\n");
basis_print(basis);
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);
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);
gsl_matrix_transpose(inv);
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);
printf("Real-space cell (nm):\n");
printf("%12.8f %12.8f %12.8f\n", ax*1e9, ay*1e9, az*1e9);
printf("%12.8f %12.8f %12.8f\n", bx*1e9, by*1e9, bz*1e9);
printf("%12.8f %12.8f %12.8f\n", cx*1e9, cy*1e9, cz*1e9);
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);
printf("Cell parameters:\n");
cell_print(&cell);
return cell;
}
|
