<commit_before>'''
The microstructure module provide elementary classes to describe a 
crystallographic granular microstructure such as mostly present in 
metallic materials.

It contains several classes which are used to describe a microstructure 
composed of several grains, each one having its own crystallographic 
orientation:

 * :py:class:`~pymicro.crystal.microstructure.Microstructure`
 * :py:class:`~pymicro.crystal.microstructure.Grain`
 * :py:class:`~pymicro.crystal.microstructure.Orientation`
'''
import numpy as np
import os, vtk
from matplotlib import pyplot as plt, colors, cm
from xml.dom.minidom import Document, parse

class Orientation:
  '''Crystallographic orientation class.

  This follows the passive rotation definition which means that it brings 
  the sample coordinate system into coincidence with the crystal coordinate 
  system. Then one may express a vector :math:`V_c` in the crystal coordinate system 
  from the vector in the sample coordinate system :math:`V_s` by:
  
  .. math::

    V_c = B.V_s
  
  and inversely (because :math:`B^{-1}=B^T`):
  
  .. math::

    V_s = B^T.V_c
  '''

  def __init__(self, matrix):
    '''Initialization from the 9 components of the orientation matrix.'''
    g = np.array(matrix, dtype=np.float64).reshape((3, 3))
    self._matrix = g
    self.euler = Orientation.OrientationMatrix2Euler(g)
    self.rod = Orientation.OrientationMatrix2Rodrigues(g)

  def orientation_matrix(self):
    '''Returns the orientation matrix in the form of a 3x3 numpy array.'''
    return self._matrix

  def __repr__(self):
    '''Provide a string representation of the class.'''
    s = 'Crystal Orientation'
    s += '\norientation matrix = %s' % self._matrix.view()
    s += '\nEuler angles (degrees) = (%8.3f,%8.3f,%8.3f)' % (self.phi1(), self.Phi(), self.phi2())
    s += '\nRodrigues vector = %s' % self.OrientationMatrix2Rodrigues(self._matrix)
    return s

  @staticmethod
  def cube():
    '''Create the particular crystal orienation called Cube and which 
    corresponds to euler angle (0,0,0).'''
    return Orientation.from_euler((0., 0., 0.))

  @staticmethod
  def brass():
    '''Create the particular crystal orienation called Brass and which 
    corresponds to euler angle (35.264,45,0).'''
    return Orientation.from_euler((35.264, 45., 0.))

  @staticmethod
  def copper():
    '''Create the particular crystal orienation called Copper and which 
    corresponds to euler angle (90,35.264,45).'''
    return Orientation.from_euler((90., 35.264, 45.))

  @staticmethod
  def s3():
    '''Create the particular crystal orienation called S3 and which 
    corresponds to euler angle (59,37,63).'''
    return Orientation.from_euler((58.980, 36.699, 63.435))

  @staticmethod
  def goss():
    '''Create the particular crystal orienation called Goss and which 
    corresponds to euler angle (0,45,0).'''
    return Orientation.from_euler((0., 45., 0.))

  @staticmethod
  def shear():
    '''Create the particular crystal orienation called shear and which 
    corresponds to euler angle (45,0,0).'''
    return Orientation.from_euler((45., 0., 0.))

  @staticmethod
  def misorientation_MacKenzie(psi):
    '''Return the fraction of the misorientations corresponding to the 
    given :math:`\\psi` angle in the refenrece solution derived By MacKenzie in 
    his 1958 paper.
    
    :param psi: the misorientation angle in radians.
    :returns: the value in the cummulative distribution corresponding to psi.
    '''
    from math import sqrt, sin, cos, tan, pi, acos
    psidg = 180*psi/pi
    if psidg >= 0 and psidg <= 45:
      p = 2./15*(1 - cos(psi))
    elif psidg > 45 and psidg <= 60:
      p = 2./15*(3*(sqrt(2)-1)*sin(psi) - 2*(1 - cos(psi)))
    elif psidg > 60 and psidg <= 60.72:
      p = 2./15*((3*(sqrt(2)-1) + 4./sqrt(3))*sin(psi) - 6.*(1 - cos(psi)))
    elif psidg > 60.72 and psidg <= 62.8:
      X = (sqrt(2)-1)/(1 - (sqrt(2)-1)**2/tan(0.5*psi)**2)**0.5
      Y = (sqrt(2)-1)**2/((3 - 1/tan(0.5*psi)**2)**0.5)
      p = (2./15)*((3*(sqrt(2)-1)+4/sqrt(3))*sin(psi) - 6*(1-cos(psi))) \
        - 8./(5*pi)*(2*(sqrt(2)-1)*acos(X/tan(0.5*psi)) + 1./sqrt(3)*acos(Y/tan(0.5*psi)))*sin(psi) \
        + 8./(5*pi)*(2*acos((sqrt(2)+1)*X/sqrt(2)) + acos((sqrt(2)+1)*Y/sqrt(2)))*(1-cos(psi))
    else:
      p = 0.
    return p
    
  @staticmethod
  def misorientation_axis_from_delta(delta):
    '''Compute the misorientation axis from the misorientation matrix. 
    
    :param delta: The 3x3 misorientation matrix.
    :returns: the misorientation axis (normalised vector).
    '''
    n = np.array([delta[1,2] - delta[2,1], delta[2,0] - delta[0,2], delta[0,1] - delta[1,0]])
    n /= np.sqrt((delta[1,2] - delta[2,1])**2 + (delta[2,0] - delta[0,2])**2 + (delta[0,1] - delta[1,0])**2)
    return n

  def misorientation_axis(self, orientation):
    '''Compute the misorientation axis with another crystal orientation. 
    This vector is by definition common to both crystalline orientations.
    
    :param orientation: an instance of :py:class:`~pymicro.crystal.microstructure.Orientation` class.
    :returns: the misorientation axis (normalised vector).
    '''
    delta = np.dot(self.orientation_matrix(), orientation.orientation_matrix().T)
    return Orientation.misorientation_axis_from_delta(delta)
    
  @staticmethod
  def misorientation_angle_from_delta(delta):
    '''Compute the misorientation angle from the misorientation matrix.
    
    Compute the angle assocated with this misorientation matrix :math:`\\Delta g`. It is defined as :math:`\\omega = \\arccos(\\text{trace}(\\Delta g)/2-1)`.
    
    .. note::
    
      This does not account for the crystal symmetries. If you want to 
      find the disorientation between two orientations, use the 
      :py:meth:`~pymicro.crystal.microstructure.Orientation.disorientation` 
      method.
    
    :param delta: The 3x3 misorientation matrix.
    :returns float: the misorientation angle in radians.
    '''
    cw = 0.5*(delta.trace() - 1)
    omega = np.arccos(cw)
    return omega

  def disorientation(self, orientation, crystal_structure='cubic'):
    '''Compute the disorientation another crystal orientation.
    
    Considering all the possible crystal symmetries, the disorientation 
    is defined as the combination of the minimum misorientation angle 
    and the misorientation axis lying in the fundamental zone, which 
    can be used to bring the two lattices into coincidence.
    
    :param orientation: an instance of :py:class:`~pymicro.crystal.microstructure.Orientation` class desribing the other crystal orientation from which to compute the angle.
    :param str crystal_structure: a string describing the crystal structure, 'cubic' by default.
    :returns tuple: the misorientation angle in radians, the axis as a numpy vector (crystal coordinates), the axis as a numpy vector (sample coordinates).
    '''
    the_angle = np.pi
    from pymicro.crystal.lattice import Lattice
    symmetries = Lattice.symmetry(crystal_structure)
    (gA, gB) = (self.orientation_matrix(), orientation.orientation_matrix()) # nicknames
    for (g1, g2) in [(gA, gB), (gB, gA)]:
      for j in range(symmetries.shape[0]):
        sym_j = symmetries[j]
        oj = np.dot(sym_j, g1) # the crystal symmetry operator is left applied
        for i in range(symmetries.shape[0]):
          sym_i = symmetries[i]
          oi = np.dot(sym_i, g2)
          delta = np.dot(oi, oj.T)
          mis_angle = Orientation.misorientation_angle_from_delta(delta)
          if mis_angle < the_angle:
            # now compute the misorientation axis, should check if it lies in the fundamental zone
            mis_axis = Orientation.misorientation_axis_from_delta(delta)
            # here we have np.dot(oi.T, mis_axis) = np.dot(oj.T, mis_axis)
            #print(mis_axis, mis_angle*180/np.pi, np.dot(oj.T, mis_axis))
            the_angle = mis_angle
            the_axis = mis_axis
            the_axis_xyz = np.dot(oi.T, the_axis)
    return (the_angle, the_axis, the_axis_xyz)

  def phi1(self):
    '''Convenience methode to expose the first Euler angle.'''
    return self.euler[0]

  def Phi(self):
    '''Convenience methode to expose the second Euler angle.'''
    return self.euler[1]

  def phi2(self):
    '''Convenience methode to expose the third Euler angle.'''
    return self.euler[2]

  def compute_XG_angle(self, hkl, omega, verbose=False):
    '''Compute the angle between the scattering vector :math:`\mathbf{G_{l}}` 
    and :math:`\mathbf{-X}` the X-ray unit vector at a given angular position :math:`\\omega`.
    
    A given hkl plane defines the scattering vector :math:`\mathbf{G_{hkl}}` by 
    the miller indices in the reciprocal space. It is expressed in the 
    cartesian coordinate system by :math:`\mathbf{B}.\mathbf{G_{hkl}}` and in the 
    laboratory coordinate system accounting for the crystal orientation 
    by :math:`\mathbf{g}^{-1}.\mathbf{B}.\mathbf{G_{hkl}}`.
    
    The crystal is assumed to be placed on a rotation stage around the 
    laboratory vertical axis. The scattering vector can finally be 
    written as :math:`\mathbf{G_l}=\mathbf{\\Omega}.\mathbf{g}^{-1}.\mathbf{B}.\mathbf{G_{hkl}}`.
    The X-rays unit vector is :math:`\mathbf{X}=[1, 0, 0]`. So the computed angle 
    is :math:`\\alpha=acos(-\mathbf{X}.\mathbf{G_l}/||\mathbf{G_l}||`
    
    The Bragg condition is fulfilled when :math:`\\alpha=\pi/2-\\theta_{Bragg}`
    
    :param hkl: the hkl plane, an instance of :py:class:`~pymicro.crystal.lattice.HklPlane` 
    :param omega: the angle of rotation of the crystal around the laboratory vertical axis.
    :param bool verbose: activate verbose mode (False by default).
    :return float: the angle between :math:`-\mathbf{X}` and :math:`\mathbf{G_{l}}` in degrees.
    '''
    X = np.array([1., 0., 0.])
    Bt = self.orientation_matrix().transpose()
    Gc = hkl.scattering_vector()
    Gs = Bt.dot(Gc) # in the cartesian sample CS
    omegar = omega*np.pi/180
    R = np.array([[np.cos(omegar), -np.sin(omegar), 0], [np.sin(omegar), np.cos(omegar), 0], [0, 0, 1]])
    Gl = R.dot(Gs)
    alpha = np.arccos(np.dot(-X, Gl)/np.linalg.norm(Gl))*180/np.pi
    if verbose:
      print('scattering vector in the crystal CS', Gc)
      print('scattering vector in the sample CS', Gs)
      print('scattering vector in the laboratory CS (including Omega rotation)', Gl)
      print('angle (deg) between -X and G', alpha)
    return alpha

  def dct_omega_angles(self, hkl, lambda_keV, verbose=False):
    '''Compute the two omega angles which satisfy the Bragg condition.
    
    For a given crystal orientation sitting on a vertical rotation axis, 
    there is exactly two omega positions in :math:`[0, 2\pi]` for which 
    a particular hkl reflexion will fulfil Bragg's law.
    
    According to the Bragg's law, a crystallographic plane of a given 
    grain will be in diffracting condition if:
    
    .. math::
    
       \sin\\theta=-[\mathbf{\Omega}.\mathbf{g}^{-1}\mathbf{G_c}]_1
    
    with :math:`\mathbf{\Omega}` the matrix associated with the rotation 
    axis:

    .. math::
    
       \mathbf{\Omega}=\\begin{array}
                       \cos\omega & -\sin\omega & 0 \\\
                       \sin\omega & \cos\omega  & 0 \\\
                       0          & 0           & 1 \\\
                       \\end{array}

    This method solves the associated second order equation to return 
    the two corresponding omega angles.
    
    :param hkl: The given cristallographic plane :py:class:`~pymicro.crystal.lattice.HklPlane`
    :param float lambda_keV: The X-rays energy expressed in keV
    :param bool verbose: Verbose mode (False by default)
    :returns tuple: :math:`(\omega_1, \omega_2)` the two values of the \
    rotation angle around the vertical axis (in degrees).    
    '''
    (h, k, l) = hkl.miller_indices()
    (a, b, c) = hkl._lattice._lengths
    theta = hkl.bragg_angle(lambda_keV, verbose=verbose)
    lambda_nm = 1.2398 / lambda_keV

    gt = self.orientation_matrix() # our B (here called gt) corresponds to g^{-1} in Poulsen 2004
    A = h*gt[0,0] + k*gt[1,0] + l*gt[2,0]
    B = -h*gt[0,1] - k*gt[1,1] - l*gt[2,1]
    C = -2*a*np.sin(theta)**2 / lambda_nm # the minus sign comes from the main equation
    Delta = 4*(A**2 + B**2 - C**2)
    if Delta < 0:
      raise ValueError('Delta < 0')
    t1 = (B - 0.5 * np.sqrt(Delta)) / (A + C)
    t2 = (B + 0.5 * np.sqrt(Delta)) / (A + C)
    w1 = 2 * np.arctan(t1) * 180. / np.pi
    w2 = 2 * np.arctan(t2) * 180. / np.pi
    if w1 < 0: w1 += 360.
    if w2 < 0: w2 += 360.
    if verbose:
      print('the two omega values in degrees fulfilling the Bragg condition are (%.1f, %.1f)' % (w1, w2))
    return (w1, w2)
  
  def rotating_crystal(self, hkl, lambda_keV, omega_step=0.5, display=True, verbose=False):
    from pymicro.xray.xray_utils import lambda_keV_to_nm
    lambda_nm = lambda_keV_to_nm(lambda_keV)
    X = np.array([1., 0., 0.])/lambda_nm
    print 'magnitude of X', np.linalg.norm(X)
    Bt = self.orientation_matrix().transpose()
    (h, k, l) = hkl.miller_indices()
    theta = hkl.bragg_angle(lambda_keV)*180./np.pi
    print('bragg angle for %d%d%d reflection is %.1f' % (h, k, l, theta))
    Gc = hkl.scattering_vector()
    Gs = Bt.dot(Gc)
    alphas = []
    twothetas = []
    magnitude_K = []
    omegas = np.linspace(0.0, 360.0, num=360.0/omega_step, endpoint=False)
    for omega in omegas:
      print('\n** COMPUTING AT OMEGA=%03.1f deg' % omega)
      # prepare rotation matrix
      omegar = omega*np.pi/180
      R = np.array([[np.cos(omegar), -np.sin(omegar), 0], [np.sin(omegar), np.cos(omegar), 0], [0, 0, 1]])
      #R = R.dot(Rlt).dot(Rut) # with tilts
      Gl = R.dot(Gs)
      print('scattering vector in laboratory CS', Gl)
      n = R.dot(Bt.dot(hkl.normal()))
      print 'plane normal:', hkl.normal()
      print R
      print 'rotated plane normal:', n, ' with a norm of', np.linalg.norm(n)
      G = n/hkl.interplanar_spacing() # here G == N
      print 'G vector:', G, ' with a norm of', np.linalg.norm(G)
      K = X + G
      print 'X + G vector', K
      magnitude_K.append(np.linalg.norm(K))
      print 'magnitude of K', np.linalg.norm(K)
      alpha = np.arccos(np.dot(-X, G)/(np.linalg.norm(-X)*np.linalg.norm(G)))*180/np.pi
      print 'angle between -X and G', alpha
      alphas.append(alpha)
      twotheta = np.arccos(np.dot(K, X)/(np.linalg.norm(K)*np.linalg.norm(X)))*180/np.pi
      print 'angle (deg) between K and X', twotheta
      twothetas.append(twotheta)
    print 'min alpha angle is ', min(alphas)
    
    # compute omega_1 and omega_2 to verify graphically
    (w1, w2) = self.dct_omega_angles(hkl, lambda_keV, verbose=False)

    # gather the results in a single figure
    fig = plt.figure(figsize=(12,10))
    fig.add_subplot(311)
    plt.title('Looking for (%d%d%d) Bragg reflexions' % (h, k, l))
    plt.plot(omegas, alphas, 'k-')
    plt.xlim(0,360)
    plt.ylim(0,180)
    plt.xticks(np.arange(0,390, 30))
    # add bragg condition
    plt.axhline(90-theta, xmin=0, xmax=360, linewidth=2)
    plt.annotate('$\pi/2-\\theta_{Bragg}$', xycoords='data', xy=(360,90-theta), horizontalalignment='left', verticalalignment='center', fontsize=16)
    # add omega solutions
    plt.axvline(w1 + 180, ymin=0, ymax=180, linewidth=2, linestyle='dashed', color='gray')
    plt.axvline(w2 + 180, ymin=0, ymax=180, linewidth=2, linestyle='dashed', color='gray')
    plt.annotate('$\\omega_1$', xycoords='data', xy=(w1+180,0), horizontalalignment='center', verticalalignment='bottom', fontsize=16)
    plt.annotate('$\\omega_2$', xycoords='data', xy=(w2+180,0), horizontalalignment='center', verticalalignment='bottom', fontsize=16)
    plt.ylabel(r'Angle between $-X$ and $\mathbf{G}$')
    fig.add_subplot(312)
    plt.plot(omegas, twothetas, 'k-')
    plt.xlim(0,360)
    #plt.ylim(0,180)
    plt.xticks(np.arange(0,390, 30))
    plt.axhline(2*theta, xmin=0, xmax=360, linewidth=2)
    plt.annotate('$2\\theta_{Bragg}$', xycoords='data', xy=(360,2*theta), horizontalalignment='left', verticalalignment='center', fontsize=16)
    plt.axvline(w1 + 180, linewidth=2, linestyle='dashed', color='gray')
    plt.axvline(w2 + 180, linewidth=2, linestyle='dashed', color='gray')
    plt.ylabel('Angle between $X$ and $K$')
    fig.add_subplot(313)
    plt.plot(omegas, magnitude_K, 'k-')
    plt.xlim(0,360)
    plt.axhline(np.linalg.norm(X), xmin=0, xmax=360, linewidth=2)
    plt.annotate('$1/\\lambda$', xycoords='data', xy=(360,1/lambda_nm), horizontalalignment='left', verticalalignment='center', fontsize=16)
    plt.axvline(w1 + 180, linewidth=2, linestyle='dashed', color='gray')
    plt.axvline(w2 + 180, linewidth=2, linestyle='dashed', color='gray')
    plt.xlabel(r'Angle of rotation $\omega$')
    plt.ylabel(r'Magnitude of $X+G$ (nm$^{-1}$)')
    plt.subplots_adjust(top=0.925,bottom=0.05,left=0.1,right=0.9)
    if display:
      plt.show()
    else:
      plt.savefig('rotating_crystal_plot_%d%d%d.pdf' % (h, k, l))

  def topotomo_tilts(self, hkl, verbose=False):
    '''Compute the tilts for topotomography alignment.
    
    :param hkl: the hkl plane, an instance of :py:class:`~pymicro.crystal.lattice.HklPlane` 
    :param bool verbose: activate verbose mode (False by default).
    :returns tuple: (ut, lt) the two values of tilts to apply (in radians).
    '''
    Bt = self.orientation_matrix().transpose()
    Gc = hkl.scattering_vector()
    Gs = Bt.dot(Gc) # in the cartesian sample CS
    # find topotomo tilts
    ut = np.arctan(-Gs[0]/Gs[2])
    lt = np.arctan(Gs[1]/(-Gs[0]*np.sin(ut) + Gs[2]*np.cos(ut)))
    if verbose:
      print('up tilt (samry) should be %.3f' % (ut*180/np.pi))
      print('low tilt (samrx) should be %.3f' % (lt*180/np.pi))
    return (ut, lt)

  def to_xml(self, doc):
    '''
    Returns an XML representation of the Orientation instance.
    '''
    orientation = doc.createElement('Orientation')
    orientation_phi1 = doc.createElement('phi1')
    orientation_phi1_text = doc.createTextNode('%f' % self.phi1())
    orientation_phi1.appendChild(orientation_phi1_text)
    orientation.appendChild(orientation_phi1)
    orientation_Phi = doc.createElement('Phi')
    orientation_Phi_text = doc.createTextNode('%f' % self.Phi())
    orientation_Phi.appendChild(orientation_Phi_text)
    orientation.appendChild(orientation_Phi)
    orientation_phi2 = doc.createElement('phi2')
    orientation_phi2_text = doc.createTextNode('%f' % self.phi2())
    orientation_phi2.appendChild(orientation_phi2_text)
    orientation.appendChild(orientation_phi2)
    return orientation

  @staticmethod
  def from_xml(orientation_node):
    orientation_phi1 = orientation_node.childNodes[0]
    orientation_Phi = orientation_node.childNodes[1]
    orientation_phi2 = orientation_node.childNodes[2]
    phi1 = float(orientation_phi1.childNodes[0].nodeValue)
    Phi = float(orientation_Phi.childNodes[0].nodeValue)
    phi2 = float(orientation_phi2.childNodes[0].nodeValue)
    orientation = Orientation.from_euler(np.array([phi1, Phi, phi2]))
    return orientation

  @staticmethod
  def from_euler(euler):
    g = Orientation.Euler2OrientationMatrix(euler)
    o = Orientation(g)
    return o

  @staticmethod
  def from_rodrigues(rod):
    g = Orientation.Rodrigues2OrientationMatrix(rod)
    o = Orientation(g)
    return o
    
  #@staticmethod
  #def from_Zrot(axes):
    #g = Orientation.Zrot2OrientationMatrix()
    #o = Orientation(g)
    #return o
    
    
  @staticmethod
  def Euler2OrientationMatrix(euler):
    '''
    Compute the orientation matrix associated with the 3 Euler angles 
    (given in degrees).
    '''
    (rphi1, rPhi, rphi2) = np.radians(euler)
    c1 = np.cos(rphi1)
    s1 = np.sin(rphi1)
    c = np.cos(rPhi)
    s = np.sin(rPhi)
    c2 = np.cos(rphi2)
    s2 = np.sin(rphi2)

    # rotation matrix B
    b11 = c1*c2 - s1*s2*c
    b12 = s1*c2 + c1*s2*c
    b13 = s2*s
    b21 = -c1*s2 - s1*c2*c
    b22 = -s1*s2 + c1*c2*c
    b23 = c2*s
    b31 = s1*s
    b32 = -c1*s
    b33 = c
    B = np.array([[b11, b12, b13], [b21, b22, b23], [b31, b32, b33]])
    return B
  
  @staticmethod
  def Zrot2OrientationMatrix(x1=None,x2=None,x3=None):
    '''Compute the orientation matrix  from the rotated coordinates given in the 
       .inp file for Zebulon's computations
       Need at least two vectors to compute cross product
       
       Still need some tests to validate this function
    '''
    
    if (x1 is None and x2 is None):
      raise NameError('Need at least two vectors to compute the matrix')
    elif (x1==None and x3==None):
      raise NameError('Need at least two vectors to compute the matrix')
    elif (x3==None and x2==None):
      raise NameError('Need at least two vectors to compute the matrix')
    
    if x1 == None:
      x1 = np.cross(x2,x3)
    elif x2 == None:
      x2 = np.cross(x3,x1)
    elif x3 == None:
      x3 = np.cross(x1,x2)

    x1 = x1/np.linalg.norm(x1)
    x2 = x2/np.linalg.norm(x2)
    x3 = x3/np.linalg.norm(x3)
      
    B = np.array([x1, x2, x3]).transpose()
    return B
    
  @staticmethod
  def OrientationMatrix2EulerSF(g):
    '''
    Compute the Euler angles (in degrees) from the orientation matrix
    in a similar way as done in MandelCrystal.c
    '''
    tol = 0.1
    r = np.zeros(9, dtype=np.float64) # double precision here
    # Z-set order for tensor is 11 22 33 12 23 13 21 32 31
    r[0] = g[0, 0]
    r[1] = g[1, 1]
    r[2] = g[2, 2]
    r[3] = g[0, 1]
    r[4] = g[1, 2]
    r[5] = g[0, 2]
    r[6] = g[1, 0]
    r[7] = g[2, 1]
    r[8] = g[2, 0]
    phi = np.arccos(r[2])
    if (phi == 0.):
      phi2 = 0.
      phi1 = np.arcsin(r[6])
      if (abs(np.cos(phi1) - r[0]) > tol): phi1 = np.pi - phi1
    else:
      x2 = r[5]/np.sin(phi)
      x1 = r[8]/np.sin(phi);
      if (x1 > 1.): x1=1.
      if (x2 > 1.): x2=1.
      if (x1 < -1.): x1=-1.
      if (x2 < -1.): x2=-1.
      phi2 = np.arcsin(x2)
      phi1 = np.arcsin(x1)
      if (abs(np.cos(phi2) * np.sin(phi) - r[7]) > tol): phi2 = np.pi - phi2
      if (abs(np.cos(phi1) * np.sin(phi) + r[4]) > tol): phi1 = np.pi - phi1
    return np.degrees(np.array([phi1, phi, phi2]))
  
  @staticmethod
  def OrientationMatrix2Euler(g, eps=0.00000001):
    '''
    Compute the Euler angles (in degrees) from the orientation matrix.
    '''
    # compute rodrigues vector
    R = Orientation.OrientationMatrix2Rodrigues(g, eps=eps)
    if np.abs(g[2,2] - 1) < eps:
      phi1 = 0.5*np.arctan(g[0,1] / g[0,0])
      Phi = 0.0
      phi2 = phi1
    else:
      phi1 = np.arctan(-g[2,0] / g[2,1])
      Phi = np.arccos(g[2,2])
      phi2 = np.arctan(g[0,2] / g[1,2])
    # conditions determined by <NAME>, fev 2012
    if phi1 < 0:
      if phi2 > 0:
        if R[0] > 0:
          phi1 += 2*np.pi
        else:
          phi1 += np.pi
          phi2 += np.pi
      else: # phi2 < 0
        if R[0] > 0:
          phi1 += 2*np.pi
          phi2 += 2*np.pi
        else:
          phi1 += np.pi
          phi2 += np.pi
    else: # phi1 > 0
      if phi2 > 0:
        if R[0] > 0:
          pass # nothing special here
        else:
          phi1 += np.pi
          phi2 += np.pi
      else: # phi2 < 0
        if R[0] > 0:
          pass # nothing special here
        else:
          phi1 += np.pi
          phi2 += np.pi
    return np.degrees(np.array([phi1, Phi, phi2]))

  @staticmethod
  def OrientationMatrix2Rodrigues(g, eps=0.00000001):
    '''
    Compute the rodrigues vector from the orientation matrix.
    '''
    t = g.trace() + 1
    if np.abs(t) < eps:
      return np.zeros(3)
    else:
      r1 = (g[1,2] - g[2,1]) / t
      r2 = (g[2,0] - g[0,2]) / t
      r3 = (g[0,1] - g[1,0]) / t
    return np.array([r1, r2, r3])

  @staticmethod
  def Rodrigues2OrientationMatrix(rod, eps=0.000001):
    '''
    Compute the orientation matrix from the rodrigues vector.
    '''
    r = np.linalg.norm(rod)
    I = np.diagflat(np.ones(3))
    if r < eps:
      return I
    else:
      theta = 2*np.arctan(r)
      n = rod / r
      omega = np.array([[0.0, n[2], -n[1]], [-n[2], 0.0, n[0]], [n[1], -n[0], 0.0]])
      return I + np.sin(theta) * omega + (1 - np.cos(theta)) * omega.dot(omega)

  @staticmethod
  def Euler2Quaternion(euler):
    '''
    Compute the quaternion from the 3 euler angles (in degrees)
    '''
    (phi1, Phi, phi2) = np.radians(euler)
    q0 = np.cos(0.5 * (phi1 + phi2)) * np.cos(0.5*Phi)
    q1 = np.cos(0.5 * (phi1 - phi2)) * np.sin(0.5*Phi)
    q2 = np.sin(0.5 * (phi1 - phi2)) * np.sin(0.5*Phi)
    q3 = np.sin(0.5 * (phi1 + phi2)) * np.cos(0.5*Phi)
    return np.array([q0, q1, q2, q3])

  @staticmethod
  def Euler2Rodrigues(euler):
    '''
    Compute the rodrigues vector from the 3 euler angles (in degrees)
    '''
    (phi1, Phi, phi2) = np.radians(euler)
    a = 0.5 * (phi1 - phi2)
    b = 0.5 * (phi1 + phi2)
    r1 = np.tan(0.5*Phi) * np.cos(a) / np.cos(b)
    r2 = np.tan(0.5*Phi) * np.sin(a) / np.cos(b)
    r3 = np.tan(b)
    return np.array([r1, r2, r3])
    
  @staticmethod
  def read_euler_txt(txt_path, data_type='euler'):
    '''
    Read a set of grain orientations from a text file.
    
    The text file must be organised in columns, the three euler angles 
    (or the three rodrigues components) must correspond to the first 
    three columns (the other are ignored).
    
    :param str txt_path: path to the text file containing the orientations.
    :param str data_type: 'euler' (default) or 'rodrigues'
    :returns dict: a dictionary with the line number and the corresponding orientation.
    '''
    data = np.genfromtxt(txt_path)
    size = len(data)
    orientations = []
    for i in range(size):
      phi1 = data[i, 0]
      Phi = data[i, 1]
      phi2 = data[i, 2]
      angles = np.array([float(phi1), float(Phi), float(phi2)])
      if data_type == 'euler':
        orientations.append([i+1, Orientation.from_euler(angles)])
      elif data_type == 'rodrigues':
        orientations.append([i+1, Orientation.from_rodrigues(angles)])
    return dict(orientations)

  @staticmethod
  def read_euler_from_zset_inp(inp_path):
    '''Read a set of grain orientations from a z-set input file.
    
    In z-set input files, the orientation data may be specified 
    either using the rotation of two vector, euler angles or 
    rodrigues components directly. For instance the following 
    lines are extracted from a polycrystalline calculation file 
    using the rotation keyword:

    ::

     **elset elset1 *file au.mat *integration theta_method_a 1.0 1.e-9 150 *rotation x1 0.438886 -1.028805 0.197933 x3 1.038339 0.893172 1.003888
     **elset elset2 *file au.mat *integration theta_method_a 1.0 1.e-9 150 *rotation x1 0.178825 -0.716937 1.043300 x3 0.954345 0.879145 1.153101
     **elset elset3 *file au.mat *integration theta_method_a 1.0 1.e-9 150 *rotation x1 -0.540479 -0.827319 1.534062 x3 1.261700 1.284318 1.004174
     **elset elset4 *file au.mat *integration theta_method_a 1.0 1.e-9 150 *rotation x1 -0.941278 0.700996 0.034552 x3 1.000816 1.006824 0.885212
     **elset elset5 *file au.mat *integration theta_method_a 1.0 1.e-9 150 *rotation x1 -2.383786 0.479058 -0.488336 x3 0.899545 0.806075 0.984268

    :param str inp_path: the path to the ascii file to read.
    :returns dict: a dictionary of the orientations associated with the elset names.
    '''
    inp = open(inp_path)
    lines = inp.readlines()
    for i,line in enumerate(lines):  
      if line.lstrip().startswith('***material'):
        break
    euler_lines = []
    for j,line in enumerate(lines[i+1:]):
      # read until next *** block
      if line.lstrip().startswith('***'):
        break
      if (not line.lstrip().startswith('%') and line.find('**elset')>=0):
        euler_lines.append(line)
    euler = []
    for l in euler_lines:
      tokens = l.split()
      elset = tokens[tokens.index('**elset') + 1]
      irot = tokens.index('*rotation')
      if tokens[irot+1] == 'x1':
        x1 = np.empty(3, dtype=float)
        x1[0] = float(tokens[irot + 2])
        x1[1] = float(tokens[irot + 3])
        x1[2] = float(tokens[irot + 4])
        x3 = np.empty(3, dtype=float)
        x3[0] = float(tokens[irot + 6])
        x3[1] = float(tokens[irot + 7])
        x3[2] = float(tokens[irot + 8])
        euler.append([elset, Orientation.Zrot2OrientationMatrix(x1=x1, x3=x3)])
      else: # euler angles
        phi1 = tokens[irot + 1]
        Phi = tokens[irot + 2]
        phi2 = tokens[irot + 3]
        angles = np.array([float(phi1), float(Phi), float(phi2)])
        euler.append([elset, Orientation.from_euler(angles)])
    return dict(euler)
  
  def slip_system_orientation_tensor(self, s):
    '''Compute the orientation strain tensor m^s for this :py:class:`~pymicro.crystal.microstructure.Orientation` 
    and the given slip system.
    
    :param s: an instance of :py:class:`~pymicro.crystal.lattice.SlipSystem`

    .. math::

      M^s_{ij} = \left(l^s_i.n^s_j)
    '''
    Bt = self.orientation_matrix().transpose()
    plane = s.get_slip_plane()
    n_rot = np.dot(Bt, plane.normal())
    slip = s.get_slip_direction()
    l_rot = np.dot(Bt, slip.direction())
    return np.outer(l_rot, n_rot)

  def slip_system_orientation_strain_tensor(self, s):
    '''Compute the orientation strain tensor m^s for this :py:class:`~pymicro.crystal.microstructure.Orientation` 
    and the given slip system.
    
    :param s: an instance of :py:class:`~pymicro.crystal.lattice.SlipSystem`

    .. math::

      m^s_{ij} = \\frac{1}{2}\left(l^s_i.n^s_j + l^s_j.n^s_i)
    '''
    Bt = self.orientation_matrix().transpose()
    plane = s.get_slip_plane()
    n_rot = np.dot(Bt, plane.normal())
    slip = s.get_slip_direction()
    l_rot = np.dot(Bt, slip.direction())
    m = 0.5 * (np.outer(l_rot, n_rot) + np.outer(n_rot, l_rot))
    return m

  def slip_system_orientation_rotation_tensor(self, s):
    '''Compute the orientation rotation tensor q^s for this :py:class:`~pymicro.crystal.microstructure.Orientation` 
    and the given slip system.
    
    :param s: an instance of :py:class:`~pymicro.crystal.lattice.SlipSystem`

    .. math::

      q^s_{ij} = \\frac{1}{2}\left(l^s_i.n^s_j - l^s_j.n^s_i)
    '''
    Bt = self.orientation_matrix().transpose()
    plane = s.get_slip_plane()
    n_rot = np.dot(Bt, plane.normal())
    slip = s.get_slip_direction()
    l_rot = np.dot(Bt, slip.direction())
    q = 0.5 * (np.outer(l_rot, n_rot) - np.outer(n_rot, l_rot))
    return q

  def schmid_factor(self, slip_system, load_direction=[0., 0., 1]):
    '''Compute the Schmid factor for this crystal orientation and the 
    given slip system.

    :param slip_system: a slip system instance.
    :param load_direction: a unit vector describing the loading direction (default: vertical axis [0, 0, 1]).
    :returns float: a number between 0 ad 0.5.
    '''
    plane = slip_system.get_slip_plane()
    Bt = self.orientation_matrix().transpose()
    n_rot = np.dot(Bt, plane.normal()) # plane.normal() is a unit vector
    slip = slip_system.get_slip_direction().direction()
    slip_rot = np.dot(Bt, slip)
    #print n_rot, load_direction, slip_rot
    SF = np.abs(np.dot(n_rot, load_direction) * np.dot(slip_rot, load_direction))
    return SF

  def compute_all_schmid_factors(self, slip_systems, load_direction=[0., 0., 1], verbose=False):
    '''Compute all Schmid factors for this crystal orientation and the 
    given list of slip systems.

    :param slip_systems: a list of the slip system from which to compute the Schmid factor values.
    :param load_direction: a unit vector describing the loading direction (default: vertical axis [0, 0, 1]).
    :param bool verbose: activate verbose mode.
    :returns list: a list of the schmid factors.
    '''
    SF_list = []
    for ss in slip_systems:
      print ss
      sf = self.schmid_factor(ss, load_direction)
      if verbose:
        print 'Slip system: %s, Schmid factor is %.3f' % (ss, sf)
      SF_list.append(sf)
    return SF_list

class Grain:
  '''
  Class defining a crystallographic grain.

  A grain has its own crystallographic orientation.
  An optional id for the grain may be specified.
  The position field is the center of mass of the grain in world coordinates.
  The volume of the grain is expressed in pixel/voxel unit.
  '''

  def __init__(self, grain_id, grain_orientation):
    self.id = grain_id
    self.orientation = grain_orientation
    self.position = (0, 0, 0)
    self.volume = 0 # warning not implemented
    self.vtkmesh = None

  def __repr__(self):
    '''Provide a string representation of the class.'''
    s = '%s\n * id = %d\n' % (self.__class__.__name__, self.id)
    s += ' * %s\n' % (self.orientation)
    s += ' * position (%f, %f, %f)\n' % (self.position)
    s += ' * has vtk mesh ? %s\n' % (self.vtkmesh != None)
    return s
  
  def schmid_factor(self, slip_system, load_direction=[0., 0., 1]):
    '''Compute the Schmid factor of this grain for the given slip system.

    **Parameters**:
    
    *slip_system*: a slip system instance.
    
    *load_direction*: a unit vector describing the loading direction.
    
    **Returns**
    
    The Schmid factor of this grain for the given slip system.
    '''<commit_msg>$message<commit_after>'''
The microstructure module provide elementary classes to describe a 
crystallographic granular microstructure such as mostly present in 
metallic materials.

It contains several classes which are used to describe a microstructure 
composed of several grains, each one having its own crystallographic 
orientation:

 * :py:class:`~pymicro.crystal.microstructure.Microstructure`
 * :py:class:`~pymicro.crystal.microstructure.Grain`
 * :py:class:`~pymicro.crystal.microstructure.Orientation`
'''
import numpy as np
import os, vtk
from matplotlib import pyplot as plt, colors, cm
from xml.dom.minidom import Document, parse

class Orientation:
  '''Crystallographic orientation class.

  This follows the passive rotation definition which means that it brings 
  the sample coordinate system into coincidence with the crystal coordinate 
  system. Then one may express a vector :math:`V_c` in the crystal coordinate system 
  from the vector in the sample coordinate system :math:`V_s` by:
  
  .. math::

    V_c = B.V_s
  
  and inversely (because :math:`B^{-1}=B^T`):
  
  .. math::

    V_s = B^T.V_c
  '''

  def __init__(self, matrix):
    '''Initialization from the 9 components of the orientation matrix.'''
    g = np.array(matrix, dtype=np.float64).reshape((3, 3))
    self._matrix = g
    self.euler = Orientation.OrientationMatrix2Euler(g)
    self.rod = Orientation.OrientationMatrix2Rodrigues(g)

  def orientation_matrix(self):
    '''Returns the orientation matrix in the form of a 3x3 numpy array.'''
    return self._matrix

  def __repr__(self):
    '''Provide a string representation of the class.'''
    s = 'Crystal Orientation'
    s += '\norientation matrix = %s' % self._matrix.view()
    s += '\nEuler angles (degrees) = (%8.3f,%8.3f,%8.3f)' % (self.phi1(), self.Phi(), self.phi2())
    s += '\nRodrigues vector = %s' % self.OrientationMatrix2Rodrigues(self._matrix)
    return s

  @staticmethod
  def cube():
    '''Create the particular crystal orienation called Cube and which 
    corresponds to euler angle (0,0,0).'''
    return Orientation.from_euler((0., 0., 0.))

  @staticmethod
  def brass():
    '''Create the particular crystal orienation called Brass and which 
    corresponds to euler angle (35.264,45,0).'''
    return Orientation.from_euler((35.264, 45., 0.))

  @staticmethod
  def copper():
    '''Create the particular crystal orienation called Copper and which 
    corresponds to euler angle (90,35.264,45).'''
    return Orientation.from_euler((90., 35.264, 45.))

  @staticmethod
  def s3():
    '''Create the particular crystal orienation called S3 and which 
    corresponds to euler angle (59,37,63).'''
    return Orientation.from_euler((58.980, 36.699, 63.435))

  @staticmethod
  def goss():
    '''Create the particular crystal orienation called Goss and which 
    corresponds to euler angle (0,45,0).'''
    return Orientation.from_euler((0., 45., 0.))

  @staticmethod
  def shear():
    '''Create the particular crystal orienation called shear and which 
    corresponds to euler angle (45,0,0).'''
    return Orientation.from_euler((45., 0., 0.))

  @staticmethod
  def misorientation_MacKenzie(psi):
    '''Return the fraction of the misorientations corresponding to the 
    given :math:`\\psi` angle in the refenrece solution derived By MacKenzie in 
    his 1958 paper.
    
    :param psi: the misorientation angle in radians.
    :returns: the value in the cummulative distribution corresponding to psi.
    '''
    from math import sqrt, sin, cos, tan, pi, acos
    psidg = 180*psi/pi
    if psidg >= 0 and psidg <= 45:
      p = 2./15*(1 - cos(psi))
    elif psidg > 45 and psidg <= 60:
      p = 2./15*(3*(sqrt(2)-1)*sin(psi) - 2*(1 - cos(psi)))
    elif psidg > 60 and psidg <= 60.72:
      p = 2./15*((3*(sqrt(2)-1) + 4./sqrt(3))*sin(psi) - 6.*(1 - cos(psi)))
    elif psidg > 60.72 and psidg <= 62.8:
      X = (sqrt(2)-1)/(1 - (sqrt(2)-1)**2/tan(0.5*psi)**2)**0.5
      Y = (sqrt(2)-1)**2/((3 - 1/tan(0.5*psi)**2)**0.5)
      p = (2./15)*((3*(sqrt(2)-1)+4/sqrt(3))*sin(psi) - 6*(1-cos(psi))) \
        - 8./(5*pi)*(2*(sqrt(2)-1)*acos(X/tan(0.5*psi)) + 1./sqrt(3)*acos(Y/tan(0.5*psi)))*sin(psi) \
        + 8./(5*pi)*(2*acos((sqrt(2)+1)*X/sqrt(2)) + acos((sqrt(2)+1)*Y/sqrt(2)))*(1-cos(psi))
    else:
      p = 0.
    return p
    
  @staticmethod
  def misorientation_axis_from_delta(delta):
    '''Compute the misorientation axis from the misorientation matrix. 
    
    :param delta: The 3x3 misorientation matrix.
    :returns: the misorientation axis (normalised vector).
    '''
    n = np.array([delta[1,2] - delta[2,1], delta[2,0] - delta[0,2], delta[0,1] - delta[1,0]])
    n /= np.sqrt((delta[1,2] - delta[2,1])**2 + (delta[2,0] - delta[0,2])**2 + (delta[0,1] - delta[1,0])**2)
    return n

  def misorientation_axis(self, orientation):
    '''Compute the misorientation axis with another crystal orientation. 
    This vector is by definition common to both crystalline orientations.
    
    :param orientation: an instance of :py:class:`~pymicro.crystal.microstructure.Orientation` class.
    :returns: the misorientation axis (normalised vector).
    '''
    delta = np.dot(self.orientation_matrix(), orientation.orientation_matrix().T)
    return Orientation.misorientation_axis_from_delta(delta)
    
  @staticmethod
  def misorientation_angle_from_delta(delta):
    '''Compute the misorientation angle from the misorientation matrix.
    
    Compute the angle assocated with this misorientation matrix :math:`\\Delta g`. It is defined as :math:`\\omega = \\arccos(\\text{trace}(\\Delta g)/2-1)`.
    
    .. note::
    
      This does not account for the crystal symmetries. If you want to 
      find the disorientation between two orientations, use the 
      :py:meth:`~pymicro.crystal.microstructure.Orientation.disorientation` 
      method.
    
    :param delta: The 3x3 misorientation matrix.
    :returns float: the misorientation angle in radians.
    '''
    cw = 0.5*(delta.trace() - 1)
    omega = np.arccos(cw)
    return omega

  def disorientation(self, orientation, crystal_structure='cubic'):
    '''Compute the disorientation another crystal orientation.
    
    Considering all the possible crystal symmetries, the disorientation 
    is defined as the combination of the minimum misorientation angle 
    and the misorientation axis lying in the fundamental zone, which 
    can be used to bring the two lattices into coincidence.
    
    :param orientation: an instance of :py:class:`~pymicro.crystal.microstructure.Orientation` class desribing the other crystal orientation from which to compute the angle.
    :param str crystal_structure: a string describing the crystal structure, 'cubic' by default.
    :returns tuple: the misorientation angle in radians, the axis as a numpy vector (crystal coordinates), the axis as a numpy vector (sample coordinates).
    '''
    the_angle = np.pi
    from pymicro.crystal.lattice import Lattice
    symmetries = Lattice.symmetry(crystal_structure)
    (gA, gB) = (self.orientation_matrix(), orientation.orientation_matrix()) # nicknames
    for (g1, g2) in [(gA, gB), (gB, gA)]:
      for j in range(symmetries.shape[0]):
        sym_j = symmetries[j]
        oj = np.dot(sym_j, g1) # the crystal symmetry operator is left applied
        for i in range(symmetries.shape[0]):
          sym_i = symmetries[i]
          oi = np.dot(sym_i, g2)
          delta = np.dot(oi, oj.T)
          mis_angle = Orientation.misorientation_angle_from_delta(delta)
          if mis_angle < the_angle:
            # now compute the misorientation axis, should check if it lies in the fundamental zone
            mis_axis = Orientation.misorientation_axis_from_delta(delta)
            # here we have np.dot(oi.T, mis_axis) = np.dot(oj.T, mis_axis)
            #print(mis_axis, mis_angle*180/np.pi, np.dot(oj.T, mis_axis))
            the_angle = mis_angle
            the_axis = mis_axis
            the_axis_xyz = np.dot(oi.T, the_axis)
    return (the_angle, the_axis, the_axis_xyz)

  def phi1(self):
    '''Convenience methode to expose the first Euler angle.'''
    return self.euler[0]

  def Phi(self):
    '''Convenience methode to expose the second Euler angle.'''
    return self.euler[1]

  def phi2(self):
    '''Convenience methode to expose the third Euler angle.'''
    return self.euler[2]

  def compute_XG_angle(self, hkl, omega, verbose=False):
    '''Compute the angle between the scattering vector :math:`\mathbf{G_{l}}` 
    and :math:`\mathbf{-X}` the X-ray unit vector at a given angular position :math:`\\omega`.
    
    A given hkl plane defines the scattering vector :math:`\mathbf{G_{hkl}}` by 
    the miller indices in the reciprocal space. It is expressed in the 
    cartesian coordinate system by :math:`\mathbf{B}.\mathbf{G_{hkl}}` and in the 
    laboratory coordinate system accounting for the crystal orientation 
    by :math:`\mathbf{g}^{-1}.\mathbf{B}.\mathbf{G_{hkl}}`.
    
    The crystal is assumed to be placed on a rotation stage around the 
    laboratory vertical axis. The scattering vector can finally be 
    written as :math:`\mathbf{G_l}=\mathbf{\\Omega}.\mathbf{g}^{-1}.\mathbf{B}.\mathbf{G_{hkl}}`.
    The X-rays unit vector is :math:`\mathbf{X}=[1, 0, 0]`. So the computed angle 
    is :math:`\\alpha=acos(-\mathbf{X}.\mathbf{G_l}/||\mathbf{G_l}||`
    
    The Bragg condition is fulfilled when :math:`\\alpha=\pi/2-\\theta_{Bragg}`
    
    :param hkl: the hkl plane, an instance of :py:class:`~pymicro.crystal.lattice.HklPlane` 
    :param omega: the angle of rotation of the crystal around the laboratory vertical axis.
    :param bool verbose: activate verbose mode (False by default).
    :return float: the angle between :math:`-\mathbf{X}` and :math:`\mathbf{G_{l}}` in degrees.
    '''
    X = np.array([1., 0., 0.])
    Bt = self.orientation_matrix().transpose()
    Gc = hkl.scattering_vector()
    Gs = Bt.dot(Gc) # in the cartesian sample CS
    omegar = omega*np.pi/180
    R = np.array([[np.cos(omegar), -np.sin(omegar), 0], [np.sin(omegar), np.cos(omegar), 0], [0, 0, 1]])
    Gl = R.dot(Gs)
    alpha = np.arccos(np.dot(-X, Gl)/np.linalg.norm(Gl))*180/np.pi
    if verbose:
      print('scattering vector in the crystal CS', Gc)
      print('scattering vector in the sample CS', Gs)
      print('scattering vector in the laboratory CS (including Omega rotation)', Gl)
      print('angle (deg) between -X and G', alpha)
    return alpha

  def dct_omega_angles(self, hkl, lambda_keV, verbose=False):
    '''Compute the two omega angles which satisfy the Bragg condition.
    
    For a given crystal orientation sitting on a vertical rotation axis, 
    there is exactly two omega positions in :math:`[0, 2\pi]` for which 
    a particular hkl reflexion will fulfil Bragg's law.
    
    According to the Bragg's law, a crystallographic plane of a given 
    grain will be in diffracting condition if:
    
    .. math::
    
       \sin\\theta=-[\mathbf{\Omega}.\mathbf{g}^{-1}\mathbf{G_c}]_1
    
    with :math:`\mathbf{\Omega}` the matrix associated with the rotation 
    axis:

    .. math::
    
       \mathbf{\Omega}=\\begin{array}
                       \cos\omega & -\sin\omega & 0 \\\
                       \sin\omega & \cos\omega  & 0 \\\
                       0          & 0           & 1 \\\
                       \\end{array}

    This method solves the associated second order equation to return 
    the two corresponding omega angles.
    
    :param hkl: The given cristallographic plane :py:class:`~pymicro.crystal.lattice.HklPlane`
    :param float lambda_keV: The X-rays energy expressed in keV
    :param bool verbose: Verbose mode (False by default)
    :returns tuple: :math:`(\omega_1, \omega_2)` the two values of the \
    rotation angle around the vertical axis (in degrees).    
    '''
    (h, k, l) = hkl.miller_indices()
    (a, b, c) = hkl._lattice._lengths
    theta = hkl.bragg_angle(lambda_keV, verbose=verbose)
    lambda_nm = 1.2398 / lambda_keV

    gt = self.orientation_matrix() # our B (here called gt) corresponds to g^{-1} in Poulsen 2004
    A = h*gt[0,0] + k*gt[1,0] + l*gt[2,0]
    B = -h*gt[0,1] - k*gt[1,1] - l*gt[2,1]
    C = -2*a*np.sin(theta)**2 / lambda_nm # the minus sign comes from the main equation
    Delta = 4*(A**2 + B**2 - C**2)
    if Delta < 0:
      raise ValueError('Delta < 0')
    t1 = (B - 0.5 * np.sqrt(Delta)) / (A + C)
    t2 = (B + 0.5 * np.sqrt(Delta)) / (A + C)
    w1 = 2 * np.arctan(t1) * 180. / np.pi
    w2 = 2 * np.arctan(t2) * 180. / np.pi
    if w1 < 0: w1 += 360.
    if w2 < 0: w2 += 360.
    if verbose:
      print('the two omega values in degrees fulfilling the Bragg condition are (%.1f, %.1f)' % (w1, w2))
    return (w1, w2)
  
  def rotating_crystal(self, hkl, lambda_keV, omega_step=0.5, display=True, verbose=False):
    from pymicro.xray.xray_utils import lambda_keV_to_nm
    lambda_nm = lambda_keV_to_nm(lambda_keV)
    X = np.array([1., 0., 0.])/lambda_nm
    print 'magnitude of X', np.linalg.norm(X)
    Bt = self.orientation_matrix().transpose()
    (h, k, l) = hkl.miller_indices()
    theta = hkl.bragg_angle(lambda_keV)*180./np.pi
    print('bragg angle for %d%d%d reflection is %.1f' % (h, k, l, theta))
    Gc = hkl.scattering_vector()
    Gs = Bt.dot(Gc)
    alphas = []
    twothetas = []
    magnitude_K = []
    omegas = np.linspace(0.0, 360.0, num=360.0/omega_step, endpoint=False)
    for omega in omegas:
      print('\n** COMPUTING AT OMEGA=%03.1f deg' % omega)
      # prepare rotation matrix
      omegar = omega*np.pi/180
      R = np.array([[np.cos(omegar), -np.sin(omegar), 0], [np.sin(omegar), np.cos(omegar), 0], [0, 0, 1]])
      #R = R.dot(Rlt).dot(Rut) # with tilts
      Gl = R.dot(Gs)
      print('scattering vector in laboratory CS', Gl)
      n = R.dot(Bt.dot(hkl.normal()))
      print 'plane normal:', hkl.normal()
      print R
      print 'rotated plane normal:', n, ' with a norm of', np.linalg.norm(n)
      G = n/hkl.interplanar_spacing() # here G == N
      print 'G vector:', G, ' with a norm of', np.linalg.norm(G)
      K = X + G
      print 'X + G vector', K
      magnitude_K.append(np.linalg.norm(K))
      print 'magnitude of K', np.linalg.norm(K)
      alpha = np.arccos(np.dot(-X, G)/(np.linalg.norm(-X)*np.linalg.norm(G)))*180/np.pi
      print 'angle between -X and G', alpha
      alphas.append(alpha)
      twotheta = np.arccos(np.dot(K, X)/(np.linalg.norm(K)*np.linalg.norm(X)))*180/np.pi
      print 'angle (deg) between K and X', twotheta
      twothetas.append(twotheta)
    print 'min alpha angle is ', min(alphas)
    
    # compute omega_1 and omega_2 to verify graphically
    (w1, w2) = self.dct_omega_angles(hkl, lambda_keV, verbose=False)

    # gather the results in a single figure
    fig = plt.figure(figsize=(12,10))
    fig.add_subplot(311)
    plt.title('Looking for (%d%d%d) Bragg reflexions' % (h, k, l))
    plt.plot(omegas, alphas, 'k-')
    plt.xlim(0,360)
    plt.ylim(0,180)
    plt.xticks(np.arange(0,390, 30))
    # add bragg condition
    plt.axhline(90-theta, xmin=0, xmax=360, linewidth=2)
    plt.annotate('$\pi/2-\\theta_{Bragg}$', xycoords='data', xy=(360,90-theta), horizontalalignment='left', verticalalignment='center', fontsize=16)
    # add omega solutions
    plt.axvline(w1 + 180, ymin=0, ymax=180, linewidth=2, linestyle='dashed', color='gray')
    plt.axvline(w2 + 180, ymin=0, ymax=180, linewidth=2, linestyle='dashed', color='gray')
    plt.annotate('$\\omega_1$', xycoords='data', xy=(w1+180,0), horizontalalignment='center', verticalalignment='bottom', fontsize=16)
    plt.annotate('$\\omega_2$', xycoords='data', xy=(w2+180,0), horizontalalignment='center', verticalalignment='bottom', fontsize=16)
    plt.ylabel(r'Angle between $-X$ and $\mathbf{G}$')
    fig.add_subplot(312)
    plt.plot(omegas, twothetas, 'k-')
    plt.xlim(0,360)
    #plt.ylim(0,180)
    plt.xticks(np.arange(0,390, 30))
    plt.axhline(2*theta, xmin=0, xmax=360, linewidth=2)
    plt.annotate('$2\\theta_{Bragg}$', xycoords='data', xy=(360,2*theta), horizontalalignment='left', verticalalignment='center', fontsize=16)
    plt.axvline(w1 + 180, linewidth=2, linestyle='dashed', color='gray')
    plt.axvline(w2 + 180, linewidth=2, linestyle='dashed', color='gray')
    plt.ylabel('Angle between $X$ and $K$')
    fig.add_subplot(313)
    plt.plot(omegas, magnitude_K, 'k-')
    plt.xlim(0,360)
    plt.axhline(np.linalg.norm(X), xmin=0, xmax=360, linewidth=2)
    plt.annotate('$1/\\lambda$', xycoords='data', xy=(360,1/lambda_nm), horizontalalignment='left', verticalalignment='center', fontsize=16)
    plt.axvline(w1 + 180, linewidth=2, linestyle='dashed', color='gray')
    plt.axvline(w2 + 180, linewidth=2, linestyle='dashed', color='gray')
    plt.xlabel(r'Angle of rotation $\omega$')
    plt.ylabel(r'Magnitude of $X+G$ (nm$^{-1}$)')
    plt.subplots_adjust(top=0.925,bottom=0.05,left=0.1,right=0.9)
    if display:
      plt.show()
    else:
      plt.savefig('rotating_crystal_plot_%d%d%d.pdf' % (h, k, l))

  def topotomo_tilts(self, hkl, verbose=False):
    '''Compute the tilts for topotomography alignment.
    
    :param hkl: the hkl plane, an instance of :py:class:`~pymicro.crystal.lattice.HklPlane` 
    :param bool verbose: activate verbose mode (False by default).
    :returns tuple: (ut, lt) the two values of tilts to apply (in radians).
    '''
    Bt = self.orientation_matrix().transpose()
    Gc = hkl.scattering_vector()
    Gs = Bt.dot(Gc) # in the cartesian sample CS
    # find topotomo tilts
    ut = np.arctan(-Gs[0]/Gs[2])
    lt = np.arctan(Gs[1]/(-Gs[0]*np.sin(ut) + Gs[2]*np.cos(ut)))
    if verbose:
      print('up tilt (samry) should be %.3f' % (ut*180/np.pi))
      print('low tilt (samrx) should be %.3f' % (lt*180/np.pi))
    return (ut, lt)

  def to_xml(self, doc):
    '''
    Returns an XML representation of the Orientation instance.
    '''
    orientation = doc.createElement('Orientation')
    orientation_phi1 = doc.createElement('phi1')
    orientation_phi1_text = doc.createTextNode('%f' % self.phi1())
    orientation_phi1.appendChild(orientation_phi1_text)
    orientation.appendChild(orientation_phi1)
    orientation_Phi = doc.createElement('Phi')
    orientation_Phi_text = doc.createTextNode('%f' % self.Phi())
    orientation_Phi.appendChild(orientation_Phi_text)
    orientation.appendChild(orientation_Phi)
    orientation_phi2 = doc.createElement('phi2')
    orientation_phi2_text = doc.createTextNode('%f' % self.phi2())
    orientation_phi2.appendChild(orientation_phi2_text)
    orientation.appendChild(orientation_phi2)
    return orientation

  @staticmethod
  def from_xml(orientation_node):
    orientation_phi1 = orientation_node.childNodes[0]
    orientation_Phi = orientation_node.childNodes[1]
    orientation_phi2 = orientation_node.childNodes[2]
    phi1 = float(orientation_phi1.childNodes[0].nodeValue)
    Phi = float(orientation_Phi.childNodes[0].nodeValue)
    phi2 = float(orientation_phi2.childNodes[0].nodeValue)
    orientation = Orientation.from_euler(np.array([phi1, Phi, phi2]))
    return orientation

  @staticmethod
  def from_euler(euler):
    g = Orientation.Euler2OrientationMatrix(euler)
    o = Orientation(g)
    return o

  @staticmethod
  def from_rodrigues(rod):
    g = Orientation.Rodrigues2OrientationMatrix(rod)
    o = Orientation(g)
    return o
    
  #@staticmethod
  #def from_Zrot(axes):
    #g = Orientation.Zrot2OrientationMatrix()
    #o = Orientation(g)
    #return o
    
    
  @staticmethod
  def Euler2OrientationMatrix(euler):
    '''
    Compute the orientation matrix associated with the 3 Euler angles 
    (given in degrees).
    '''
    (rphi1, rPhi, rphi2) = np.radians(euler)
    c1 = np.cos(rphi1)
    s1 = np.sin(rphi1)
    c = np.cos(rPhi)
    s = np.sin(rPhi)
    c2 = np.cos(rphi2)
    s2 = np.sin(rphi2)

    # rotation matrix B
    b11 = c1*c2 - s1*s2*c
    b12 = s1*c2 + c1*s2*c
    b13 = s2*s
    b21 = -c1*s2 - s1*c2*c
    b22 = -s1*s2 + c1*c2*c
    b23 = c2*s
    b31 = s1*s
    b32 = -c1*s
    b33 = c
    B = np.array([[b11, b12, b13], [b21, b22, b23], [b31, b32, b33]])
    return B
  
  @staticmethod
  def Zrot2OrientationMatrix(x1=None,x2=None,x3=None):
    '''Compute the orientation matrix  from the rotated coordinates given in the 
       .inp file for Zebulon's computations
       Need at least two vectors to compute cross product
       
       Still need some tests to validate this function
    '''
    
    if (x1 is None and x2 is None):
      raise NameError('Need at least two vectors to compute the matrix')
    elif (x1==None and x3==None):
      raise NameError('Need at least two vectors to compute the matrix')
    elif (x3==None and x2==None):
      raise NameError('Need at least two vectors to compute the matrix')
    
    if x1 == None:
      x1 = np.cross(x2,x3)
    elif x2 == None:
      x2 = np.cross(x3,x1)
    elif x3 == None:
      x3 = np.cross(x1,x2)

    x1 = x1/np.linalg.norm(x1)
    x2 = x2/np.linalg.norm(x2)
    x3 = x3/np.linalg.norm(x3)
      
    B = np.array([x1, x2, x3]).transpose()
    return B
    
  @staticmethod
  def OrientationMatrix2EulerSF(g):
    '''
    Compute the Euler angles (in degrees) from the orientation matrix
    in a similar way as done in MandelCrystal.c
    '''
    tol = 0.1
    r = np.zeros(9, dtype=np.float64) # double precision here
    # Z-set order for tensor is 11 22 33 12 23 13 21 32 31
    r[0] = g[0, 0]
    r[1] = g[1, 1]
    r[2] = g[2, 2]
    r[3] = g[0, 1]
    r[4] = g[1, 2]
    r[5] = g[0, 2]
    r[6] = g[1, 0]
    r[7] = g[2, 1]
    r[8] = g[2, 0]
    phi = np.arccos(r[2])
    if (phi == 0.):
      phi2 = 0.
      phi1 = np.arcsin(r[6])
      if (abs(np.cos(phi1) - r[0]) > tol): phi1 = np.pi - phi1
    else:
      x2 = r[5]/np.sin(phi)
      x1 = r[8]/np.sin(phi);
      if (x1 > 1.): x1=1.
      if (x2 > 1.): x2=1.
      if (x1 < -1.): x1=-1.
      if (x2 < -1.): x2=-1.
      phi2 = np.arcsin(x2)
      phi1 = np.arcsin(x1)
      if (abs(np.cos(phi2) * np.sin(phi) - r[7]) > tol): phi2 = np.pi - phi2
      if (abs(np.cos(phi1) * np.sin(phi) + r[4]) > tol): phi1 = np.pi - phi1
    return np.degrees(np.array([phi1, phi, phi2]))
  
  @staticmethod
  def OrientationMatrix2Euler(g, eps=0.00000001):
    '''
    Compute the Euler angles (in degrees) from the orientation matrix.
    '''
    # compute rodrigues vector
    R = Orientation.OrientationMatrix2Rodrigues(g, eps=eps)
    if np.abs(g[2,2] - 1) < eps:
      phi1 = 0.5*np.arctan(g[0,1] / g[0,0])
      Phi = 0.0
      phi2 = phi1
    else:
      phi1 = np.arctan(-g[2,0] / g[2,1])
      Phi = np.arccos(g[2,2])
      phi2 = np.arctan(g[0,2] / g[1,2])
    # conditions determined by <NAME>, fev 2012
    if phi1 < 0:
      if phi2 > 0:
        if R[0] > 0:
          phi1 += 2*np.pi
        else:
          phi1 += np.pi
          phi2 += np.pi
      else: # phi2 < 0
        if R[0] > 0:
          phi1 += 2*np.pi
          phi2 += 2*np.pi
        else:
          phi1 += np.pi
          phi2 += np.pi
    else: # phi1 > 0
      if phi2 > 0:
        if R[0] > 0:
          pass # nothing special here
        else:
          phi1 += np.pi
          phi2 += np.pi
      else: # phi2 < 0
        if R[0] > 0:
          pass # nothing special here
        else:
          phi1 += np.pi
          phi2 += np.pi
    return np.degrees(np.array([phi1, Phi, phi2]))

  @staticmethod
  def OrientationMatrix2Rodrigues(g, eps=0.00000001):
    '''
    Compute the rodrigues vector from the orientation matrix.
    '''
    t = g.trace() + 1
    if np.abs(t) < eps:
      return np.zeros(3)
    else:
      r1 = (g[1,2] - g[2,1]) / t
      r2 = (g[2,0] - g[0,2]) / t
      r3 = (g[0,1] - g[1,0]) / t
    return np.array([r1, r2, r3])

  @staticmethod
  def Rodrigues2OrientationMatrix(rod, eps=0.000001):
    '''
    Compute the orientation matrix from the rodrigues vector.
    '''
    r = np.linalg.norm(rod)
    I = np.diagflat(np.ones(3))
    if r < eps:
      return I
    else:
      theta = 2*np.arctan(r)
      n = rod / r
      omega = np.array([[0.0, n[2], -n[1]], [-n[2], 0.0, n[0]], [n[1], -n[0], 0.0]])
      return I + np.sin(theta) * omega + (1 - np.cos(theta)) * omega.dot(omega)

  @staticmethod
  def Euler2Quaternion(euler):
    '''
    Compute the quaternion from the 3 euler angles (in degrees)
    '''
    (phi1, Phi, phi2) = np.radians(euler)
    q0 = np.cos(0.5 * (phi1 + phi2)) * np.cos(0.5*Phi)
    q1 = np.cos(0.5 * (phi1 - phi2)) * np.sin(0.5*Phi)
    q2 = np.sin(0.5 * (phi1 - phi2)) * np.sin(0.5*Phi)
    q3 = np.sin(0.5 * (phi1 + phi2)) * np.cos(0.5*Phi)
    return np.array([q0, q1, q2, q3])

  @staticmethod
  def Euler2Rodrigues(euler):
    '''
    Compute the rodrigues vector from the 3 euler angles (in degrees)
    '''
    (phi1, Phi, phi2) = np.radians(euler)
    a = 0.5 * (phi1 - phi2)
    b = 0.5 * (phi1 + phi2)
    r1 = np.tan(0.5*Phi) * np.cos(a) / np.cos(b)
    r2 = np.tan(0.5*Phi) * np.sin(a) / np.cos(b)
    r3 = np.tan(b)
    return np.array([r1, r2, r3])
    
  @staticmethod
  def read_euler_txt(txt_path, data_type='euler'):
    '''
    Read a set of grain orientations from a text file.
    
    The text file must be organised in columns, the three euler angles 
    (or the three rodrigues components) must correspond to the first 
    three columns (the other are ignored).
    
    :param str txt_path: path to the text file containing the orientations.
    :param str data_type: 'euler' (default) or 'rodrigues'
    :returns dict: a dictionary with the line number and the corresponding orientation.
    '''
    data = np.genfromtxt(txt_path)
    size = len(data)
    orientations = []
    for i in range(size):
      phi1 = data[i, 0]
      Phi = data[i, 1]
      phi2 = data[i, 2]
      angles = np.array([float(phi1), float(Phi), float(phi2)])
      if data_type == 'euler':
        orientations.append([i+1, Orientation.from_euler(angles)])
      elif data_type == 'rodrigues':
        orientations.append([i+1, Orientation.from_rodrigues(angles)])
    return dict(orientations)

  @staticmethod
  def read_euler_from_zset_inp(inp_path):
    '''Read a set of grain orientations from a z-set input file.
    
    In z-set input files, the orientation data may be specified 
    either using the rotation of two vector, euler angles or 
    rodrigues components directly. For instance the following 
    lines are extracted from a polycrystalline calculation file 
    using the rotation keyword:

    ::

     **elset elset1 *file au.mat *integration theta_method_a 1.0 1.e-9 150 *rotation x1 0.438886 -1.028805 0.197933 x3 1.038339 0.893172 1.003888
     **elset elset2 *file au.mat *integration theta_method_a 1.0 1.e-9 150 *rotation x1 0.178825 -0.716937 1.043300 x3 0.954345 0.879145 1.153101
     **elset elset3 *file au.mat *integration theta_method_a 1.0 1.e-9 150 *rotation x1 -0.540479 -0.827319 1.534062 x3 1.261700 1.284318 1.004174
     **elset elset4 *file au.mat *integration theta_method_a 1.0 1.e-9 150 *rotation x1 -0.941278 0.700996 0.034552 x3 1.000816 1.006824 0.885212
     **elset elset5 *file au.mat *integration theta_method_a 1.0 1.e-9 150 *rotation x1 -2.383786 0.479058 -0.488336 x3 0.899545 0.806075 0.984268

    :param str inp_path: the path to the ascii file to read.
    :returns dict: a dictionary of the orientations associated with the elset names.
    '''
    inp = open(inp_path)
    lines = inp.readlines()
    for i,line in enumerate(lines):  
      if line.lstrip().startswith('***material'):
        break
    euler_lines = []
    for j,line in enumerate(lines[i+1:]):
      # read until next *** block
      if line.lstrip().startswith('***'):
        break
      if (not line.lstrip().startswith('%') and line.find('**elset')>=0):
        euler_lines.append(line)
    euler = []
    for l in euler_lines:
      tokens = l.split()
      elset = tokens[tokens.index('**elset') + 1]
      irot = tokens.index('*rotation')
      if tokens[irot+1] == 'x1':
        x1 = np.empty(3, dtype=float)
        x1[0] = float(tokens[irot + 2])
        x1[1] = float(tokens[irot + 3])
        x1[2] = float(tokens[irot + 4])
        x3 = np.empty(3, dtype=float)
        x3[0] = float(tokens[irot + 6])
        x3[1] = float(tokens[irot + 7])
        x3[2] = float(tokens[irot + 8])
        euler.append([elset, Orientation.Zrot2OrientationMatrix(x1=x1, x3=x3)])
      else: # euler angles
        phi1 = tokens[irot + 1]
        Phi = tokens[irot + 2]
        phi2 = tokens[irot + 3]
        angles = np.array([float(phi1), float(Phi), float(phi2)])
        euler.append([elset, Orientation.from_euler(angles)])
    return dict(euler)
  
  def slip_system_orientation_tensor(self, s):
    '''Compute the orientation strain tensor m^s for this :py:class:`~pymicro.crystal.microstructure.Orientation` 
    and the given slip system.
    
    :param s: an instance of :py:class:`~pymicro.crystal.lattice.SlipSystem`

    .. math::

      M^s_{ij} = \left(l^s_i.n^s_j)
    '''
    Bt = self.orientation_matrix().transpose()
    plane = s.get_slip_plane()
    n_rot = np.dot(Bt, plane.normal())
    slip = s.get_slip_direction()
    l_rot = np.dot(Bt, slip.direction())
    return np.outer(l_rot, n_rot)

  def slip_system_orientation_strain_tensor(self, s):
    '''Compute the orientation strain tensor m^s for this :py:class:`~pymicro.crystal.microstructure.Orientation` 
    and the given slip system.
    
    :param s: an instance of :py:class:`~pymicro.crystal.lattice.SlipSystem`

    .. math::

      m^s_{ij} = \\frac{1}{2}\left(l^s_i.n^s_j + l^s_j.n^s_i)
    '''
    Bt = self.orientation_matrix().transpose()
    plane = s.get_slip_plane()
    n_rot = np.dot(Bt, plane.normal())
    slip = s.get_slip_direction()
    l_rot = np.dot(Bt, slip.direction())
    m = 0.5 * (np.outer(l_rot, n_rot) + np.outer(n_rot, l_rot))
    return m

  def slip_system_orientation_rotation_tensor(self, s):
    '''Compute the orientation rotation tensor q^s for this :py:class:`~pymicro.crystal.microstructure.Orientation` 
    and the given slip system.
    
    :param s: an instance of :py:class:`~pymicro.crystal.lattice.SlipSystem`

    .. math::

      q^s_{ij} = \\frac{1}{2}\left(l^s_i.n^s_j - l^s_j.n^s_i)
    '''
    Bt = self.orientation_matrix().transpose()
    plane = s.get_slip_plane()
    n_rot = np.dot(Bt, plane.normal())
    slip = s.get_slip_direction()
    l_rot = np.dot(Bt, slip.direction())
    q = 0.5 * (np.outer(l_rot, n_rot) - np.outer(n_rot, l_rot))
    return q

  def schmid_factor(self, slip_system, load_direction=[0., 0., 1]):
    '''Compute the Schmid factor for this crystal orientation and the 
    given slip system.

    :param slip_system: a slip system instance.
    :param load_direction: a unit vector describing the loading direction (default: vertical axis [0, 0, 1]).
    :returns float: a number between 0 ad 0.5.
    '''
    plane = slip_system.get_slip_plane()
    Bt = self.orientation_matrix().transpose()
    n_rot = np.dot(Bt, plane.normal()) # plane.normal() is a unit vector
    slip = slip_system.get_slip_direction().direction()
    slip_rot = np.dot(Bt, slip)
    #print n_rot, load_direction, slip_rot
    SF = np.abs(np.dot(n_rot, load_direction) * np.dot(slip_rot, load_direction))
    return SF

  def compute_all_schmid_factors(self, slip_systems, load_direction=[0., 0., 1], verbose=False):
    '''Compute all Schmid factors for this crystal orientation and the 
    given list of slip systems.

    :param slip_systems: a list of the slip system from which to compute the Schmid factor values.
    :param load_direction: a unit vector describing the loading direction (default: vertical axis [0, 0, 1]).
    :param bool verbose: activate verbose mode.
    :returns list: a list of the schmid factors.
    '''
    SF_list = []
    for ss in slip_systems:
      print ss
      sf = self.schmid_factor(ss, load_direction)
      if verbose:
        print 'Slip system: %s, Schmid factor is %.3f' % (ss, sf)
      SF_list.append(sf)
    return SF_list

class Grain:
  '''
  Class defining a crystallographic grain.

  A grain has its own crystallographic orientation.
  An optional id for the grain may be specified.
  The position field is the center of mass of the grain in world coordinates.
  The volume of the grain is expressed in pixel/voxel unit.
  '''

  def __init__(self, grain_id, grain_orientation):
    self.id = grain_id
    self.orientation = grain_orientation
    self.position = np.array([0., 0., 0.])
    self.volume = 0 # warning not implemented
    self.vtkmesh = None

  def __repr__(self):
    '''Provide a string representation of the class.'''
    s = '%s\n * id = %d\n' % (self.__class__.__name__, self.id)
    s += ' * %s\n' % (self.orientation)
    s += ' * position %s\n' % np.array_str(self.position)
    s += ' * has vtk mesh ? %s\n' % (self.vtkmesh != None)
    return s
  
  def schmid_factor(self, slip_system, load_direction=[0., 0., 1]):
    '''Compute the Schmid factor of this grain for the given slip system.

    **Parameters**:
    
    *slip_system*: a slip system instance.
    
    *load_direction*: a unit vector describing the loading direction.
    
    **Returns**
    
    The Schmid factor of this grain for the given slip system.
    '''
{'fraction_non_alphanumeric': 0.1396629551609991, 'fraction_numerical': 0.03930183569064099, 'mean_word_length': 2.9419335705812575, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 39, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 6}, 'pii_count': 0, 'substrings_counts': 4, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '630452', 'n_tokens_mistral': 25503, 'n_tokens_neox': 23553, 'n_words': 6314}