documenting (in accordance with new prospector rules)

This commit is contained in:
Martin Diehl 2019-09-03 16:52:28 -07:00
parent a428a924eb
commit 3657f81c59
1 changed files with 257 additions and 122 deletions

View File

@ -10,35 +10,44 @@ from .quaternion import P
#################################################################################################### ####################################################################################################
class Rotation: class Rotation:
u""" u"""
Orientation stored as unit quaternion. Orientation stored with functionality for conversion to different representations.
Following: D Rowenhorst et al. Consistent representations of and conversions between 3D rotations References
10.1088/0965-0393/23/8/083501 ----------
Convention 1: coordinate frames are right-handed D. Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015
Convention 2: a rotation angle ω is taken to be positive for a counterclockwise rotation https://doi.org/10.1088/0965-0393/23/8/083501
when viewing from the end point of the rotation axis towards the origin
Convention 3: rotations will be interpreted in the passive sense Conventions
-----------
Convention 1: Coordinate frames are right-handed.
Convention 2: A rotation angle ω is taken to be positive for a counterclockwise rotation
when viewing from the end point of the rotation axis towards the origin.
Convention 3: Rotations will be interpreted in the passive sense.
Convention 4: Euler angle triplets are implemented using the Bunge convention, Convention 4: Euler angle triplets are implemented using the Bunge convention,
with the angular ranges as [0, 2π],[0, π],[0, 2π] with the angular ranges as [0, 2π],[0, π],[0, 2π].
Convention 5: the rotation angle ω is limited to the interval [0, π] Convention 5: The rotation angle ω is limited to the interval [0, π].
Convention 6: P = -1 (as default) Convention 6: P = -1 (as default).
q is the real part, p = (x, y, z) are the imaginary parts.
Usage
-----
Vector "a" (defined in coordinate system "A") is passively rotated Vector "a" (defined in coordinate system "A") is passively rotated
resulting in new coordinates "b" when expressed in system "B". resulting in new coordinates "b" when expressed in system "B".
b = Q * a b = Q * a
b = np.dot(Q.asMatrix(),a) b = np.dot(Q.asMatrix(),a)
""" """
__slots__ = ['quaternion'] __slots__ = ['quaternion']
def __init__(self,quaternion = np.array([1.0,0.0,0.0,0.0])): def __init__(self,quaternion = np.array([1.0,0.0,0.0,0.0])):
""" """
Initializes to identity unless specified Initializes to identity unless specified.
Parameters
----------
quaternion : numpy.ndarray, optional
Unit quaternion that follows the conventions. Use .fromQuaternion to perform a sanity check.
If a quaternion is given, it needs to comply with the convection. Use .fromQuaternion
to check the input.
""" """
if isinstance(quaternion,Quaternion): if isinstance(quaternion,Quaternion):
self.quaternion = quaternion.copy() self.quaternion = quaternion.copy()
@ -46,14 +55,14 @@ class Rotation:
self.quaternion = Quaternion(q=quaternion[0],p=quaternion[1:4]) self.quaternion = Quaternion(q=quaternion[0],p=quaternion[1:4])
def __copy__(self): def __copy__(self):
"""Copy""" """Copy."""
return self.__class__(self.quaternion) return self.__class__(self.quaternion)
copy = __copy__ copy = __copy__
def __repr__(self): def __repr__(self):
"""Value in selected representation""" """Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles."""
return '\n'.join([ return '\n'.join([
'{}'.format(self.quaternion), '{}'.format(self.quaternion),
'Matrix:\n{}'.format( '\n'.join(['\t'.join(list(map(str,self.asMatrix()[i,:]))) for i in range(3)]) ), 'Matrix:\n{}'.format( '\n'.join(['\t'.join(list(map(str,self.asMatrix()[i,:]))) for i in range(3)]) ),
@ -62,9 +71,18 @@ class Rotation:
def __mul__(self, other): def __mul__(self, other):
""" """
Multiplication Multiplication.
Parameters
----------
other : numpy.ndarray or Rotation
Vector, second or fourth order tensor, or rotation object that is rotated.
Todo
----
Document details active/passive)
considere rotation of (3,3,3,3)-matrix
Rotation: Details needed (active/passive), rotation of (3,3,3,3)-matrix should be considered
""" """
if isinstance(other, Rotation): # rotate a rotation if isinstance(other, Rotation): # rotate a rotation
return self.__class__(self.quaternion * other.quaternion).standardize() return self.__class__(self.quaternion * other.quaternion).standardize()
@ -104,32 +122,48 @@ class Rotation:
def inverse(self): def inverse(self):
"""In-place inverse rotation/backward rotation""" """In-place inverse rotation/backward rotation."""
self.quaternion.conjugate() self.quaternion.conjugate()
return self return self
def inversed(self): def inversed(self):
"""Inverse rotation/backward rotation""" """Inverse rotation/backward rotation."""
return self.copy().inverse() return self.copy().inverse()
def standardize(self): def standardize(self):
"""In-place quaternion representation with positive q""" """In-place quaternion representation with positive q."""
if self.quaternion.q < 0.0: self.quaternion.homomorph() if self.quaternion.q < 0.0: self.quaternion.homomorph()
return self return self
def standardized(self): def standardized(self):
"""Quaternion representation with positive q""" """Quaternion representation with positive q."""
return self.copy().standardize() return self.copy().standardize()
def misorientation(self,other): def misorientation(self,other):
"""Misorientation""" """
Get Misorientation.
Parameters
----------
other : Rotation
Rotation to which the misorientation is computed.
"""
return other*self.inversed() return other*self.inversed()
def average(self,other): def average(self,other):
"""Calculate the average rotation""" """
Calculate the average rotation.
Parameters
----------
other : Rotation
Rotation from which the average is rotated.
"""
return Rotation.fromAverage([self,other]) return Rotation.fromAverage([self,other])
@ -137,41 +171,75 @@ class Rotation:
# convert to different orientation representations (numpy arrays) # convert to different orientation representations (numpy arrays)
def asQuaternion(self): def asQuaternion(self):
"""Unit quaternion: (q, [p_1, p_2, p_3])""" """Unit quaternion: (q, p_1, p_2, p_3)."""
return self.quaternion.asArray() return self.quaternion.asArray()
def asEulers(self, def asEulers(self,
degrees = False): degrees = False):
"""Bunge-Euler angles: (φ_1, ϕ, φ_2)""" """
Bunge-Euler angles: (φ_1, ϕ, φ_2).
Parameters
----------
degrees : bool, optional
return angles in degrees.
"""
eu = qu2eu(self.quaternion.asArray()) eu = qu2eu(self.quaternion.asArray())
if degrees: eu = np.degrees(eu) if degrees: eu = np.degrees(eu)
return eu return eu
def asAxisAngle(self, def asAxisAngle(self,
degrees = False): degrees = False):
"""Axis-angle pair: ([n_1, n_2, n_3], ω)""" """
Axis angle pair: ([n_1, n_2, n_3], ω).
Parameters
----------
degrees : bool, optional
return rotation angle in degrees.
"""
ax = qu2ax(self.quaternion.asArray()) ax = qu2ax(self.quaternion.asArray())
if degrees: ax[3] = np.degrees(ax[3]) if degrees: ax[3] = np.degrees(ax[3])
return ax return ax
def asMatrix(self): def asMatrix(self):
"""Rotation matrix""" """Rotation matrix."""
return qu2om(self.quaternion.asArray()) return qu2om(self.quaternion.asArray())
def asRodrigues(self,vector=False): def asRodrigues(self,
"""Rodrigues-Frank vector: ([n_1, n_2, n_3], tan(ω/2))""" vector=False):
"""
Rodrigues-Frank vector: ([n_1, n_2, n_3], tan(ω/2)).
Parameters
----------
vector : bool, optional
return as array of length 3, i.e. scale the unit vector giving the rotation axis.
"""
ro = qu2ro(self.quaternion.asArray()) ro = qu2ro(self.quaternion.asArray())
return ro[:3]*ro[3] if vector else ro return ro[:3]*ro[3] if vector else ro
def asHomochoric(self): def asHomochoric(self):
"""Homochoric vector: (h_1, h_2, h_3)""" """Homochoric vector: (h_1, h_2, h_3)."""
return qu2ho(self.quaternion.asArray()) return qu2ho(self.quaternion.asArray())
def asCubochoric(self): def asCubochoric(self):
"""Cubochoric vector: (c_1, c_2, c_3)."""
return qu2cu(self.quaternion.asArray()) return qu2cu(self.quaternion.asArray())
def asM(self): def asM(self):
"""Intermediate representation supporting quaternion averaging (see F. Landis Markley et al.)""" """
Intermediate representation supporting quaternion averaging.
References
----------
F. Landis Markley et al., Journal of Guidance, Control, and Dynamics 30(4):1193-1197, 2007
https://doi.org/10.2514/1.28949
"""
return self.quaternion.asM() return self.quaternion.asM()
@ -299,14 +367,20 @@ class Rotation:
rotations, rotations,
weights = []): weights = []):
""" """
Average rotation Average rotation.
ref: F. Landis Markley, Yang Cheng, John Lucas Crassidis, and Yaakov Oshman. References
Averaging Quaternions, ----------
Journal of Guidance, Control, and Dynamics, Vol. 30, No. 4 (2007), pp. 1193-1197. F. Landis Markley et al., Journal of Guidance, Control, and Dynamics 30(4):1193-1197, 2007
doi: 10.2514/1.28949 https://doi.org/10.2514/1.28949
Parameters
----------
rotations : list of Rotations
Rotations to average from
weights : list of floats, optional
Weights for each rotation used for averaging
usage: input a list of rotations and optional weights
""" """
if not all(isinstance(item, Rotation) for item in rotations): if not all(isinstance(item, Rotation) for item in rotations):
raise TypeError("Only instances of Rotation can be averaged.") raise TypeError("Only instances of Rotation can be averaged.")
@ -339,42 +413,78 @@ class Rotation:
# ****************************************************************************************** # ******************************************************************************************
class Symmetry: class Symmetry:
""" """
Symmetry operations for lattice systems Symmetry operations for lattice systems.
References
----------
https://en.wikipedia.org/wiki/Crystal_system https://en.wikipedia.org/wiki/Crystal_system
""" """
lattices = [None,'orthorhombic','tetragonal','hexagonal','cubic',] lattices = [None,'orthorhombic','tetragonal','hexagonal','cubic',]
def __init__(self, symmetry = None): def __init__(self, symmetry = None):
if isinstance(symmetry, str) and symmetry.lower() in Symmetry.lattices: """
self.lattice = symmetry.lower() Symmetry Definition.
else:
self.lattice = None Parameters
----------
symmetry : str, optional
label of the crystal system
"""
if symmetry is not None and symmetry.lower() not in Symmetry.lattices:
raise KeyError('Symmetry/crystal system "{}" is unknown'.format(symmetry))
self.lattice = symmetry.lower() if isinstance(symmetry,str) else symmetry
def __copy__(self): def __copy__(self):
"""Copy""" """Copy."""
return self.__class__(self.lattice) return self.__class__(self.lattice)
copy = __copy__ copy = __copy__
def __repr__(self): def __repr__(self):
"""Readable string""" """Readable string."""
return '{}'.format(self.lattice) return '{}'.format(self.lattice)
def __eq__(self, other): def __eq__(self, other):
"""Equal to other""" """
Equal to other.
Parameters
----------
other : Symmetry
Symmetry to check for equality.
"""
return self.lattice == other.lattice return self.lattice == other.lattice
def __neq__(self, other): def __neq__(self, other):
"""Not equal to other""" """
Not Equal to other.
Parameters
----------
other : Symmetry
Symmetry to check for inequality.
"""
return not self.__eq__(other) return not self.__eq__(other)
def __cmp__(self,other): def __cmp__(self,other):
"""Linear ordering""" """
Linear ordering.
Parameters
----------
other : Symmetry
Symmetry to check for for order.
"""
myOrder = Symmetry.lattices.index(self.lattice) myOrder = Symmetry.lattices.index(self.lattice)
otherOrder = Symmetry.lattices.index(other.lattice) otherOrder = Symmetry.lattices.index(other.lattice)
return (myOrder > otherOrder) - (myOrder < otherOrder) return (myOrder > otherOrder) - (myOrder < otherOrder)
@ -458,7 +568,7 @@ class Symmetry:
def inFZ(self,rodrigues): def inFZ(self,rodrigues):
""" """
Check whether given Rodrigues vector falls into fundamental zone of own symmetry. Check whether given Rodriques-Frank vector falls into fundamental zone of own symmetry.
Fundamental zone in Rodrigues space is point symmetric around origin. Fundamental zone in Rodrigues space is point symmetric around origin.
""" """
@ -491,11 +601,13 @@ class Symmetry:
def inDisorientationSST(self,rodrigues): def inDisorientationSST(self,rodrigues):
""" """
Check whether given Rodrigues vector (of misorientation) falls into standard stereographic triangle of own symmetry. Check whether given Rodriques-Frank vector (of misorientation) falls into standard stereographic triangle of own symmetry.
References
----------
A. Heinz and P. Neumann, Acta Crystallographica Section A 47:780-789, 1991
https://doi.org/10.1107/S0108767391006864
Determination of disorientations follow the work of A. Heinz and P. Neumann:
Representation of Orientation and Disorientation Data for Cubic, Hexagonal, Tetragonal and Orthorhombic Crystals
Acta Cryst. (1991). A47, 780-789
""" """
if (len(rodrigues) != 3): if (len(rodrigues) != 3):
raise ValueError('Input is not a Rodriques-Frank vector.\n') raise ValueError('Input is not a Rodriques-Frank vector.\n')
@ -611,11 +723,15 @@ class Symmetry:
# ****************************************************************************************** # ******************************************************************************************
class Lattice: class Lattice:
""" """
Lattice system Lattice system.
Currently, this contains only a mapping from Bravais lattice to symmetry Currently, this contains only a mapping from Bravais lattice to symmetry
and orientation relationships. It could include twin and slip systems. and orientation relationships. It could include twin and slip systems.
References
----------
https://en.wikipedia.org/wiki/Bravais_lattice https://en.wikipedia.org/wiki/Bravais_lattice
""" """
lattices = { lattices = {
@ -628,12 +744,21 @@ class Lattice:
def __init__(self, lattice): def __init__(self, lattice):
"""
New lattice of given type.
Parameters
----------
lattice : str
Bravais lattice.
"""
self.lattice = lattice self.lattice = lattice
self.symmetry = Symmetry(self.lattices[lattice]['symmetry']) self.symmetry = Symmetry(self.lattices[lattice]['symmetry'])
def __repr__(self): def __repr__(self):
"""Report basic lattice information""" """Report basic lattice information."""
return 'Bravais lattice {} ({} symmetry)'.format(self.lattice,self.symmetry) return 'Bravais lattice {} ({} symmetry)'.format(self.lattice,self.symmetry)
@ -878,7 +1003,7 @@ class Lattice:
'NW':self.NW, 'Pitsch': self.Pitsch, 'Bain':self.Bain} 'NW':self.NW, 'Pitsch': self.Pitsch, 'Bain':self.Bain}
try: try:
relationship = models[model] relationship = models[model]
except: except KeyError :
raise KeyError('Orientation relationship "{}" is unknown'.format(model)) raise KeyError('Orientation relationship "{}" is unknown'.format(model))
if self.lattice not in relationship['mapping']: if self.lattice not in relationship['mapping']:
@ -909,19 +1034,29 @@ class Lattice:
class Orientation: class Orientation:
""" """
Crystallographic orientation Crystallographic orientation.
A crystallographic orientation contains a rotation and a lattice A crystallographic orientation contains a rotation and a lattice.
""" """
__slots__ = ['rotation','lattice'] __slots__ = ['rotation','lattice']
def __repr__(self): def __repr__(self):
"""Report lattice type and orientation""" """Report lattice type and orientation."""
return self.lattice.__repr__()+'\n'+self.rotation.__repr__() return self.lattice.__repr__()+'\n'+self.rotation.__repr__()
def __init__(self, rotation, lattice): def __init__(self, rotation, lattice):
"""
New orientation from rotation and lattice.
Parameters
----------
rotation : Rotation
Rotation specifying the lattice orientation.
lattice : Lattice
Lattice type of the crystal.
"""
if isinstance(lattice, Lattice): if isinstance(lattice, Lattice):
self.lattice = lattice self.lattice = lattice
else: else:
@ -971,7 +1106,7 @@ class Orientation:
return self.lattice.symmetry.inFZ(self.rotation.asRodrigues(vector=True)) return self.lattice.symmetry.inFZ(self.rotation.asRodrigues(vector=True))
def equivalentOrientations(self,members=[]): def equivalentOrientations(self,members=[]):
"""List of orientations which are symmetrically equivalent""" """List of orientations which are symmetrically equivalent."""
try: try:
iter(members) # asking for (even empty) list of members? iter(members) # asking for (even empty) list of members?
except TypeError: except TypeError:
@ -981,12 +1116,12 @@ class Orientation:
for q in self.lattice.symmetry.symmetryOperations(members)] # yes, return list of rotations for q in self.lattice.symmetry.symmetryOperations(members)] # yes, return list of rotations
def relatedOrientations(self,model): def relatedOrientations(self,model):
"""List of orientations related by the given orientation relationship""" """List of orientations related by the given orientation relationship."""
r = self.lattice.relationOperations(model) r = self.lattice.relationOperations(model)
return [self.__class__(self.rotation*o,r['lattice']) for o in r['rotations']] return [self.__class__(self.rotation*o,r['lattice']) for o in r['rotations']]
def reduced(self): def reduced(self):
"""Transform orientation to fall into fundamental zone according to symmetry""" """Transform orientation to fall into fundamental zone according to symmetry."""
for me in self.equivalentOrientations(): for me in self.equivalentOrientations():
if self.lattice.symmetry.inFZ(me.rotation.asRodrigues(vector=True)): break if self.lattice.symmetry.inFZ(me.rotation.asRodrigues(vector=True)): break
@ -996,7 +1131,7 @@ class Orientation:
axis, axis,
proper = False, proper = False,
SST = True): SST = True):
"""Axis rotated according to orientation (using crystal symmetry to ensure location falls into SST)""" """Axis rotated according to orientation (using crystal symmetry to ensure location falls into SST)."""
if SST: # pole requested to be within SST if SST: # pole requested to be within SST
for i,o in enumerate(self.equivalentOrientations()): # test all symmetric equivalent quaternions for i,o in enumerate(self.equivalentOrientations()): # test all symmetric equivalent quaternions
pole = o.rotation*axis # align crystal direction to axis pole = o.rotation*axis # align crystal direction to axis
@ -1008,7 +1143,7 @@ class Orientation:
def IPFcolor(self,axis): def IPFcolor(self,axis):
"""TSL color of inverse pole figure for given axis""" """TSL color of inverse pole figure for given axis."""
color = np.zeros(3,'d') color = np.zeros(3,'d')
for o in self.equivalentOrientations(): for o in self.equivalentOrientations():
@ -1023,7 +1158,7 @@ class Orientation:
def fromAverage(cls, def fromAverage(cls,
orientations, orientations,
weights = []): weights = []):
"""Create orientation from average of list of orientations""" """Create orientation from average of list of orientations."""
if not all(isinstance(item, Orientation) for item in orientations): if not all(isinstance(item, Orientation) for item in orientations):
raise TypeError("Only instances of Orientation can be averaged.") raise TypeError("Only instances of Orientation can be averaged.")
@ -1039,7 +1174,7 @@ class Orientation:
def average(self,other): def average(self,other):
"""Calculate the average rotation""" """Calculate the average rotation."""
return Orientation.fromAverage([self,other]) return Orientation.fromAverage([self,other])
@ -1080,10 +1215,10 @@ def isone(a):
def iszero(a): def iszero(a):
return np.isclose(a,0.0,atol=1.0e-12,rtol=0.0) return np.isclose(a,0.0,atol=1.0e-12,rtol=0.0)
#---------- quaternion ---------- #---------- Quaternion ----------
def qu2om(qu): def qu2om(qu):
"""Quaternion to orientation matrix""" """Quaternion to rotation matrix."""
qq = qu[0]**2-(qu[1]**2 + qu[2]**2 + qu[3]**2) qq = qu[0]**2-(qu[1]**2 + qu[2]**2 + qu[3]**2)
om = np.diag(qq + 2.0*np.array([qu[1],qu[2],qu[3]])**2) om = np.diag(qq + 2.0*np.array([qu[1],qu[2],qu[3]])**2)
@ -1097,7 +1232,7 @@ def qu2om(qu):
def qu2eu(qu): def qu2eu(qu):
"""Quaternion to Euler angles""" """Quaternion to Bunge-Euler angles."""
q03 = qu[0]**2+qu[3]**2 q03 = qu[0]**2+qu[3]**2
q12 = qu[1]**2+qu[2]**2 q12 = qu[1]**2+qu[2]**2
chi = np.sqrt(q03*q12) chi = np.sqrt(q03*q12)
@ -1117,7 +1252,7 @@ def qu2eu(qu):
def qu2ax(qu): def qu2ax(qu):
""" """
Quaternion to axis angle Quaternion to axis angle pair.
Modified version of the original formulation, should be numerically more stable Modified version of the original formulation, should be numerically more stable
""" """
@ -1134,7 +1269,7 @@ def qu2ax(qu):
def qu2ro(qu): def qu2ro(qu):
"""Quaternion to Rodrigues vector""" """Quaternion to Rodriques-Frank vector."""
if iszero(qu[0]): if iszero(qu[0]):
ro = [qu[1], qu[2], qu[3], np.inf] ro = [qu[1], qu[2], qu[3], np.inf]
else: else:
@ -1146,7 +1281,7 @@ def qu2ro(qu):
def qu2ho(qu): def qu2ho(qu):
"""Quaternion to homochoric""" """Quaternion to homochoric vector."""
omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) # avoid numerical difficulties omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) # avoid numerical difficulties
if iszero(omega): if iszero(omega):
@ -1160,15 +1295,15 @@ def qu2ho(qu):
def qu2cu(qu): def qu2cu(qu):
"""Quaternion to cubochoric""" """Quaternion to cubochoric vector."""
return ho2cu(qu2ho(qu)) return ho2cu(qu2ho(qu))
#---------- orientation matrix ---------- #---------- Rotation matrix ----------
def om2qu(om): def om2qu(om):
""" """
Orientation matrix to quaternion Rotation matrix to quaternion.
The original formulation (direct conversion) had (numerical?) issues The original formulation (direct conversion) had (numerical?) issues
""" """
@ -1176,7 +1311,7 @@ def om2qu(om):
def om2eu(om): def om2eu(om):
"""Orientation matrix to Euler angles""" """Rotation matrix to Bunge-Euler angles."""
if abs(om[2,2]) < 1.0: if abs(om[2,2]) < 1.0:
zeta = 1.0/np.sqrt(1.0-om[2,2]**2) zeta = 1.0/np.sqrt(1.0-om[2,2]**2)
eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta), eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta),
@ -1191,7 +1326,7 @@ def om2eu(om):
def om2ax(om): def om2ax(om):
"""Orientation matrix to axis angle""" """Rotation matrix to axis angle pair."""
ax=np.empty(4) ax=np.empty(4)
# first get the rotation angle # first get the rotation angle
@ -1212,24 +1347,24 @@ def om2ax(om):
def om2ro(om): def om2ro(om):
"""Orientation matrix to Rodriques vector""" """Rotation matrix to Rodriques-Frank vector."""
return eu2ro(om2eu(om)) return eu2ro(om2eu(om))
def om2ho(om): def om2ho(om):
"""Orientation matrix to homochoric""" """Rotation matrix to homochoric vector."""
return ax2ho(om2ax(om)) return ax2ho(om2ax(om))
def om2cu(om): def om2cu(om):
"""Orientation matrix to cubochoric""" """Rotation matrix to cubochoric vector."""
return ho2cu(om2ho(om)) return ho2cu(om2ho(om))
#---------- Euler angles ---------- #---------- Bunge-Euler angles ----------
def eu2qu(eu): def eu2qu(eu):
"""Euler angles to quaternion""" """Bunge-Euler angles to quaternion."""
ee = 0.5*eu ee = 0.5*eu
cPhi = np.cos(ee[1]) cPhi = np.cos(ee[1])
sPhi = np.sin(ee[1]) sPhi = np.sin(ee[1])
@ -1242,7 +1377,7 @@ def eu2qu(eu):
def eu2om(eu): def eu2om(eu):
"""Euler angles to orientation matrix""" """Bunge-Euler angles to rotation matrix."""
c = np.cos(eu) c = np.cos(eu)
s = np.sin(eu) s = np.sin(eu)
@ -1255,7 +1390,7 @@ def eu2om(eu):
def eu2ax(eu): def eu2ax(eu):
"""Euler angles to axis angle""" """Bunge-Euler angles to axis angle pair."""
t = np.tan(eu[1]*0.5) t = np.tan(eu[1]*0.5)
sigma = 0.5*(eu[0]+eu[2]) sigma = 0.5*(eu[0]+eu[2])
delta = 0.5*(eu[0]-eu[2]) delta = 0.5*(eu[0]-eu[2])
@ -1266,7 +1401,7 @@ def eu2ax(eu):
if iszero(alpha): if iszero(alpha):
ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) ax = np.array([ 0.0, 0.0, 1.0, 0.0 ])
else: else:
ax = -P/tau * np.array([ t*np.cos(delta), t*np.sin(delta), np.sin(sigma) ]) # passive axis-angle pair so a minus sign in front ax = -P/tau * np.array([ t*np.cos(delta), t*np.sin(delta), np.sin(sigma) ]) # passive axis angle pair so a minus sign in front
ax = np.append(ax,alpha) ax = np.append(ax,alpha)
if alpha < 0.0: ax *= -1.0 # ensure alpha is positive if alpha < 0.0: ax *= -1.0 # ensure alpha is positive
@ -1274,8 +1409,8 @@ def eu2ax(eu):
def eu2ro(eu): def eu2ro(eu):
"""Euler angles to Rodrigues vector""" """Bunge-Euler angles to Rodriques-Frank vector."""
ro = eu2ax(eu) # convert to axis angle representation ro = eu2ax(eu) # convert to axis angle pair representation
if ro[3] >= np.pi: # Differs from original implementation. check convention 5 if ro[3] >= np.pi: # Differs from original implementation. check convention 5
ro[3] = np.inf ro[3] = np.inf
elif iszero(ro[3]): elif iszero(ro[3]):
@ -1287,19 +1422,19 @@ def eu2ro(eu):
def eu2ho(eu): def eu2ho(eu):
"""Euler angles to homochoric""" """Bunge-Euler angles to homochoric vector."""
return ax2ho(eu2ax(eu)) return ax2ho(eu2ax(eu))
def eu2cu(eu): def eu2cu(eu):
"""Euler angles to cubochoric""" """Bunge-Euler angles to cubochoric vector."""
return ho2cu(eu2ho(eu)) return ho2cu(eu2ho(eu))
#---------- axis angle ---------- #---------- Axis angle pair ----------
def ax2qu(ax): def ax2qu(ax):
"""Axis angle to quaternion""" """Axis angle pair to quaternion."""
if iszero(ax[3]): if iszero(ax[3]):
qu = np.array([ 1.0, 0.0, 0.0, 0.0 ]) qu = np.array([ 1.0, 0.0, 0.0, 0.0 ])
else: else:
@ -1311,7 +1446,7 @@ def ax2qu(ax):
def ax2om(ax): def ax2om(ax):
"""Axis angle to orientation matrix""" """Axis angle pair to rotation matrix."""
c = np.cos(ax[3]) c = np.cos(ax[3])
s = np.sin(ax[3]) s = np.sin(ax[3])
omc = 1.0-c omc = 1.0-c
@ -1326,12 +1461,12 @@ def ax2om(ax):
def ax2eu(ax): def ax2eu(ax):
"""Orientation matrix to Euler angles""" """Rotation matrix to Bunge Euler angles."""
return om2eu(ax2om(ax)) return om2eu(ax2om(ax))
def ax2ro(ax): def ax2ro(ax):
"""Axis angle to Rodrigues vector""" """Axis angle pair to Rodriques-Frank vector."""
if iszero(ax[3]): if iszero(ax[3]):
ro = [ 0.0, 0.0, P, 0.0 ] ro = [ 0.0, 0.0, P, 0.0 ]
else: else:
@ -1344,36 +1479,36 @@ def ax2ro(ax):
def ax2ho(ax): def ax2ho(ax):
"""Axis angle to homochoric""" """Axis angle pair to homochoric vector."""
f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0) f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0)
ho = ax[0:3] * f ho = ax[0:3] * f
return ho return ho
def ax2cu(ax): def ax2cu(ax):
"""Axis angle to cubochoric""" """Axis angle pair to cubochoric vector."""
return ho2cu(ax2ho(ax)) return ho2cu(ax2ho(ax))
#---------- Rodrigues--Frank ---------- #---------- Rodrigues-Frank vector ----------
def ro2qu(ro): def ro2qu(ro):
"""Rodrigues vector to quaternion""" """Rodriques-Frank vector to quaternion."""
return ax2qu(ro2ax(ro)) return ax2qu(ro2ax(ro))
def ro2om(ro): def ro2om(ro):
"""Rodgrigues vector to orientation matrix""" """Rodgrigues-Frank vector to rotation matrix."""
return ax2om(ro2ax(ro)) return ax2om(ro2ax(ro))
def ro2eu(ro): def ro2eu(ro):
"""Rodrigues vector to orientation matrix""" """Rodriques-Frank vector to Bunge-Euler angles."""
return om2eu(ro2om(ro)) return om2eu(ro2om(ro))
def ro2ax(ro): def ro2ax(ro):
"""Rodrigues vector to axis angle""" """Rodriques-Frank vector to axis angle pair."""
ta = ro[3] ta = ro[3]
if iszero(ta): if iszero(ta):
@ -1389,7 +1524,7 @@ def ro2ax(ro):
def ro2ho(ro): def ro2ho(ro):
"""Rodrigues vector to homochoric""" """Rodriques-Frank vector to homochoric vector."""
if iszero(np.sum(ro[0:3]**2.0)): if iszero(np.sum(ro[0:3]**2.0)):
ho = [ 0.0, 0.0, 0.0 ] ho = [ 0.0, 0.0, 0.0 ]
else: else:
@ -1400,29 +1535,29 @@ def ro2ho(ro):
def ro2cu(ro): def ro2cu(ro):
"""Rodrigues vector to cubochoric""" """Rodriques-Frank vector to cubochoric vector."""
return ho2cu(ro2ho(ro)) return ho2cu(ro2ho(ro))
#---------- homochoric ---------- #---------- Homochoric vector----------
def ho2qu(ho): def ho2qu(ho):
"""Homochoric to quaternion""" """Homochoric vector to quaternion."""
return ax2qu(ho2ax(ho)) return ax2qu(ho2ax(ho))
def ho2om(ho): def ho2om(ho):
"""Homochoric to orientation matrix""" """Homochoric vector to rotation matrix."""
return ax2om(ho2ax(ho)) return ax2om(ho2ax(ho))
def ho2eu(ho): def ho2eu(ho):
"""Homochoric to Euler angles""" """Homochoric vector to Bunge-Euler angles."""
return ax2eu(ho2ax(ho)) return ax2eu(ho2ax(ho))
def ho2ax(ho): def ho2ax(ho):
"""Homochoric to axis angle""" """Homochoric vector to axis angle pair."""
tfit = np.array([+1.0000000000018852, -0.5000000002194847, tfit = np.array([+1.0000000000018852, -0.5000000002194847,
-0.024999992127593126, -0.003928701544781374, -0.024999992127593126, -0.003928701544781374,
-0.0008152701535450438, -0.0002009500426119712, -0.0008152701535450438, -0.0002009500426119712,
@ -1448,42 +1583,42 @@ def ho2ax(ho):
def ho2ro(ho): def ho2ro(ho):
"""Axis angle to Rodriques vector""" """Axis angle pair to Rodriques-Frank vector."""
return ax2ro(ho2ax(ho)) return ax2ro(ho2ax(ho))
def ho2cu(ho): def ho2cu(ho):
"""Homochoric to cubochoric""" """Homochoric vector to cubochoric vector."""
return Lambert.BallToCube(ho) return Lambert.BallToCube(ho)
#---------- cubochoric ---------- #---------- Cubochoric ----------
def cu2qu(cu): def cu2qu(cu):
"""Cubochoric to quaternion""" """Cubochoric vector to quaternion."""
return ho2qu(cu2ho(cu)) return ho2qu(cu2ho(cu))
def cu2om(cu): def cu2om(cu):
"""Cubochoric to orientation matrix""" """Cubochoric vector to rotation matrix."""
return ho2om(cu2ho(cu)) return ho2om(cu2ho(cu))
def cu2eu(cu): def cu2eu(cu):
"""Cubochoric to Euler angles""" """Cubochoric vector to Bunge-Euler angles."""
return ho2eu(cu2ho(cu)) return ho2eu(cu2ho(cu))
def cu2ax(cu): def cu2ax(cu):
"""Cubochoric to axis angle""" """Cubochoric vector to axis angle pair."""
return ho2ax(cu2ho(cu)) return ho2ax(cu2ho(cu))
def cu2ro(cu): def cu2ro(cu):
"""Cubochoric to Rodrigues vector""" """Cubochoric vector to Rodriques-Frank vector."""
return ho2ro(cu2ho(cu)) return ho2ro(cu2ho(cu))
def cu2ho(cu): def cu2ho(cu):
"""Cubochoric to homochoric""" """Cubochoric vector to homochoric vector."""
return Lambert.CubeToBall(cu) return Lambert.CubeToBall(cu)