tascalc.py 6.38 KB
 Tobias WEBER committed Aug 02, 2018 1 ``````# `````` Tobias WEBER committed Jun 29, 2019 2 ``````# calculates TAS angles from rlu `````` Tobias WEBER committed Aug 02, 2018 3 4 5 6 7 8 9 10 11 ``````# @author Tobias Weber # @date 1-aug-18 # @license see 'LICENSE' file # import numpy as np import numpy.linalg as la `````` Tobias WEBER committed Aug 04, 2018 12 ``````use_scipy = False `````` Tobias WEBER committed Oct 22, 2018 13 14 15 `````` # ----------------------------------------------------------------------------- # choose an a3 convention `````` Tobias WEBER committed Oct 25, 2018 16 ``````#a3_offs = 0. # for sics: a3 is angle between ki and orient1 `````` Tobias WEBER committed Oct 22, 2018 17 ``````#a3_offs = np.pi/2. # for takin: Q along orient1 => a3:=a4/2 `````` Tobias WEBER committed Oct 25, 2018 18 ``````a3_offs = np.pi # for nomad: a3 is angle between ki and orient1 + 180 deg `````` Tobias WEBER committed Oct 24, 2018 19 20 21 22 `````` def set_a3_offs(offs): global a3_offs a3_offs = offs `````` Tobias WEBER committed Oct 22, 2018 23 ``````# ----------------------------------------------------------------------------- `````` Tobias WEBER committed Aug 04, 2018 24 25 `````` `````` Tobias WEBER committed Aug 04, 2018 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) `````` Tobias WEBER committed Aug 04, 2018 38 ``````# get metric from crystal B matrix `````` Tobias WEBER committed Feb 26, 2019 39 ``````# basis vectors are in the columns of B, i.e. the second index `````` Tobias WEBER committed Aug 04, 2018 40 ``````def get_metric(B): `````` Tobias WEBER committed Mar 04, 2019 41 42 `````` #return np.einsum("ij,ik -> jk", B, B) return np.dot(np.transpose(B), B) `````` Tobias WEBER committed Aug 04, 2018 43 44 `````` `````` Tobias WEBER committed Aug 04, 2018 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) `````` Tobias WEBER committed Aug 04, 2018 52 `````` metric_inv = la.inv(get_metric(B)) `````` Tobias WEBER committed Aug 04, 2018 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) `````` Tobias WEBER committed Aug 04, 2018 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 committed Sep 04, 2018 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 committed Aug 04, 2018 69 70 71 ``````# ----------------------------------------------------------------------------- `````` Tobias WEBER committed Aug 04, 2018 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 `````` Tobias WEBER committed Jul 31, 2019 82 ``````k2_to_E = 1./E_to_k2 `````` Tobias WEBER committed Aug 04, 2018 83 ``````# ----------------------------------------------------------------------------- `````` Tobias WEBER committed Aug 02, 2018 84 85 `````` `````` Tobias WEBER committed Aug 04, 2018 86 87 ``````# ----------------------------------------------------------------------------- # mono (or ana) k -> A1 & A2 angles (or A5 & A6) `````` Tobias WEBER committed Aug 02, 2018 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 committed Aug 04, 2018 94 ``````# a1 angle (or a5) -> mono (or ana) k `````` Tobias WEBER committed Aug 04, 2018 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 committed Aug 04, 2018 103 ``````# scattering angle a4 `````` Tobias WEBER committed Aug 02, 2018 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 committed Aug 04, 2018 109 ``````# get |Q| from ki, kf and a4 `````` Tobias WEBER committed Aug 04, 2018 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 committed Aug 04, 2018 117 ``````# angle enclosed by ki and Q `````` Tobias WEBER committed Feb 25, 2019 118 ``````def get_psi(ki, kf, Q, sense=1.): `````` Tobias WEBER committed Aug 02, 2018 119 `````` c = (ki**2. + Q**2. - kf**2.) / (2.*ki*Q) `````` Tobias WEBER committed Feb 25, 2019 120 `````` return sense*np.arccos(c) `````` Tobias WEBER committed Aug 02, 2018 121 122 `````` `````` Tobias WEBER committed Aug 04, 2018 123 ``````# crystallographic A matrix converting fractional to lab coordinates `````` Tobias WEBER committed Aug 04, 2018 124 ``````# see: https://de.wikipedia.org/wiki/Fraktionelle_Koordinaten `````` Tobias WEBER committed Aug 02, 2018 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]) `````` Tobias WEBER committed Feb 28, 2019 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 committed Feb 26, 2019 139 `````` # the real-space basis vectors form the columns of the A matrix `````` Tobias WEBER committed Aug 02, 2018 140 141 142 `````` return np.transpose(np.array([a, b, c])) `````` Tobias WEBER committed Aug 04, 2018 143 ``````# crystallographic B matrix converting rlu to 1/A `````` Tobias WEBER committed Feb 26, 2019 144 ``````# the reciprocal-space basis vectors form the columns of the B matrix `````` Tobias WEBER committed Aug 02, 2018 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 committed Oct 26, 2018 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 committed Feb 25, 2019 165 `````` `````` Tobias WEBER committed Aug 04, 2018 166 ``````# a3 & a4 angles `````` Tobias WEBER committed Feb 25, 2019 167 ``````def get_a3a4(ki, kf, Q_rlu, orient_rlu, orient_up_rlu, B, sense_sample=1.): `````` Tobias WEBER committed Aug 04, 2018 168 `````` metric = get_metric(B) `````` Tobias WEBER committed Aug 02, 2018 169 `````` `````` Tobias WEBER committed Aug 04, 2018 170 `````` # angle xi between Q and orientation reflex `````` Tobias WEBER committed Sep 04, 2018 171 `````` xi = angle(Q_rlu, orient_rlu, metric) `````` Tobias WEBER committed Aug 04, 2018 172 173 174 175 `````` # sign of xi if dot(cross(orient_rlu, Q_rlu, B), orient_up_rlu, metric) < 0.: xi = -xi `````` Tobias WEBER committed Aug 02, 2018 176 `````` `````` Tobias WEBER committed Oct 25, 2018 177 `````` # length of Q `````` Tobias WEBER committed Sep 04, 2018 178 179 `````` Qlen = np.sqrt(dot(Q_rlu, Q_rlu, metric)) `````` Tobias WEBER committed Oct 25, 2018 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 committed Aug 04, 2018 184 `````` # angle psi enclosed by ki and Q `````` Tobias WEBER committed Feb 25, 2019 185 `````` psi = get_psi(ki, kf, Qlen, sense_sample) `````` Tobias WEBER committed Aug 02, 2018 186 `````` `````` Tobias WEBER committed Aug 04, 2018 187 `````` a3 = - psi - xi + a3_offs `````` Tobias WEBER committed Aug 02, 2018 188 189 `````` a4 = get_a4(ki, kf, Qlen) `````` Tobias WEBER committed Jul 07, 2019 190 `````` #print("xi = " + str(xi/np.pi*180.) + ", psi = " + str(psi/np.pi*180.) + ", offs = " + str(a3_offs/np.pi*180.)) `````` Tobias WEBER committed Oct 25, 2018 191 `````` return [a3, a4, dist_Q_plane] `````` Tobias WEBER committed Aug 02, 2018 192 193 `````` `````` Tobias WEBER committed Feb 25, 2019 194 ``````def get_hkl(ki, kf, a3, Qlen, orient_rlu, orient_up_rlu, B, sense_sample=1.): `````` Tobias WEBER committed Aug 04, 2018 195 196 `````` B_inv = la.inv(B) `````` Tobias WEBER committed Aug 04, 2018 197 `````` # angle enclosed by ki and Q `````` Tobias WEBER committed Feb 25, 2019 198 `````` psi = get_psi(ki, kf, Qlen, sense_sample) `````` Tobias WEBER committed Aug 04, 2018 199 `````` `````` Tobias WEBER committed Aug 04, 2018 200 `````` # angle between Q and orientation reflex `````` Tobias WEBER committed Aug 04, 2018 201 202 `````` xi = - a3 + a3_offs - psi `````` Tobias WEBER committed Aug 04, 2018 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 committed Aug 04, 2018 205 206 207 208 209 210 211 `````` Q_rlu = np.dot(B_inv, Q_lab) return Q_rlu # ----------------------------------------------------------------------------- # ----------------------------------------------------------------------------- `````` Tobias WEBER committed Aug 04, 2018 212 ``````# get ki from kf and energy transfer `````` Tobias WEBER committed Aug 02, 2018 213 214 215 216 ``````def get_ki(kf, E): return np.sqrt(kf**2. + E_to_k2*E) `````` Tobias WEBER committed Aug 04, 2018 217 ``````# get kf from ki and energy transfer `````` Tobias WEBER committed Aug 02, 2018 218 219 220 221 ``````def get_kf(ki, E): return np.sqrt(ki**2. - E_to_k2*E) `````` Tobias WEBER committed Aug 04, 2018 222 ``````# get energy transfer from ki and kf `````` Tobias WEBER committed Aug 04, 2018 223 224 225 ``````def get_E(ki, kf): return (ki**2. - kf**2.) / E_to_k2 # -----------------------------------------------------------------------------``````