-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathdiffraction.py
256 lines (179 loc) · 8.65 KB
/
diffraction.py
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
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
"""
The following code follows the treatment in "Theory of X-Ray Diffraction in Crystals" by William H. Zachariasen
"""
import numpy as np
from scipy import constants
class Crystal(object):
"""
This class stores crystal data read from a data file (e.g. Si.txt, Ge.txt, etc.)
"""
# TODO: asymmetry angle not needed here, used in Diffraction
def __init__(self, crystal_type, asymmetry_angle):
"""
Initializer.
:param crystal_type: <class 'str'> string specifying which crystal is being used
:param asymmetry_angle: <class 'float'> angle in degrees
"""
self._crystal_type = crystal_type
self._asymmetry_angle = asymmetry_angle * np.pi / 180 # converting into radians
"""
The files are structured in the following way (each row has 11 elements):
(0, 0)-(0, 4) -> 5x{a_i} | (0, 5) -> c | (0, 6)-(0, 10) -> 5x{b_i} goal:f0 calculation
(1, 0)-(1, 5) -> a b c alpha beta gamma | (1, 6)-(1, 10) -> 5x0 goal:unit cell volume calculation
(2-9, 0-10) -> atomic_number fraction x y z 6x0 goal:structure factor calculation
"""
self._data_array = np.loadtxt("{}.txt".format(self._crystal_type), dtype='float')
self._a = [] # parameters a_i (x5)
self._b = [] # parameters b_i (x5)
self._unit_cell = [] # unit cell sides (x3) and angles (x3)
self._atom_coordinates = np.zeros((8, 3)) # each row holds the x, y, z coordinates of one atom in the unit cell
# there are 8 atoms in the silicon unit cell
for i in np.arange(self._data_array.shape[0]):
for j in np.arange(self._data_array.shape[1]):
if i == 0:
if j < 5:
self._a.append(self._data_array[(i, j)])
elif 5 < j < self._data_array.shape[1]:
self._b.append(self._data_array[(i, j)])
elif j == 5:
self._c = self._data_array[(i, j)]
elif i == 1 and j < 6:
self._unit_cell.append(self._data_array[(i, j)])
else:
if j in [2, 3, 4]:
self._atom_coordinates[(i - 2, j - 2)] = self._data_array[(i, j)]
def unit_cell_volume(self):
"""
calculates the unit cell volume.
Conventions: alpha is the angle between b and c, beta between a and c, gamma between a and b,
with a, b being on the horizontal plane.
Volume = a * b * c * sin(alpha) * sin(beta) * sin(gamma)
:return: unit cell volume in Angstroms^3
"""
a = self._unit_cell[0]
b = self._unit_cell[1]
c = self._unit_cell[2]
alpha = self._unit_cell[3] * np.pi / 180 # converting into radians
beta = self._unit_cell[4] * np.pi / 180 # converting into radians
gamma = self._unit_cell[5] * np.pi / 180 # converting into radians
return a * b * c * np.sin(alpha) * np.sin(beta) * np.sin(gamma) # Angstroms ^ 3
class Diffraction(Crystal):
"""
This class defines a diffraction geometry for the crystal and the incoming x-ray beam.
In particular, this class is initialized by specifying the wavelength of the beam and then this is used to calculate
the Bragg angle, the b parameter and the structure factors.
"""
# todo: not here: global variables?? use "import scipy.constants as constants" and use "constants.c" etc
e = constants.elementary_charge # charge of the proton
c = constants.speed_of_light # speed of light in vacuum
m = constants.electron_mass # mass of the electron
h = constants.Planck # Planck's constant
def __init__(self, crystal_type, asymmetry_angle, miller_h, miller_k, miller_l, wavelength, polarization):
"""
:param miller_h: < class 'int'> Miller index H
:param miller_k: < class 'int'> Miller index K
:param miller_l: < class 'int'> Miller index L
:param wavelength: < class 'float'> wavelength of the x-ray beam in keV
:param polarization: < class 'bool'> True = normal polarization (s), False = parallel polarization (p)
"""
super(Diffraction, self).__init__(crystal_type, asymmetry_angle)
self._wavelength = self.h * self.c * 1e+7 / (self.e * wavelength) # conversion keV -> Angstrom
self._polarization = polarization
self._miller_h = miller_h
self._miller_k = miller_k
self._miller_l = miller_l
def d_spacing(self):
"""
calculates the spacing between the Miller planes in the special case of a cubic unit cell
:return: d spacing in Angstroms
"""
return self._unit_cell[0] / np.sqrt(self._miller_h ** 2 + self._miller_k ** 2 + self._miller_l ** 2)
def bragg_angle(self):
"""
calculates the Bragg angle using Bragg's law: lambda = 2 * d_spacing * sin(theta_Bragg) (Zachariasen [3.5])
:return: Bragg angle in radians
"""
d_spacing = self.d_spacing() # in Angstroms
return np.arcsin(0.5 * self._wavelength / d_spacing)
def parameter_f_0(self):
"""
calculates f0 following the 5-gaussian interpolation by D. Waasmaier & A. Kirfel
:return: atomic scattering power, f0 @ specified wavelength
"""
f = self._c
s = np.sin(self.bragg_angle()) / self._wavelength # in Angstrom ^ -1
for i in np.arange(5):
f += self._a[i] * np.exp(-self._b[i] * s ** 2)
return f
def parameter_b(self):
"""
calculates b = gamma_0 / gamma_h (Zachariasen [3.115])
:return: b
"""
gamma_h = - np.sin(self.bragg_angle())
gamma_0 = np.sin(self.bragg_angle() + self._asymmetry_angle)
return gamma_0 / gamma_h
def parameter_k(self):
"""
k = 1 for normal polarization; k = |cos(2*theta_Bragg)| for parallel polarization (Zachariasen [3.141])
:return: k
"""
if self._polarization:
return 1
if not self._polarization:
return np.absolute(np.cos(2*self.bragg_angle()))
def parameter_psi_0(self):
"""
psi_0 = num_atoms_in_unit_cell * f_0 (Zachariasen [3.95])
:return: psi_0
"""
f_0 = self.parameter_f_0()
v = self.unit_cell_volume() * 1e-30 # Angstrom ^3 -> m ^3
wavelength = self._wavelength * 1e-10 # Angstrom -> metre
e_radius = 2.8179403267e-15 # cgs units
psi_0 = - e_radius * wavelength ** 2 / (v * np.pi) * f_0
return (self._data_array.shape[0] - 2) * psi_0 # this only holds for crystals with atoms all of the same kind
def parameter_f_h(self):
"""
F_H = f_0 * SUM(i) { exp[ (i*2*pi) * (h*a_i + k*b_i + l*c_i) ] }
where a_i, b_i, c_i are the coordinates of the i-th atom in the unit cell (Zachariasen [3.52])
and h, k, l are the Miller indices.
:return: F_H
"""
f_h = 0 + 0j
for i in np.arange(self._data_array.shape[0] - 2):
j = 0
f_h += np.exp(1j * 2 * np.pi * (self._miller_h * self._atom_coordinates[(i, j)]
+ self._miller_k * self._atom_coordinates[(i, j + 1)]
+ self._miller_l * self._atom_coordinates[(i, j + 2)]))
return f_h * self.parameter_f_0()
def parameter_psi_h(self):
"""
calculates psi_H starting from F_H, the wavelength and the unit cell volume (Zachariasen [3.95])
:return: psi_H
"""
v = self.unit_cell_volume() * 1e-30 # Angstrom ^3 -> m ^3
wavelength = self._wavelength * 1e-10 # Angstrom -> metre
f_h = self.parameter_f_h()
e_radius = 2.8179403267e-15 # cgs units
return - e_radius * (wavelength ** 2 / (v * np.pi)) * f_h
def darwin_profile_theta(self,theta_diff):
num = theta_diff.size
theta_bragg = self.bragg_angle()
alpha = 2 * theta_diff * 1e-6 * np.sin(2 * theta_bragg) # Zachariasen [3.116]
b = self.parameter_b()
k = self.parameter_k()
psi_0 = np.ones(num) * self.parameter_psi_0()
psi_h = self.parameter_psi_h()
y = ((0.5 * (1-b) * psi_0) + (0.5 * b * alpha)) / (np.sqrt(np.absolute(b)) * k * np.absolute(psi_h))
reflectivity = np.ones(num)
for i in range(num): # Zachariasen [3.192b]
if y[i] < -1:
reflectivity[i] = (-y[i] - np.sqrt(y[i] ** 2 - 1)) ** 2
elif y[i] > 1:
reflectivity[i] = (y[i] - np.sqrt(y[i] ** 2 - 1)) ** 2
else:
y[i] = 1
return reflectivity
def darwin_profile_y(self,y_array):
pass