tascalc.py 6.38 KB
Newer Older
1
#
2
# calculates TAS angles from rlu
3
4
5
6
7
8
9
10
11
# @author Tobias Weber <tweber@ill.fr>
# @date 1-aug-18
# @license see 'LICENSE' file
#

import numpy as np
import numpy.linalg as la


Tobias WEBER's avatar
Tobias WEBER committed
12
use_scipy = False
Tobias WEBER's avatar
Tobias WEBER committed
13
14
15

# -----------------------------------------------------------------------------
# choose an a3 convention
Tobias WEBER's avatar
Tobias WEBER committed
16
#a3_offs = 0.			# for sics: a3 is angle between ki and orient1
Tobias WEBER's avatar
comment    
Tobias WEBER committed
17
#a3_offs = np.pi/2.		# for takin: Q along orient1 => a3:=a4/2
Tobias WEBER's avatar
Tobias WEBER committed
18
a3_offs = np.pi 		# for nomad: a3 is angle between ki and orient1 + 180 deg
Tobias WEBER's avatar
Tobias WEBER committed
19
20
21
22

def set_a3_offs(offs):
	global a3_offs
	a3_offs = offs
Tobias WEBER's avatar
Tobias WEBER committed
23
# -----------------------------------------------------------------------------
Tobias WEBER's avatar
Tobias WEBER committed
24
25


Tobias WEBER's avatar
Tobias WEBER committed
26
27
28
29
30
31
32
33
34
35
36
37
# -----------------------------------------------------------------------------
# rotate a vector around an axis using Rodrigues' formula
# see: https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula
def rotate(_axis, vec, phi):
	axis = _axis / la.norm(_axis)

	s = np.sin(phi)
	c = np.cos(phi)

	return c*vec + (1.-c)*np.dot(vec, axis)*axis + s*np.cross(axis, vec)


38
# get metric from crystal B matrix
Tobias WEBER's avatar
Tobias WEBER committed
39
# basis vectors are in the columns of B, i.e. the second index
40
def get_metric(B):
Tobias WEBER's avatar
Tobias WEBER committed
41
42
	#return np.einsum("ij,ik -> jk", B, B)
	return np.dot(np.transpose(B), B)
43
44


Tobias WEBER's avatar
Tobias WEBER committed
45
46
47
48
49
50
51
# cross product in fractional coordinates
def cross(a, b, B):
	# levi-civita in fractional coordinates
	def levi(i,j,k, B):
		M = np.array([B[:,i], B[:,j], B[:,k]])
		return la.det(M)

52
	metric_inv = la.inv(get_metric(B))
Tobias WEBER's avatar
Tobias WEBER committed
53
54
	eps = [[[ levi(i,j,k, B) for k in range(0,3) ] for j in range(0,3) ] for i in range(0,3) ]
	return np.einsum("ijk,j,k,li -> l", eps, a, b, metric_inv)
55
56
57
58
59


# dot product in fractional coordinates
def dot(a, b, metric):
	return np.dot(a, np.dot(metric, b))
Tobias WEBER's avatar
Tobias WEBER committed
60
61
62
63
64
65
66
67
68


# angle between peaks in fractional coordinates
def angle(a, b, metric):
	len_a = np.sqrt(dot(a, a, metric))
	len_b = np.sqrt(dot(b, b, metric))

	c = dot(a, b, metric) / (len_a * len_b)
	return np.arccos(c)
Tobias WEBER's avatar
Tobias WEBER committed
69
70
71
# -----------------------------------------------------------------------------


Tobias WEBER's avatar
Tobias WEBER committed
72
73
74
75
76
77
78
79
80
81
# -----------------------------------------------------------------------------
if use_scipy:
	import scipy as sp
	import scipy.constants as co

	hbar_in_meVs = co.Planck/co.elementary_charge*1000./2./np.pi
	E_to_k2 = 2.*co.neutron_mass/hbar_in_meVs**2. / co.elementary_charge*1000. * 1e-20
else:
	E_to_k2 = 0.482596406464	# calculated with scipy, using the formula above

82
k2_to_E = 1./E_to_k2
Tobias WEBER's avatar
Tobias WEBER committed
83
# -----------------------------------------------------------------------------
84
85


Tobias WEBER's avatar
Tobias WEBER committed
86
87
# -----------------------------------------------------------------------------
# mono (or ana) k  ->  A1 & A2 angles (or A5 & A6)
88
89
90
91
92
93
def get_a1a2(k, d):
	s = np.pi/(d*k)
	a1 = np.arcsin(s)
	return [a1, 2.*a1]


Tobias WEBER's avatar
Tobias WEBER committed
94
# a1 angle (or a5)  ->  mono (or ana) k
Tobias WEBER's avatar
Tobias WEBER committed
95
96
97
98
99
100
101
102
def get_monok(theta, d):
	s = np.sin(theta)
	k = np.pi/(d*s)
	return k
# -----------------------------------------------------------------------------


# -----------------------------------------------------------------------------
Tobias WEBER's avatar
Tobias WEBER committed
103
# scattering angle a4
104
105
106
107
108
def get_a4(ki, kf, Q):
	c = (ki**2. + kf**2. - Q**2.) / (2.*ki*kf)
	return np.arccos(c)


Tobias WEBER's avatar
Tobias WEBER committed
109
# get |Q| from ki, kf and a4
Tobias WEBER's avatar
Tobias WEBER committed
110
111
112
113
114
115
116
def get_Q(ki, kf, a4):
	c = np.cos(a4)
	return np.sqrt(ki**2. + kf**2. - c*(2.*ki*kf))
# -----------------------------------------------------------------------------


# -----------------------------------------------------------------------------
Tobias WEBER's avatar
Tobias WEBER committed
117
# angle enclosed by ki and Q
Tobias WEBER's avatar
Tobias WEBER committed
118
def get_psi(ki, kf, Q, sense=1.):
119
	c = (ki**2. + Q**2. - kf**2.) / (2.*ki*Q)
Tobias WEBER's avatar
Tobias WEBER committed
120
	return sense*np.arccos(c)
121
122


Tobias WEBER's avatar
Tobias WEBER committed
123
# crystallographic A matrix converting fractional to lab coordinates
Tobias WEBER's avatar
Tobias WEBER committed
124
# see: https://de.wikipedia.org/wiki/Fraktionelle_Koordinaten
125
126
127
128
129
130
131
132
133
134
def get_A(lattice, angles):
	cs = np.cos(angles)
	s2 = np.sin(angles[2])

	a = lattice[0] * np.array([1, 0, 0])
	b = lattice[1] * np.array([cs[2], s2, 0])
	c = lattice[2] * np.array([cs[1], \
		(cs[0]-cs[1]*cs[2]) / s2, \
		(np.sqrt(1. - np.dot(cs,cs) + 2.*cs[0]*cs[1]*cs[2])) / s2])

135
136
137
138
	# testing equality with own derivation
	#print((np.sqrt(1. - np.dot(cs,cs) + 2.*cs[0]*cs[1]*cs[2])) / s2)
	#print(np.sqrt(1. - cs[1]*cs[1] - ((cs[0] - cs[2]*cs[1])/s2)**2.))

Tobias WEBER's avatar
Tobias WEBER committed
139
	# the real-space basis vectors form the columns of the A matrix
140
141
142
	return np.transpose(np.array([a, b, c]))


Tobias WEBER's avatar
Tobias WEBER committed
143
# crystallographic B matrix converting rlu to 1/A
Tobias WEBER's avatar
Tobias WEBER committed
144
# the reciprocal-space basis vectors form the columns of the B matrix
145
146
147
148
149
150
def get_B(lattice, angles):
	A = get_A(lattice, angles)
	B = 2.*np.pi * np.transpose(la.inv(A))
	return B


Tobias WEBER's avatar
Tobias WEBER committed
151
152
153
154
155
156
157
158
159
160
161
162
163
164
# UB orientation matrix
def get_UB(B, orient1_rlu, orient2_rlu, orientup_rlu):
	orient1_invA = np.dot(B, orient1_rlu)
	orient2_invA = np.dot(B, orient2_rlu)
	orientup_invA = np.dot(B, orientup_rlu)

	orient1_invA = orient1_invA / la.norm(orient1_invA)
	orient2_invA = orient2_invA / la.norm(orient2_invA)
	orientup_invA = orientup_invA / la.norm(orientup_invA)

	U_invA = np.array([orient1_invA, orient2_invA, orientup_invA])
	UB = np.dot(U_invA, B)
	return UB

Tobias WEBER's avatar
Tobias WEBER committed
165

Tobias WEBER's avatar
Tobias WEBER committed
166
# a3 & a4 angles
Tobias WEBER's avatar
Tobias WEBER committed
167
def get_a3a4(ki, kf, Q_rlu, orient_rlu, orient_up_rlu, B, sense_sample=1.):
168
	metric = get_metric(B)
169

Tobias WEBER's avatar
Tobias WEBER committed
170
	# angle xi between Q and orientation reflex
Tobias WEBER's avatar
Tobias WEBER committed
171
	xi = angle(Q_rlu, orient_rlu, metric)
172
173
174
175

	# sign of xi
	if dot(cross(orient_rlu, Q_rlu, B), orient_up_rlu, metric) < 0.:
		xi = -xi
176

Tobias WEBER's avatar
Tobias WEBER committed
177
	# length of Q
Tobias WEBER's avatar
Tobias WEBER committed
178
179
	Qlen = np.sqrt(dot(Q_rlu, Q_rlu, metric))

Tobias WEBER's avatar
Tobias WEBER committed
180
181
182
183
	# distance to plane
	up_len = np.sqrt(dot(orient_up_rlu, orient_up_rlu, metric))
	dist_Q_plane = dot(Q_rlu, orient_up_rlu, metric) / up_len

Tobias WEBER's avatar
Tobias WEBER committed
184
	# angle psi enclosed by ki and Q
Tobias WEBER's avatar
Tobias WEBER committed
185
	psi = get_psi(ki, kf, Qlen, sense_sample)
186

Tobias WEBER's avatar
Tobias WEBER committed
187
	a3 = - psi - xi + a3_offs
188
189
	a4 = get_a4(ki, kf, Qlen)

190
	#print("xi = " + str(xi/np.pi*180.) + ", psi = " + str(psi/np.pi*180.) + ", offs = " + str(a3_offs/np.pi*180.))
Tobias WEBER's avatar
Tobias WEBER committed
191
	return [a3, a4, dist_Q_plane]
192
193


Tobias WEBER's avatar
Tobias WEBER committed
194
def get_hkl(ki, kf, a3, Qlen, orient_rlu, orient_up_rlu, B, sense_sample=1.):
Tobias WEBER's avatar
Tobias WEBER committed
195
196
	B_inv = la.inv(B)

Tobias WEBER's avatar
Tobias WEBER committed
197
	# angle enclosed by ki and Q
Tobias WEBER's avatar
Tobias WEBER committed
198
	psi = get_psi(ki, kf, Qlen, sense_sample)
Tobias WEBER's avatar
Tobias WEBER committed
199

Tobias WEBER's avatar
Tobias WEBER committed
200
	# angle between Q and orientation reflex
Tobias WEBER's avatar
Tobias WEBER committed
201
202
	xi = - a3 + a3_offs - psi

Tobias WEBER's avatar
Tobias WEBER committed
203
204
	Q_lab = rotate(np.dot(B, orient_up_rlu), np.dot(B, orient_rlu*Qlen), xi)
	Q_lab *= Qlen / la.norm(Q_lab)
Tobias WEBER's avatar
Tobias WEBER committed
205
206
207
208
209
210
211
	Q_rlu = np.dot(B_inv, Q_lab)

	return Q_rlu
# -----------------------------------------------------------------------------


# -----------------------------------------------------------------------------
Tobias WEBER's avatar
Tobias WEBER committed
212
# get ki from kf and energy transfer
213
214
215
216
def get_ki(kf, E):
	return np.sqrt(kf**2. + E_to_k2*E)


Tobias WEBER's avatar
Tobias WEBER committed
217
# get kf from ki and energy transfer
218
219
220
221
def get_kf(ki, E):
	return np.sqrt(ki**2. - E_to_k2*E)


Tobias WEBER's avatar
Tobias WEBER committed
222
# get energy transfer from ki and kf
Tobias WEBER's avatar
Tobias WEBER committed
223
224
225
def get_E(ki, kf):
	return (ki**2. - kf**2.) / E_to_k2
# -----------------------------------------------------------------------------