ADD: Goniometer class

BETA FEATURE: Ability to defined a goniometer to create a UB matrix for crystal alignment

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,25 +55,109 @@ 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
     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.]