diff --git a/neutronpy/lattice.py b/neutronpy/lattice.py index cd40966..95cd9c0 100644 --- a/neutronpy/lattice.py +++ b/neutronpy/lattice.py @@ -21,12 +21,31 @@ class Lattice(object): Attributes ---------- + a + b + c + astar + bstar + cstar + alpha + beta + gamma + alphastar + betastar + gammastar + alphastar_rad + betastar_rad + gammastar_rad abg_rad + reciprocal_abc + reciprocal_abg + reciprocal_abg_rad lattice_type volume reciprocal_volume G Gstar + B Methods ------- @@ -36,17 +55,101 @@ class Lattice(object): get_angle_between_planes ''' - + def __init__(self, abc, abg): self.abc = np.array(abc) self.abg = np.array(abg) + @property + def a(self): + return self.abc[0] + + @property + def b(self): + return self.abc[1] + + @property + def c(self): + return self.abc[2] + + @property + def alpha(self): + return self.abg[0] + + @property + def beta(self): + return self.abg[1] + + @property + def gamma(self): + return self.abg[2] + + @property + def alpha_rad(self): + return self.abg_rad[0] + + @property + def beta_rad(self): + return self.abg_rad[1] + + @property + def gamma_rad(self): + return self.abg_rad[2] + + @property + def astar(self): + return self.b * self.c * np.sin(self.alpha_rad) / self.volume + + @property + def bstar(self): + return self.a * self.c * np.sin(self.beta_rad) / self.volume + + @property + def cstar(self): + return self.a * self.b * np.sin(self.gamma) / self.volume + + @property + def alphastar_rad(self): + return np.arccos((np.cos(self.beta_rad) * np.cos(self.gamma_rad) - np.cos(self.alpha_rad)) / (np.sin(self.beta_rad) * np.sin(self.gamma_rad))) + + @property + def betastar_rad(self): + return np.arccos((np.cos(self.alpha_rad) * np.cos(self.gamma_rad) - np.cos(self.beta_rad)) / (np.sin(self.alpha_rad) * np.sin(self.gamma_rad))) + + @property + def gammastar_rad(self): + return np.arccos((np.cos(self.alpha_rad) * np.cos(self.beta_rad) - np.cos(self.gamma_rad)) / (np.sin(self.alpha_rad) * np.sin(self.beta_rad))) + + @property + def alphastar(self): + return np.rad2deg(self.alphastar_rad) + + @property + def betastar(self): + return np.rad2deg(self.betastar_rad) + + @property + def gammastar(self): + return np.rad2deg(self.gammastar_rad) + + @property + def reciprocal_abc(self): + return [self.astar, self.bstar, self.cstar] + + @property + def reciprocal_abg(self): + return [self.alphastar, self.betastar, self.gammastar] + + @property + def reciprocal_abg_rad(self): + return [self.alphastar_rad, self.betastar_rad, self.gammastar_rad] + @property def abg_rad(self): r'''Lattice angles in radians ''' - + return np.deg2rad(self.abg) @property @@ -54,7 +157,7 @@ def lattice_type(self): r'''Type of lattice determined by the provided lattice constants and angles ''' - + if len(np.unique(self.abc)) == 3 and len(np.unique(self.abg)) == 3: return 'triclinic' elif len(np.unique(self.abc)) == 3 and self.abg[1] != 90 and np.all(self.abg[:3:2] == 90): @@ -107,6 +210,17 @@ def Gstar(self): return np.linalg.inv(self.G) + + @property + def B(self): + r'''Cartesian basis matrix in reciprocal units such that B*B.T = Gstar + + ''' + + return np.matrix([[self.astar, self.bstar * np.cos(self.gammastar_rad), self.cstar * np.cos(self.betastar_rad)], + [0, self.bstar * np.sin(self.gammastar_rad), -self.cstar * np.sin(self.betastar_rad) * np.cos(self.alpha_rad)], + [0, 0, 1. / self.c]]) + def get_d_spacing(self, hkl): u'''Returns the d-spacing of a given reciprocal lattice vector. @@ -143,12 +257,12 @@ def get_angle_between_planes(self, v1, v2): The angle between v1 and v2 in degrees ''' - + v1, v2 = np.array(v1), np.array(v2) - return float(np.rad2deg(np.arccos(np.dot(np.dot(v1, self.Gstar), v2) / - np.sqrt(np.dot(np.dot(v1, self.Gstar), v1)) / - np.sqrt(np.dot(np.dot(v2, self.Gstar), v2))))) + return float(np.rad2deg(np.arccos(np.inner(np.inner(v1, self.Gstar), v2) / + np.sqrt(np.inner(np.inner(v1, self.Gstar), v1)) / + np.sqrt(np.inner(np.inner(v2, self.Gstar), v2))))) def get_two_theta(self, hkl, wavelength): u'''Returns the detector angle 2\U0001D703 for a given reciprocal @@ -168,7 +282,7 @@ def get_two_theta(self, hkl, wavelength): The angle of the detector 2\U0001D703 in degrees ''' - + return 2 * np.rad2deg(np.arcsin(wavelength / 2 / self.get_d_spacing(hkl))) def get_q(self, hkl): @@ -187,10 +301,98 @@ def get_q(self, hkl): \u212B\ :sup:`-1` ''' - + return 2 * np.pi / self.get_d_spacing(hkl) +class Goniometer(object): + r'''Defines a goniometer + + ''' + + def __init__(self, u, theta_u, v, theta_v, sgu, sgl, omega=0): + self.u = u + self.theta_u = theta_u + + self.v = v + self.theta_v = theta_v + + self.sgu = sgu + self.sgl = sgl + + self.omega = 0 + + @property + def omega_rad(self): + return self.omega_rad + + @property + def sgu_rad(self): + return np.deg2rad(self.sgu) + + @property + def sgl_rad(self): + return np.deg2rad(self.sgl) + + @property + def theta_rad(self): + return np.arctan((self.ki - self.kf * np.cos(self.phi)) / (self.kf * np.sin(self.phi))) + + @property + def theta(self): + pass + + @property + def N(self): + return np.matrix([[1, 0, 0], + [0, np.cos(self.sgu_rad), -np.sin(self.sgu_rad)], + [0, np.sin(self.sgu_rad), np.cos(self.sgu_rad)]]) + + @property + def M(self): + return np.matrix([[np.cos(self.sgl_rad), 0, np.sin(self.sgl_rad)], + [0, 1, 0], + [-np.sin(self.sgl_rad), 0, np.cos(self.sgl_rad)]]) + + @property + def Omega(self): + return np.matrix([[np.cos(self.omega_rad), -np.sin(self.omega_rad), 0], + [np.sin(self.omega_rad), np.cos(self.omega_rad), 0], + [0, 0, 1]]) + + @property + def Theta(self): + return np.matrix([[np.cos(self.theta_rad), -np.sin(self.theta_rad), 0], + [np.sin(self.theta_rad), np.cos(self.theta_rad), 0], + [0, 0, 1]]) + + @property + def T_c(self): + return np.matrix([self.u, self.v, np.cross(self.u, self.v)]).T + + @property + def T_phi(self): + return np.matrix([self.u_phi(np.deg2rad(self.theta_u), self.sgl_rad, self.sgu_rad), + self.u_phi(np.deg2rad(self.theta_v), self.sgl_rad, self.sgu_rad), + self.u_phi(np.deg2rad(0), np.deg2rad(90), np.deg2rad(0))]).T + + @property + def R(self): + return self.Omega * self.M * self.N + + @property + def U(self): + r'''Defines an orientation matrix based on supplied goniometer angles + + ''' + return self.T_phi * np.linalg.inv(self.T_c) + + def u_phi(self, omega, chi, phi): + return [np.cos(omega) * np.cos(chi) * np.cos(phi) - np.sin(omega) * np.sin(phi), + np.cos(omega) * np.cos(chi) * np.sin(phi) + np.sin(omega) * np.cos(phi), + np.cos(omega) * np.sin(chi)] + + if __name__ == "__main__": abc = [6.3496, 6.3496, 6.3496] abg = [90., 90., 90.]