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
|
/*
* crystal.c
*
* Crystallographic stuff
*
* (c) 2009 Thomas White <taw27@cam.ac.uk>
*
* Triclinator - solve nasty triclinic unit cells
*
*/
#include <math.h>
#include "crystal.h"
#include "util.h"
/* Some attempt at preserving sanity */
#define S11 s11(cell)
#define S22 s22(cell)
#define S33 s33(cell)
#define S12 s12(cell)
#define S23 s23(cell)
#define S13 s13(cell)
static double s11(Cell cell)
{
return cell.b*cell.b*cell.c*cell.c*sin(cell.al)*sin(cell.al);
}
static double s22(Cell cell)
{
return cell.a*cell.a*cell.c*cell.c*sin(cell.be)*sin(cell.be);
}
static double s33(Cell cell)
{
return cell.a*cell.a*cell.b*cell.b*sin(cell.ga)*sin(cell.ga);
}
static double s12(Cell cell)
{
return cell.a*cell.b*cell.c*cell.c*(cos(cell.al)*cos(cell.be) - cos(cell.ga));
}
static double s23(Cell cell)
{
return cell.a*cell.a*cell.b*cell.c*(cos(cell.be)*cos(cell.ga) - cos(cell.al));
}
static double s13(Cell cell)
{
return cell.a*cell.b*cell.b*cell.c*(cos(cell.ga)*cos(cell.al) - cos(cell.be));
}
static double volume(Cell cell)
{
return cell.a*cell.b*cell.c*sqrt(1 - cos(cell.al)*cos(cell.al)
- cos(cell.be)*cos(cell.be)
- cos(cell.ga)*cos(cell.ga)
+ 2*cos(cell.al)*cos(cell.be)*cos(cell.ga) );
}
double dspacing(Cell cell, double h, double k, double l)
{
double sum_S_terms, one_over_Vsq, one_over_LHS;
sum_S_terms = S11*h*h + S22*k*k + S33*l*l
+ 2*S12*h*k + 2*S23*k*l + 2*S13*h*l;
one_over_Vsq = 1 / pow(volume(cell), 2);
one_over_LHS = 1/( one_over_Vsq*sum_S_terms );
return sqrt(one_over_LHS);
}
double plane_angle(Cell cell, double h1, double k1, double l1,
double h2, double k2, double l2)
{
double sum_S_terms, d1, d2, dd_over_Vsq;
sum_S_terms = S11*h1*h2 + S22*k1*k2 * S33*l1*l2
+ (S23*(k1*l2 + k2*l1))
+ (S13*(l1*h2 + l2*h1))
+ (S12*(h1*k2 + h2*k1));
d1 = dspacing(cell, h1, k1, l1);
d2 = dspacing(cell, h2, k2, l2);
dd_over_Vsq = d1*d2 / pow(volume(cell), 2);
/* Bodge to put the axes in the right place */
return M_PI - acos(dd_over_Vsq * sum_S_terms);
}
/* Return what the measurement 'val' would have been if the cell were 'cell' */
double crystal_calc(MVal val, Cell cell)
{
if ( is_dspacing(val) ) {
return dspacing(cell, val.h1, val.k1, val.l1);
} else {
return plane_angle(cell, val.h1, val.k1, val.l1,
val.h2, val.k2, val.l2);
}
}
|
