source: trunk/GSASIIlattice.py @ 3308

Last change on this file since 3308 was 3136, checked in by vondreele, 8 years ago

make GSAS-II python 3.6 compliant & preserve python 2.7 use;changes:
do from future import division, print_function for all GSAS-II py sources
all menu items revised to be py 2.7/3.6 compliant
all wx.OPEN --> wx.FD_OPEN in file dialogs
all integer divides (typically for image pixel math) made explicit with ; ambiguous ones made floats as appropriate
all print "stuff" --> print (stuff)
all print >> pFile,'stuff' --> pFile.writeCIFtemplate('stuff')
all read file opens made explicit 'r' or 'rb'
all cPickle imports made for py2.7 or 3.6 as cPickle or _pickle; test for '2' platform.version_tuple[0] for py 2.7
define cPickleload to select load(fp) or load(fp,encoding='latin-1') for loading gpx files; provides cross compatibility between py 2.7/3.6 gpx files
make dict.keys() as explicit list(dict.keys()) as needed (NB: possible source of remaining py3.6 bugs)
make zip(a,b) as explicit list(zip(a,b)) as needed (NB: possible source of remaining py3.6 bugs)
select unichr/chr according test for '2' platform.version_tuple[0] for py 2.7 (G2pwdGUI * G2plot) for special characters
select wg.EVT_GRID_CELL_CHANGE (classic) or wg.EVT_GRID_CELL_CHANGED (phoenix) in grid Bind
maxint --> maxsize; used in random number stuff
raise Exception,"stuff" --> raise Exception("stuff")
wx 'classic' sizer.DeleteWindows?() or 'phoenix' sizer.Clear(True)
wx 'classic' SetToolTipString?(text) or 'phoenix' SetToolTip?(wx.ToolTip?(text)); define SetToolTipString?(self,text) to handle the choice in plots
status.SetFields? --> status.SetStatusText?
'classic' AddSimpleTool? or 'phoenix' self.AddTool? for plot toolbar; Bind different as well
define GetItemPydata? as it doesn't exist in wx 'phoenix'
allow python versions 2.7 & 3.6 to run GSAS-II
Bind override commented out - no logging capability (NB: remove all logging code?)
all import ContentsValidator? open filename & test if valid then close; filepointer removed from Reader
binary importers (mostly images) test for 'byte' type & convert as needed to satisfy py 3.6 str/byte rules
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 98.3 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIlattice: Unit cells*
4---------------------------
5
6Perform lattice-related computations
7
8Note that *g* is the reciprocal lattice tensor, and *G* is its inverse,
9:math:`G = g^{-1}`, where
10
11  .. math::
12
13   G = \\left( \\begin{matrix}
14   a^2 & a b\\cos\gamma & a c\\cos\\beta \\\\
15   a b\\cos\\gamma & b^2 & b c \cos\\alpha \\\\
16   a c\\cos\\beta &  b c \\cos\\alpha & c^2
17   \\end{matrix}\\right)
18
19The "*A* tensor" terms are defined as
20:math:`A = (\\begin{matrix} G_{11} & G_{22} & G_{33} & 2G_{12} & 2G_{13} & 2G_{23}\\end{matrix})` and *A* can be used in this fashion:
21:math:`d^* = \sqrt {A_1 h^2 + A_2 k^2 + A_3 l^2 + A_4 hk + A_5 hl + A_6 kl}`, where
22*d* is the d-spacing, and :math:`d^*` is the reciprocal lattice spacing,
23:math:`Q = 2 \\pi d^* = 2 \\pi / d`
24'''
25########### SVN repository information ###################
26# $Date: 2017-10-23 16:39:16 +0000 (Mon, 23 Oct 2017) $
27# $Author: vondreele $
28# $Revision: 3136 $
29# $URL: trunk/GSASIIlattice.py $
30# $Id: GSASIIlattice.py 3136 2017-10-23 16:39:16Z vondreele $
31########### SVN repository information ###################
32from __future__ import division, print_function
33import math
34import copy
35import sys
36import random as ran
37import numpy as np
38import numpy.linalg as nl
39import GSASIIpath
40import GSASIImath as G2mth
41import GSASIIspc as G2spc
42import GSASIIElem as G2elem
43GSASIIpath.SetVersionNumber("$Revision: 3136 $")
44# trig functions in degrees
45sind = lambda x: np.sin(x*np.pi/180.)
46asind = lambda x: 180.*np.arcsin(x)/np.pi
47tand = lambda x: np.tan(x*np.pi/180.)
48atand = lambda x: 180.*np.arctan(x)/np.pi
49atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
50cosd = lambda x: np.cos(x*np.pi/180.)
51acosd = lambda x: 180.*np.arccos(x)/np.pi
52rdsq2d = lambda x,p: round(1.0/np.sqrt(x),p)
53rpd = np.pi/180.
54RSQ2PI = 1./np.sqrt(2.*np.pi)
55SQ2 = np.sqrt(2.)
56RSQPI = 1./np.sqrt(np.pi)
57R2pisq = 1./(2.*np.pi**2)
58nxs = np.newaxis
59
60def sec2HMS(sec):
61    """Convert time in sec to H:M:S string
62   
63    :param sec: time in seconds
64    :return: H:M:S string (to nearest 100th second)
65   
66    """
67    H = int(sec//3600)
68    M = int(sec//60-H*60)
69    S = sec-3600*H-60*M
70    return '%d:%2d:%.2f'%(H,M,S)
71   
72def rotdMat(angle,axis=0):
73    """Prepare rotation matrix for angle in degrees about axis(=0,1,2)
74
75    :param angle: angle in degrees
76    :param axis:  axis (0,1,2 = x,y,z) about which for the rotation
77    :return: rotation matrix - 3x3 numpy array
78
79    """
80    if axis == 2:
81        return np.array([[cosd(angle),-sind(angle),0],[sind(angle),cosd(angle),0],[0,0,1]])
82    elif axis == 1:
83        return np.array([[cosd(angle),0,-sind(angle)],[0,1,0],[sind(angle),0,cosd(angle)]])
84    else:
85        return np.array([[1,0,0],[0,cosd(angle),-sind(angle)],[0,sind(angle),cosd(angle)]])
86       
87def rotdMat4(angle,axis=0):
88    """Prepare rotation matrix for angle in degrees about axis(=0,1,2) with scaling for OpenGL
89
90    :param angle: angle in degrees
91    :param axis:  axis (0,1,2 = x,y,z) about which for the rotation
92    :return: rotation matrix - 4x4 numpy array (last row/column for openGL scaling)
93
94    """
95    Mat = rotdMat(angle,axis)
96    return np.concatenate((np.concatenate((Mat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
97   
98def fillgmat(cell):
99    """Compute lattice metric tensor from unit cell constants
100
101    :param cell: tuple with a,b,c,alpha, beta, gamma (degrees)
102    :return: 3x3 numpy array
103
104    """
105    a,b,c,alp,bet,gam = cell
106    g = np.array([
107        [a*a,  a*b*cosd(gam),  a*c*cosd(bet)],
108        [a*b*cosd(gam),  b*b,  b*c*cosd(alp)],
109        [a*c*cosd(bet) ,b*c*cosd(alp),   c*c]])
110    return g
111           
112def cell2Gmat(cell):
113    """Compute real and reciprocal lattice metric tensor from unit cell constants
114
115    :param cell: tuple with a,b,c,alpha, beta, gamma (degrees)
116    :return: reciprocal (G) & real (g) metric tensors (list of two numpy 3x3 arrays)
117
118    """
119    g = fillgmat(cell)
120    G = nl.inv(g)       
121    return G,g
122
123def A2Gmat(A,inverse=True):
124    """Fill real & reciprocal metric tensor (G) from A.
125
126    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
127    :param bool inverse: if True return both G and g; else just G
128    :return: reciprocal (G) & real (g) metric tensors (list of two numpy 3x3 arrays)
129
130    """
131    G = np.array([
132        [A[0],  A[3]/2.,  A[4]/2.], 
133        [A[3]/2.,A[1],    A[5]/2.], 
134        [A[4]/2.,A[5]/2.,    A[2]]])
135    if inverse:
136        g = nl.inv(G)
137        return G,g
138    else:
139        return G
140
141def Gmat2A(G):
142    """Extract A from reciprocal metric tensor (G)
143
144    :param G: reciprocal maetric tensor (3x3 numpy array
145    :return: A = [G11,G22,G33,2*G12,2*G13,2*G23]
146
147    """
148    return [G[0][0],G[1][1],G[2][2],2.*G[0][1],2.*G[0][2],2.*G[1][2]]
149   
150def cell2A(cell):
151    """Obtain A = [G11,G22,G33,2*G12,2*G13,2*G23] from lattice parameters
152
153    :param cell: [a,b,c,alpha,beta,gamma] (degrees)
154    :return: G reciprocal metric tensor as 3x3 numpy array
155
156    """
157    G,g = cell2Gmat(cell)
158    return Gmat2A(G)
159
160def A2cell(A):
161    """Compute unit cell constants from A
162
163    :param A: [G11,G22,G33,2*G12,2*G13,2*G23] G - reciprocal metric tensor
164    :return: a,b,c,alpha, beta, gamma (degrees) - lattice parameters
165
166    """
167    G,g = A2Gmat(A)
168    return Gmat2cell(g)
169
170def Gmat2cell(g):
171    """Compute real/reciprocal lattice parameters from real/reciprocal metric tensor (g/G)
172    The math works the same either way.
173
174    :param g (or G): real (or reciprocal) metric tensor 3x3 array
175    :return: a,b,c,alpha, beta, gamma (degrees) (or a*,b*,c*,alpha*,beta*,gamma* degrees)
176
177    """
178    oldset = np.seterr('raise')
179    a = np.sqrt(max(0,g[0][0]))
180    b = np.sqrt(max(0,g[1][1]))
181    c = np.sqrt(max(0,g[2][2]))
182    alp = acosd(g[2][1]/(b*c))
183    bet = acosd(g[2][0]/(a*c))
184    gam = acosd(g[0][1]/(a*b))
185    np.seterr(**oldset)
186    return a,b,c,alp,bet,gam
187
188def invcell2Gmat(invcell):
189    """Compute real and reciprocal lattice metric tensor from reciprocal
190       unit cell constants
191       
192    :param invcell: [a*,b*,c*,alpha*, beta*, gamma*] (degrees)
193    :return: reciprocal (G) & real (g) metric tensors (list of two 3x3 arrays)
194
195    """
196    G = fillgmat(invcell)
197    g = nl.inv(G)
198    return G,g
199
200def cellDijFill(pfx,phfx,SGData,parmDict): 
201    '''Returns the filled-out reciprocal cell (A) terms
202    from the parameter dictionaries corrected for Dij.
203
204    :param str pfx: parameter prefix ("n::", where n is a phase number)
205    :param dict SGdata: a symmetry object
206    :param dict parmDict: a dictionary of parameters
207
208    :returns: A,sigA where each is a list of six terms with the A terms
209    '''
210    if SGData['SGLaue'] in ['-1',]:
211        A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A1']+parmDict[phfx+'D22'],
212             parmDict[pfx+'A2']+parmDict[phfx+'D33'],
213             parmDict[pfx+'A3']+parmDict[phfx+'D12'],parmDict[pfx+'A4']+parmDict[phfx+'D13'],
214             parmDict[pfx+'A5']+parmDict[phfx+'D23']]
215    elif SGData['SGLaue'] in ['2/m',]:
216        if SGData['SGUniq'] == 'a':
217            A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A1']+parmDict[phfx+'D22'],
218                 parmDict[pfx+'A2']+parmDict[phfx+'D33'],0,0,parmDict[pfx+'A5']+parmDict[phfx+'D23']]
219        elif SGData['SGUniq'] == 'b':
220            A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A1']+parmDict[phfx+'D22'],
221                 parmDict[pfx+'A2']+parmDict[phfx+'D33'],0,parmDict[pfx+'A4']+parmDict[phfx+'D13'],0]
222        else:
223            A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A1']+parmDict[phfx+'D22'],
224                 parmDict[pfx+'A2']+parmDict[phfx+'D33'],parmDict[pfx+'A3']+parmDict[phfx+'D12'],0,0]
225    elif SGData['SGLaue'] in ['mmm',]:
226        A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A1']+parmDict[phfx+'D22'],
227             parmDict[pfx+'A2']+parmDict[phfx+'D33'],0,0,0]
228    elif SGData['SGLaue'] in ['4/m','4/mmm']:
229        A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A0']+parmDict[phfx+'D11'],
230             parmDict[pfx+'A2']+parmDict[phfx+'D33'],0,0,0]
231    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
232        A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A0']+parmDict[phfx+'D11'],
233             parmDict[pfx+'A2']+parmDict[phfx+'D33'],parmDict[pfx+'A0']+parmDict[phfx+'D11'],0,0]
234    elif SGData['SGLaue'] in ['3R', '3mR']:
235        A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A0']+parmDict[phfx+'D11'],
236            parmDict[pfx+'A0']+parmDict[phfx+'D11'],
237            parmDict[pfx+'A3']+parmDict[phfx+'D23'],parmDict[pfx+'A3']+parmDict[phfx+'D23'],
238            parmDict[pfx+'A3']+parmDict[phfx+'D23']]
239    elif SGData['SGLaue'] in ['m3m','m3']:
240        A = [parmDict[pfx+'A0']+parmDict[phfx+'D11'],parmDict[pfx+'A0']+parmDict[phfx+'D11'],
241             parmDict[pfx+'A0']+parmDict[phfx+'D11'],0,0,0]
242    return A
243   
244def prodMGMT(G,Mat):
245    '''Transform metric tensor by matrix
246   
247    :param G: array metric tensor
248    :param Mat: array transformation matrix
249    :return: array new metric tensor
250   
251    '''
252    return np.inner(Mat,np.inner(G,Mat).T)
253   
254def TransformCell(cell,Trans):
255    '''Transform lattice parameters by matrix
256   
257    :param cell: list a,b,c,alpha,beta,gamma,(volume)
258    :param Trans: array transformation matrix
259    :return: array transformed a,b,c,alpha,beta,gamma,volume
260   
261    '''
262    newCell = np.zeros(7)
263    g = cell2Gmat(cell)[1]
264    newg = prodMGMT(g,Trans)
265    newCell[:6] = Gmat2cell(newg)
266    newCell[6] = calc_V(cell2A(newCell[:6]))
267    return newCell
268   
269def TransformXYZ(XYZ,Trans,Vec):
270    return np.inner(XYZ,Trans)+Vec
271   
272def TransformU6(U6,Trans):
273    Uij = np.inner(Trans,np.inner(U6toUij(U6),Trans))
274    return UijtoU6(Uij)
275   
276def TransformPhase(oldPhase,newPhase,Trans,Uvec,Vvec,ifMag):
277    '''Transform atoms from oldPhase to newPhase by Trans & Vec
278   
279    :param oldPhase: dict G2 phase info for old phase
280    :param newPhase: dict G2 phase info for new phase; with new cell & space group
281            atoms are from oldPhase & will be transformed
282    :param Trans: array transformation matrix
283    :param Vec: array transformation vector
284    :param ifMag: bool True if convert to magnetic phase;
285        if True all nonmagnetic atoms will be removed
286       
287    '''
288   
289    cx,ct,cs,cia = oldPhase['General']['AtomPtrs']
290    cm = 0
291    if oldPhase['General']['Type'] == 'magnetic':
292        cm = cx+4
293    oAmat,oBmat = cell2AB(oldPhase['General']['Cell'][1:7])
294    nAmat,nBmat = cell2AB(newPhase['General']['Cell'][1:7])
295    SGData = newPhase['General']['SGData']
296    invTrans = nl.inv(Trans)
297    newAtoms,atCodes = FillUnitCell(oldPhase)
298#    GSASIIpath.IPyBreak()
299    Unit =[abs(int(max(unit))-1) for unit in Trans]
300    for i,unit in enumerate(Unit):
301        if unit > 0:
302            for j in range(unit):
303                moreAtoms = copy.deepcopy(newAtoms)
304                moreCodes = []
305                for atom,code in zip(moreAtoms,atCodes):
306                    atom[cx+i] += 1.
307                    if '+' in code:
308                        cell = list(eval(code.split('+')[1]))
309                        ops = code.split('+')[0]
310                    else:
311                        cell = [0,0,0]
312                        ops = code
313                    cell[i] += 1
314                    moreCodes.append('%s+%d,%d,%d'%(ops,cell[0],cell[1],cell[2])) 
315                newAtoms += moreAtoms
316                atCodes += moreCodes
317    if ifMag:
318        cia += 3
319        cs += 3
320        newPhase['General']['Type'] = 'magnetic'
321        newPhase['General']['AtomPtrs'] = [cx,ct,cs,cia]
322        magAtoms = []
323        magatCodes = []
324        Landeg = 2.0
325        for iat,atom in enumerate(newAtoms):
326            if len(G2elem.GetMFtable([atom[ct],],[Landeg,])):
327                magAtoms.append(atom[:cx+4]+[0.,0.,0.]+atom[cx+4:])
328                magatCodes.append(atCodes[iat])
329        newAtoms = magAtoms
330        atCodes = magatCodes
331        newPhase['Draw Atoms'] = []
332    for atom in newAtoms:
333        atom[cx:cx+3] = TransformXYZ(atom[cx:cx+3]-Uvec,invTrans,Vvec)%1.
334        if atom[cia] == 'A':
335            atom[cia+2:cia+8] = TransformU6(atom[cia+2:cia+8],invTrans)
336        atom[cs:cs+2] = G2spc.SytSym(atom[cx:cx+3],SGData)[:2]
337        atom[cia+8] = ran.randint(0,sys.maxsize)
338        if cm:
339            mag = np.sqrt(np.sum(np.array(atom[cm:cm+3])**2))
340            if mag:
341                mom = np.inner(np.array(atom[cm:cm+3]),oBmat)
342                mom = np.inner(mom,invTrans.T)
343                mom = np.inner(mom,nAmat)
344                mom /= np.sqrt(np.sum(mom**2))
345                atom[cm:cm+3] = mom*mag
346    newPhase['Atoms'] = newAtoms
347    newPhase['Atoms'],atCodes = GetUnique(newPhase,atCodes)
348    newPhase['Drawing'] = []
349    newPhase['ranId'] = ran.randint(0,sys.maxsize)
350    return newPhase,atCodes
351   
352def FillUnitCell(Phase):
353    Atoms = copy.deepcopy(Phase['Atoms'])
354    atomData = []
355    atCodes = []
356    SGData = Phase['General']['SGData']
357    SpnFlp = SGData.get('SpnFlp',[])
358    Amat,Bmat = cell2AB(Phase['General']['Cell'][1:7])
359    cx,ct,cs,cia = Phase['General']['AtomPtrs']
360    cm = 0
361    if Phase['General']['Type'] == 'magnetic':
362        cm = cx+4
363    for iat,atom in enumerate(Atoms):
364        XYZ = np.array(atom[cx:cx+3])
365        xyz = XYZ%1.
366        if atom[cia] == 'A':
367            Uij = atom[cia+2:cia+8]
368            result = G2spc.GenAtom(xyz,SGData,False,Uij,True)
369            for item in result:
370                if item[0][2] >= .95: item[0][2] -= 1.
371                atom[cx:cx+3] = item[0]
372                atom[cia+2:cia+8] = item[1]
373                if cm:
374                    Opr = abs(item[2])%100
375                    M = SGData['SGOps'][Opr-1][0]
376                    opNum = G2spc.GetOpNum(item[2],SGData)
377                    mom = np.inner(np.array(atom[cm:cm+3]),Bmat)
378                    atom[cm:cm+3] = np.inner(np.inner(mom,M),Amat)*nl.det(M)*SpnFlp[opNum-1]
379                atCodes.append('%d:%s'%(iat,str(item[2])))
380                atomData.append(atom[:cia+9])  #not SS stuff
381        else:
382            result = G2spc.GenAtom(xyz,SGData,False,Move=True)
383            for item in result:
384                if item[0][2] >= .95: item[0][2] -= 1.
385                atom[cx:cx+3] = item[0]
386                if cm:
387                    Opr = abs(item[1])%100
388                    M = SGData['SGOps'][Opr-1][0]
389                    opNum = G2spc.GetOpNum(item[1],SGData)
390                    mom = np.inner(np.array(atom[cm:cm+3]),Bmat)
391                    atom[cm:cm+3] = np.inner(np.inner(mom,M),Amat)*nl.det(M)*SpnFlp[opNum-1]
392#                GSASIIpath.IPyBreak()
393                atCodes.append('%d:%s'%(iat,str(item[1])))
394                atomData.append(atom[:cia+9])  #not SS stuff
395           
396    return atomData,atCodes
397       
398def GetUnique(Phase,atCodes):
399   
400    def noDuplicate(xyzA,XYZ,Amat):
401        if True in [np.allclose(np.inner(Amat,xyzA),np.inner(Amat,xyzB),atol=0.05) for xyzB in XYZ]:
402            return False
403        return True
404
405    cx,ct,cs,cia = Phase['General']['AtomPtrs']
406    cell = Phase['General']['Cell'][1:7]
407    Amat,Bmat = cell2AB(cell)
408    SGData = Phase['General']['SGData']
409    Atoms = Phase['Atoms']
410    Ind = len(Atoms)
411    newAtoms = []
412    newAtCodes = []
413    Indx = {}
414    XYZ = {}
415    for ind in range(Ind):
416        XYZ[ind] = np.array(Atoms[ind][cx:cx+3])%1.
417        Indx[ind] = True
418    for ind in range(Ind):
419        if Indx[ind]:
420            xyz = XYZ[ind]
421            for jnd in range(Ind):
422                if ind != jnd and Indx[jnd]:                       
423                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
424                    xyzs = np.array([equiv[0] for equiv in Equiv])
425                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
426    Ind = []
427    for ind in Indx:
428        if Indx[ind]:
429            newAtoms.append(Atoms[ind])
430            newAtCodes.append(atCodes[ind])
431    return newAtoms,newAtCodes
432           
433def calc_rVsq(A):
434    """Compute the square of the reciprocal lattice volume (1/V**2) from A'
435
436    """
437    G,g = A2Gmat(A)
438    rVsq = nl.det(G)
439    if rVsq < 0:
440        return 1
441    return rVsq
442   
443def calc_rV(A):
444    """Compute the reciprocal lattice volume (V*) from A
445    """
446    return np.sqrt(calc_rVsq(A))
447   
448def calc_V(A):
449    """Compute the real lattice volume (V) from A
450    """
451    return 1./calc_rV(A)
452
453def A2invcell(A):
454    """Compute reciprocal unit cell constants from A
455    returns tuple with a*,b*,c*,alpha*, beta*, gamma* (degrees)
456    """
457    G,g = A2Gmat(A)
458    return Gmat2cell(G)
459   
460def Gmat2AB(G):
461    """Computes orthogonalization matrix from reciprocal metric tensor G
462
463    :returns: tuple of two 3x3 numpy arrays (A,B)
464
465       * A for crystal to Cartesian transformations (A*x = np.inner(A,x) = X)
466       * B (= inverse of A) for Cartesian to crystal transformation (B*X = np.inner(B,X) = x)
467
468    """
469    cellstar = Gmat2cell(G)
470    g = nl.inv(G)
471    cell = Gmat2cell(g)
472    A = np.zeros(shape=(3,3))
473    # from Giacovazzo (Fundamentals 2nd Ed.) p.75
474    A[0][0] = cell[0]                # a
475    A[0][1] = cell[1]*cosd(cell[5])  # b cos(gamma)
476    A[0][2] = cell[2]*cosd(cell[4])  # c cos(beta)
477    A[1][1] = cell[1]*sind(cell[5])  # b sin(gamma)
478    A[1][2] = -cell[2]*cosd(cellstar[3])*sind(cell[4]) # - c cos(alpha*) sin(beta)
479    A[2][2] = 1./cellstar[2]         # 1/c*
480    B = nl.inv(A)
481    return A,B
482   
483def cell2AB(cell):
484    """Computes orthogonalization matrix from unit cell constants
485
486    :param tuple cell: a,b,c, alpha, beta, gamma (degrees)
487    :returns: tuple of two 3x3 numpy arrays (A,B)
488       A for crystal to Cartesian transformations A*x = np.inner(A,x) = X
489       B (= inverse of A) for Cartesian to crystal transformation B*X = np.inner(B,X) = x
490    """
491    G,g = cell2Gmat(cell) 
492    cellstar = Gmat2cell(G)
493    A = np.zeros(shape=(3,3))
494    # from Giacovazzo (Fundamentals 2nd Ed.) p.75
495    A[0][0] = cell[0]                # a
496    A[0][1] = cell[1]*cosd(cell[5])  # b cos(gamma)
497    A[0][2] = cell[2]*cosd(cell[4])  # c cos(beta)
498    A[1][1] = cell[1]*sind(cell[5])  # b sin(gamma)
499    A[1][2] = -cell[2]*cosd(cellstar[3])*sind(cell[4]) # - c cos(alpha*) sin(beta)
500    A[2][2] = 1./cellstar[2]         # 1/c*
501    B = nl.inv(A)
502    return A,B
503   
504def HKL2SpAng(H,cell,SGData):
505    """Computes spherical coords for hkls; view along 001
506
507    :param array H: arrays of hkl
508    :param tuple cell: a,b,c, alpha, beta, gamma (degrees)
509    :param dict SGData: space group dictionary
510    :returns: arrays of r,phi,psi (radius,inclination,azimuth) about 001
511    """
512    A,B = cell2AB(cell)
513    xH = np.inner(B.T,H)
514    r = np.sqrt(np.sum(xH**2,axis=0))
515    phi = acosd(xH[2]/r)
516    psi = atan2d(xH[1],xH[0])
517    phi = np.where(phi>90.,180.-phi,phi)
518#    GSASIIpath.IPyBreak()
519    return r,phi,psi
520   
521def U6toUij(U6):
522    """Fill matrix (Uij) from U6 = [U11,U22,U33,U12,U13,U23]
523    NB: there is a non numpy version in GSASIIspc: U2Uij
524
525    :param list U6: 6 terms of u11,u22,...
526    :returns:
527        Uij - numpy [3][3] array of uij
528    """
529    U = np.array([
530        [U6[0],  U6[3],  U6[4]], 
531        [U6[3],  U6[1],  U6[5]], 
532        [U6[4],  U6[5],  U6[2]]])
533    return U
534
535def UijtoU6(U):
536    """Fill vector [U11,U22,U33,U12,U13,U23] from Uij
537    NB: there is a non numpy version in GSASIIspc: Uij2U
538    """
539    U6 = np.array([U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]])
540    return U6
541
542def betaij2Uij(betaij,G):
543    """
544    Convert beta-ij to Uij tensors
545   
546    :param beta-ij - numpy array [beta-ij]
547    :param G: reciprocal metric tensor
548    :returns: Uij: numpy array [Uij]
549    """
550    ast = np.sqrt(np.diag(G))   #a*, b*, c*
551    Mast = np.multiply.outer(ast,ast)   
552    return R2pisq*UijtoU6(U6toUij(betaij)/Mast)
553   
554def Uij2betaij(Uij,G):
555    """
556    Convert Uij to beta-ij tensors -- stub for eventual completion
557   
558    :param Uij: numpy array [Uij]
559    :param G: reciprocal metric tensor
560    :returns: beta-ij - numpy array [beta-ij]
561    """
562    pass
563   
564def cell2GS(cell):
565    ''' returns Uij to betaij conversion matrix'''
566    G,g = cell2Gmat(cell)
567    GS = G
568    GS[0][1] = GS[1][0] = math.sqrt(GS[0][0]*GS[1][1])
569    GS[0][2] = GS[2][0] = math.sqrt(GS[0][0]*GS[2][2])
570    GS[1][2] = GS[2][1] = math.sqrt(GS[1][1]*GS[2][2])
571    return GS   
572   
573def Uij2Ueqv(Uij,GS,Amat):
574    ''' returns 1/3 trace of diagonalized U matrix'''
575    U = np.multiply(U6toUij(Uij),GS)
576    U = np.inner(Amat,np.inner(U,Amat).T)
577    E,R = nl.eigh(U)
578    return np.sum(E)/3.
579       
580def CosAngle(U,V,G):
581    """ calculate cos of angle between U & V in generalized coordinates
582    defined by metric tensor G
583
584    :param U: 3-vectors assume numpy arrays, can be multiple reflections as (N,3) array
585    :param V: 3-vectors assume numpy arrays, only as (3) vector
586    :param G: metric tensor for U & V defined space assume numpy array
587    :returns:
588        cos(phi)
589    """
590    u = (U.T/np.sqrt(np.sum(np.inner(U,G)*U,axis=1))).T
591    v = V/np.sqrt(np.inner(V,np.inner(G,V)))
592    cosP = np.inner(u,np.inner(G,v))
593    return cosP
594   
595def CosSinAngle(U,V,G):
596    """ calculate sin & cos of angle between U & V in generalized coordinates
597    defined by metric tensor G
598
599    :param U: 3-vectors assume numpy arrays
600    :param V: 3-vectors assume numpy arrays
601    :param G: metric tensor for U & V defined space assume numpy array
602    :returns:
603        cos(phi) & sin(phi)
604    """
605    u = U/np.sqrt(np.inner(U,np.inner(G,U)))
606    v = V/np.sqrt(np.inner(V,np.inner(G,V)))
607    cosP = np.inner(u,np.inner(G,v))
608    sinP = np.sqrt(max(0.0,1.0-cosP**2))
609    return cosP,sinP
610   
611def criticalEllipse(prob):
612    """
613    Calculate critical values for probability ellipsoids from probability
614    """
615    if not ( 0.01 <= prob < 1.0):
616        return 1.54 
617    coeff = np.array([6.44988E-09,4.16479E-07,1.11172E-05,1.58767E-04,0.00130554,
618        0.00604091,0.0114921,-0.040301,-0.6337203,1.311582])
619    llpr = math.log(-math.log(prob))
620    return np.polyval(coeff,llpr)
621   
622def CellBlock(nCells):
623    """
624    Generate block of unit cells n*n*n on a side; [0,0,0] centered, n = 2*nCells+1
625    currently only works for nCells = 0 or 1 (not >1)
626    """
627    if nCells:
628        N = 2*nCells+1
629        N2 = N*N
630        N3 = N*N*N
631        cellArray = []
632        A = np.array(range(N3))
633        cellGen = np.array([A//N2-1,A//N%N-1,A%N-1]).T
634        for cell in cellGen:
635            cellArray.append(cell)
636        return cellArray
637    else:
638        return [0,0,0]
639       
640def CellAbsorption(ElList,Volume):
641    '''Compute unit cell absorption
642
643    :param dict ElList: dictionary of element contents including mu and
644      number of atoms be cell
645    :param float Volume: unit cell volume
646    :returns: mu-total/Volume
647    '''
648    muT = 0
649    for El in ElList:
650        muT += ElList[El]['mu']*ElList[El]['FormulaNo']
651    return muT/Volume
652   
653#Permutations and Combinations
654# Four routines: combinations,uniqueCombinations, selections & permutations
655#These taken from Python Cookbook, 2nd Edition. 19.15 p724-726
656#   
657def _combinators(_handle, items, n):
658    """ factored-out common structure of all following combinators """
659    if n==0:
660        yield [ ]
661        return
662    for i, item in enumerate(items):
663        this_one = [ item ]
664        for cc in _combinators(_handle, _handle(items, i), n-1):
665            yield this_one + cc
666def combinations(items, n):
667    """ take n distinct items, order matters """
668    def skipIthItem(items, i):
669        return items[:i] + items[i+1:]
670    return _combinators(skipIthItem, items, n)
671def uniqueCombinations(items, n):
672    """ take n distinct items, order is irrelevant """
673    def afterIthItem(items, i):
674        return items[i+1:]
675    return _combinators(afterIthItem, items, n)
676def selections(items, n):
677    """ take n (not necessarily distinct) items, order matters """
678    def keepAllItems(items, i):
679        return items
680    return _combinators(keepAllItems, items, n)
681def permutations(items):
682    """ take all items, order matters """
683    return combinations(items, len(items))
684
685#reflection generation routines
686#for these: H = [h,k,l]; A is as used in calc_rDsq; G - inv metric tensor, g - metric tensor;
687#           cell - a,b,c,alp,bet,gam in A & deg
688   
689def Pos2dsp(Inst,pos):
690    ''' convert powder pattern position (2-theta or TOF, musec) to d-spacing
691    '''
692    if 'C' in Inst['Type'][0] or 'PKS' in Inst['Type'][0]:
693        wave = G2mth.getWave(Inst)
694        return wave/(2.0*sind((pos-Inst.get('Zero',[0,0])[1])/2.0))
695    else:   #'T'OF - ignore difB
696        return TOF2dsp(Inst,pos)
697       
698def TOF2dsp(Inst,Pos):
699    ''' convert powder pattern TOF, musec to d-spacing by successive approximation
700    Pos can be numpy array
701    '''
702    def func(d,pos,Inst):       
703        return (pos-Inst['difA'][1]*d**2-Inst['Zero'][1]-Inst['difB'][1]/d)/Inst['difC'][1]
704    dsp0 = np.ones_like(Pos)
705    N = 0
706    while True:      #successive approximations
707        dsp = func(dsp0,Pos,Inst)
708        if np.allclose(dsp,dsp0,atol=0.000001):
709            return dsp
710        dsp0 = dsp
711        N += 1
712        if N > 10:
713            return dsp
714   
715def Dsp2pos(Inst,dsp):
716    ''' convert d-spacing to powder pattern position (2-theta or TOF, musec)
717    '''
718    if 'C' in Inst['Type'][0] or 'PKS' in Inst['Type'][0]:
719        wave = G2mth.getWave(Inst)
720        val = min(0.995,wave/(2.*dsp))  #set max at 168deg
721        pos = 2.0*asind(val)+Inst.get('Zero',[0,0])[1]             
722    else:   #'T'OF
723        pos = Inst['difC'][1]*dsp+Inst['Zero'][1]+Inst['difA'][1]*dsp**2+Inst.get('difB',[0,0,False])[1]/dsp
724    return pos
725   
726def getPeakPos(dataType,parmdict,dsp):
727    ''' convert d-spacing to powder pattern position (2-theta or TOF, musec)
728    '''
729    if 'C' in dataType:
730        pos = 2.0*asind(parmdict['Lam']/(2.*dsp))+parmdict['Zero']
731    else:   #'T'OF
732        pos = parmdict['difC']*dsp+parmdict['difA']*dsp**2+parmdict['difB']/dsp+parmdict['Zero']
733    return pos
734                   
735def calc_rDsq(H,A):
736    'needs doc string'
737    rdsq = H[0]*H[0]*A[0]+H[1]*H[1]*A[1]+H[2]*H[2]*A[2]+H[0]*H[1]*A[3]+H[0]*H[2]*A[4]+H[1]*H[2]*A[5]
738    return rdsq
739   
740def calc_rDsq2(H,G):
741    'needs doc string'
742    return np.inner(H,np.inner(G,H))
743   
744def calc_rDsqSS(H,A,vec):
745    'needs doc string'
746    rdsq = calc_rDsq(H[:3]+(H[3]*vec).T,A)
747    return rdsq
748       
749def calc_rDsqZ(H,A,Z,tth,lam):
750    'needs doc string'
751    rdsq = calc_rDsq(H,A)+Z*sind(tth)*2.0*rpd/lam**2
752    return rdsq
753       
754def calc_rDsqZSS(H,A,vec,Z,tth,lam):
755    'needs doc string'
756    rdsq = calc_rDsq(H[:3]+(H[3][:,np.newaxis]*vec).T,A)+Z*sind(tth)*2.0*rpd/lam**2
757    return rdsq
758       
759def calc_rDsqT(H,A,Z,tof,difC):
760    'needs doc string'
761    rdsq = calc_rDsq(H,A)+Z/difC
762    return rdsq
763       
764def calc_rDsqTSS(H,A,vec,Z,tof,difC):
765    'needs doc string'
766    rdsq = calc_rDsq(H[:3]+(H[3][:,np.newaxis]*vec).T,A)+Z/difC
767    return rdsq
768   
769def PlaneIntercepts(Amat,Bmat,H,phase,stack):
770    ''' find unit cell intercepts for a stack of hkl planes
771    '''
772    Steps = range(-1,2,2)
773    if stack:
774        Steps = range(-10,10,1)
775    Stack = []
776    Ux = np.array([[0,0],[1,0],[1,1],[0,1]])
777    for step in Steps:
778        HX = []
779        for i in [0,1,2]:
780            if H[i]:
781               h,k,l = [(i+1)%3,(i+2)%3,(i+3)%3]
782               for j in [0,1,2,3]:
783                    hx = [0,0,0]
784                    intcpt = (phase/360.+step-H[h]*Ux[j,0]-H[k]*Ux[j,1])/H[l]
785                    if 0. <= intcpt <= 1.:                       
786                        hx[h] = Ux[j,0]
787                        hx[k] = Ux[j,1]
788                        hx[l] = intcpt
789                        HX.append(hx)
790        if len(HX)> 2:
791            HX = np.array(HX)
792            DX = np.inner(HX-HX[0],Amat)
793            D = np.sqrt(np.sum(DX**2,axis=1))
794            Dsort = np.argsort(D)
795            HX = HX[Dsort]
796            DX = DX[Dsort]
797            D = D[Dsort]
798            DX[1:,:] = DX[1:,:]/D[1:,nxs]
799            A = 2.*np.ones(HX.shape[0])
800            A[1:] = [np.dot(DX[1],dx) for dx in DX[1:]]
801            HX = HX[np.argsort(A)]
802#            GSASIIpath.IPyBreak()
803            Stack.append(HX)
804    return Stack
805       
806def MaxIndex(dmin,A):
807    'needs doc string'
808    Hmax = [0,0,0]
809    try:
810        cell = A2cell(A)
811    except:
812        cell = [1,1,1,90,90,90]
813    for i in range(3):
814        Hmax[i] = int(round(cell[i]/dmin))
815    return Hmax
816   
817def transposeHKLF(transMat,Super,refList):
818    ''' Apply transformation matrix to hkl(m)
819    param: transmat: 3x3 or 4x4 array
820    param: Super: 0 or 1 for extra index
821    param: refList list of h,k,l,....
822    return: newRefs transformed list of h',k',l',,,
823    return: badRefs list of noninteger h',k',l'...
824    '''
825    newRefs = np.copy(refList)
826    badRefs = []
827    for H in newRefs:
828        newH = np.inner(transMat,H[:3+Super])
829        H[:3+Super] = np.rint(newH)
830        if not np.allclose(newH,H[:3+Super],atol=0.01):
831            badRefs.append(newH)
832    return newRefs,badRefs
833   
834def sortHKLd(HKLd,ifreverse,ifdup,ifSS=False):
835    '''sort reflection list on d-spacing; can sort in either order
836
837    :param HKLd: a list of [h,k,l,d,...];
838    :param ifreverse: True for largest d first
839    :param ifdup: True if duplicate d-spacings allowed
840    :return: sorted reflection list
841    '''
842    T = []
843    N = 3
844    if ifSS:
845        N = 4
846    for i,H in enumerate(HKLd):
847        if ifdup:
848            T.append((H[N],i))
849        else:
850            T.append(H[N])           
851    D = dict(zip(T,HKLd))
852    T.sort()
853    if ifreverse:
854        T.reverse()
855    X = []
856    okey = ''
857    for key in T: 
858        if key != okey: X.append(D[key])    #remove duplicate d-spacings
859        okey = key
860    return X
861   
862def SwapIndx(Axis,H):
863    'needs doc string'
864    if Axis in [1,-1]:
865        return H
866    elif Axis in [2,-3]:
867        return [H[1],H[2],H[0]]
868    else:
869        return [H[2],H[0],H[1]]
870       
871def Rh2Hx(Rh):
872    'needs doc string'
873    Hx = [0,0,0]
874    Hx[0] = Rh[0]-Rh[1]
875    Hx[1] = Rh[1]-Rh[2]
876    Hx[2] = np.sum(Rh)
877    return Hx
878   
879def Hx2Rh(Hx):
880    'needs doc string'
881    Rh = [0,0,0]
882    itk = -Hx[0]+Hx[1]+Hx[2]
883    if itk%3 != 0:
884        return 0        #error - not rhombohedral reflection
885    else:
886        Rh[1] = itk//3
887        Rh[0] = Rh[1]+Hx[0]
888        Rh[2] = Rh[1]-Hx[1]
889        if Rh[0] < 0:
890            for i in range(3):
891                Rh[i] = -Rh[i]
892        return Rh
893       
894def CentCheck(Cent,H):
895    'needs doc string'
896    h,k,l = H
897    if Cent == 'A' and (k+l)%2:
898        return False
899    elif Cent == 'B' and (h+l)%2:
900        return False
901    elif Cent == 'C' and (h+k)%2:
902        return False
903    elif Cent == 'I' and (h+k+l)%2:
904        return False
905    elif Cent == 'F' and ((h+k)%2 or (h+l)%2 or (k+l)%2):
906        return False
907    elif Cent == 'R' and (-h+k+l)%3:
908        return False
909    else:
910        return True
911                                   
912def GetBraviasNum(center,system):
913    """Determine the Bravais lattice number, as used in GenHBravais
914   
915    :param center: one of: 'P', 'C', 'I', 'F', 'R' (see SGLatt from GSASIIspc.SpcGroup)
916    :param system: one of 'cubic', 'hexagonal', 'tetragonal', 'orthorhombic', 'trigonal' (for R)
917      'monoclinic', 'triclinic' (see SGSys from GSASIIspc.SpcGroup)
918    :return: a number between 0 and 13
919      or throws a ValueError exception if the combination of center, system is not found (i.e. non-standard)
920
921    """
922    if center.upper() == 'F' and system.lower() == 'cubic':
923        return 0
924    elif center.upper() == 'I' and system.lower() == 'cubic':
925        return 1
926    elif center.upper() == 'P' and system.lower() == 'cubic':
927        return 2
928    elif center.upper() == 'R' and system.lower() == 'trigonal':
929        return 3
930    elif center.upper() == 'P' and system.lower() == 'hexagonal':
931        return 4
932    elif center.upper() == 'I' and system.lower() == 'tetragonal':
933        return 5
934    elif center.upper() == 'P' and system.lower() == 'tetragonal':
935        return 6
936    elif center.upper() == 'F' and system.lower() == 'orthorhombic':
937        return 7
938    elif center.upper() == 'I' and system.lower() == 'orthorhombic':
939        return 8
940    elif center.upper() == 'C' and system.lower() == 'orthorhombic':
941        return 9
942    elif center.upper() == 'P' and system.lower() == 'orthorhombic':
943        return 10
944    elif center.upper() == 'C' and system.lower() == 'monoclinic':
945        return 11
946    elif center.upper() == 'P' and system.lower() == 'monoclinic':
947        return 12
948    elif center.upper() == 'P' and system.lower() == 'triclinic':
949        return 13
950    raise ValueError('non-standard Bravais lattice center=%s, cell=%s' % (center,system))
951
952def GenHBravais(dmin,Bravais,A):
953    """Generate the positionally unique powder diffraction reflections
954     
955    :param dmin: minimum d-spacing in A
956    :param Bravais: lattice type (see GetBraviasNum). Bravais is one of:
957   
958            * 0 F cubic
959            * 1 I cubic
960            * 2 P cubic
961            * 3 R hexagonal (trigonal not rhombohedral)
962            * 4 P hexagonal
963            * 5 I tetragonal
964            * 6 P tetragonal
965            * 7 F orthorhombic
966            * 8 I orthorhombic
967            * 9 C orthorhombic
968            * 10 P orthorhombic
969            * 11 C monoclinic
970            * 12 P monoclinic
971            * 13 P triclinic
972           
973    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
974    :return: HKL unique d list of [h,k,l,d,-1] sorted with largest d first
975           
976    """
977    if Bravais in [9,11]:
978        Cent = 'C'
979    elif Bravais in [1,5,8]:
980        Cent = 'I'
981    elif Bravais in [0,7]:
982        Cent = 'F'
983    elif Bravais in [3]:
984        Cent = 'R'
985    else:
986        Cent = 'P'
987    Hmax = MaxIndex(dmin,A)
988    dminsq = 1./(dmin**2)
989    HKL = []
990    if Bravais == 13:                       #triclinic
991        for l in range(-Hmax[2],Hmax[2]+1):
992            for k in range(-Hmax[1],Hmax[1]+1):
993                hmin = 0
994                if (k < 0): hmin = 1
995                if (k ==0 and l < 0): hmin = 1
996                for h in range(hmin,Hmax[0]+1):
997                    H=[h,k,l]
998                    rdsq = calc_rDsq(H,A)
999                    if 0 < rdsq <= dminsq:
1000                        HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1001    elif Bravais in [11,12]:                #monoclinic - b unique
1002        Hmax = SwapIndx(2,Hmax)
1003        for h in range(Hmax[0]+1):
1004            for k in range(-Hmax[1],Hmax[1]+1):
1005                lmin = 0
1006                if k < 0:lmin = 1
1007                for l in range(lmin,Hmax[2]+1):
1008                    [h,k,l] = SwapIndx(-2,[h,k,l])
1009                    H = []
1010                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1011                    if H:
1012                        rdsq = calc_rDsq(H,A)
1013                        if 0 < rdsq <= dminsq:
1014                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1015                    [h,k,l] = SwapIndx(2,[h,k,l])
1016    elif Bravais in [7,8,9,10]:            #orthorhombic
1017        for h in range(Hmax[0]+1):
1018            for k in range(Hmax[1]+1):
1019                for l in range(Hmax[2]+1):
1020                    H = []
1021                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1022                    if H:
1023                        rdsq = calc_rDsq(H,A)
1024                        if 0 < rdsq <= dminsq:
1025                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1026    elif Bravais in [5,6]:                  #tetragonal
1027        for l in range(Hmax[2]+1):
1028            for k in range(Hmax[1]+1):
1029                for h in range(k,Hmax[0]+1):
1030                    H = []
1031                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1032                    if H:
1033                        rdsq = calc_rDsq(H,A)
1034                        if 0 < rdsq <= dminsq:
1035                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1036    elif Bravais in [3,4]:
1037        lmin = 0
1038        if Bravais == 3: lmin = -Hmax[2]                  #hexagonal/trigonal
1039        for l in range(lmin,Hmax[2]+1):
1040            for k in range(Hmax[1]+1):
1041                hmin = k
1042                if l < 0: hmin += 1
1043                for h in range(hmin,Hmax[0]+1):
1044                    H = []
1045                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1046                    if H:
1047                        rdsq = calc_rDsq(H,A)
1048                        if 0 < rdsq <= dminsq:
1049                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1050
1051    else:                                   #cubic
1052        for l in range(Hmax[2]+1):
1053            for k in range(l,Hmax[1]+1):
1054                for h in range(k,Hmax[0]+1):
1055                    H = []
1056                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
1057                    if H:
1058                        rdsq = calc_rDsq(H,A)
1059                        if 0 < rdsq <= dminsq:
1060                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
1061    return sortHKLd(HKL,True,False)
1062   
1063def getHKLmax(dmin,SGData,A):
1064    'finds maximum allowed hkl for given A within dmin'
1065    SGLaue = SGData['SGLaue']
1066    if SGLaue in ['3R','3mR']:        #Rhombohedral axes
1067        Hmax = [0,0,0]
1068        cell = A2cell(A)
1069        aHx = cell[0]*math.sqrt(2.0*(1.0-cosd(cell[3])))
1070        cHx = cell[0]*math.sqrt(3.0*(1.0+2.0*cosd(cell[3])))
1071        Hmax[0] = Hmax[1] = int(round(aHx/dmin))
1072        Hmax[2] = int(round(cHx/dmin))
1073        #print Hmax,aHx,cHx
1074    else:                           # all others
1075        Hmax = MaxIndex(dmin,A)
1076    return Hmax
1077   
1078def GenHLaue(dmin,SGData,A):
1079    """Generate the crystallographically unique powder diffraction reflections
1080    for a lattice and Bravais type
1081   
1082    :param dmin: minimum d-spacing
1083    :param SGData: space group dictionary with at least
1084   
1085        * 'SGLaue': Laue group symbol: one of '-1','2/m','mmm','4/m','6/m','4/mmm','6/mmm', '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
1086        * 'SGLatt': lattice centering: one of 'P','A','B','C','I','F'
1087        * 'SGUniq': code for unique monoclinic axis one of 'a','b','c' (only if 'SGLaue' is '2/m') otherwise an empty string
1088       
1089    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
1090    :return: HKL = list of [h,k,l,d] sorted with largest d first and is unique
1091            part of reciprocal space ignoring anomalous dispersion
1092           
1093    """
1094    import math
1095    SGLaue = SGData['SGLaue']
1096    SGLatt = SGData['SGLatt']
1097    SGUniq = SGData['SGUniq']
1098    #finds maximum allowed hkl for given A within dmin
1099    Hmax = getHKLmax(dmin,SGData,A)
1100       
1101    dminsq = 1./(dmin**2)
1102    HKL = []
1103    if SGLaue == '-1':                       #triclinic
1104        for l in range(-Hmax[2],Hmax[2]+1):
1105            for k in range(-Hmax[1],Hmax[1]+1):
1106                hmin = 0
1107                if (k < 0) or (k ==0 and l < 0): hmin = 1
1108                for h in range(hmin,Hmax[0]+1):
1109                    H = []
1110                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1111                    if H:
1112                        rdsq = calc_rDsq(H,A)
1113                        if 0 < rdsq <= dminsq:
1114                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1115    elif SGLaue == '2/m':                #monoclinic
1116        axisnum = 1 + ['a','b','c'].index(SGUniq)
1117        Hmax = SwapIndx(axisnum,Hmax)
1118        for h in range(Hmax[0]+1):
1119            for k in range(-Hmax[1],Hmax[1]+1):
1120                lmin = 0
1121                if k < 0:lmin = 1
1122                for l in range(lmin,Hmax[2]+1):
1123                    [h,k,l] = SwapIndx(-axisnum,[h,k,l])
1124                    H = []
1125                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1126                    if H:
1127                        rdsq = calc_rDsq(H,A)
1128                        if 0 < rdsq <= dminsq:
1129                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1130                    [h,k,l] = SwapIndx(axisnum,[h,k,l])
1131    elif SGLaue in ['mmm','4/m','6/m']:            #orthorhombic
1132        for l in range(Hmax[2]+1):
1133            for h in range(Hmax[0]+1):
1134                kmin = 1
1135                if SGLaue == 'mmm' or h ==0: kmin = 0
1136                for k in range(kmin,Hmax[1]+1):
1137                    H = []
1138                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1139                    if H:
1140                        rdsq = calc_rDsq(H,A)
1141                        if 0 < rdsq <= dminsq:
1142                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1143    elif SGLaue in ['4/mmm','6/mmm']:                  #tetragonal & hexagonal
1144        for l in range(Hmax[2]+1):
1145            for h in range(Hmax[0]+1):
1146                for k in range(h+1):
1147                    H = []
1148                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1149                    if H:
1150                        rdsq = calc_rDsq(H,A)
1151                        if 0 < rdsq <= dminsq:
1152                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1153    elif SGLaue in ['3m1','31m','3','3R','3mR']:                  #trigonals
1154        for l in range(-Hmax[2],Hmax[2]+1):
1155            hmin = 0
1156            if l < 0: hmin = 1
1157            for h in range(hmin,Hmax[0]+1):
1158                if SGLaue in ['3R','3']:
1159                    kmax = h
1160                    kmin = -int((h-1.)/2.)
1161                else:
1162                    kmin = 0
1163                    kmax = h
1164                    if SGLaue in ['3m1','3mR'] and l < 0: kmax = h-1
1165                    if SGLaue == '31m' and l < 0: kmin = 1
1166                for k in range(kmin,kmax+1):
1167                    H = []
1168                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1169                    if SGLaue in ['3R','3mR']:
1170                        H = Hx2Rh(H)
1171                    if H:
1172                        rdsq = calc_rDsq(H,A)
1173                        if 0 < rdsq <= dminsq:
1174                            HKL.append([H[0],H[1],H[2],1./math.sqrt(rdsq)])
1175    else:                                   #cubic
1176        for h in range(Hmax[0]+1):
1177            for k in range(h+1):
1178                lmin = 0
1179                lmax = k
1180                if SGLaue =='m3':
1181                    lmax = h-1
1182                    if h == k: lmax += 1
1183                for l in range(lmin,lmax+1):
1184                    H = []
1185                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
1186                    if H:
1187                        rdsq = calc_rDsq(H,A)
1188                        if 0 < rdsq <= dminsq:
1189                            HKL.append([h,k,l,1./math.sqrt(rdsq)])
1190    return sortHKLd(HKL,True,True)
1191   
1192def GenPfHKLs(nMax,SGData,A):   
1193    """Generate the unique pole figure reflections for a lattice and Bravais type.
1194    Min d-spacing=1.0A & no more than nMax returned
1195   
1196    :param nMax: maximum number of hkls returned
1197    :param SGData: space group dictionary with at least
1198   
1199        * 'SGLaue': Laue group symbol: one of '-1','2/m','mmm','4/m','6/m','4/mmm','6/mmm', '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
1200        * 'SGLatt': lattice centering: one of 'P','A','B','C','I','F'
1201        * 'SGUniq': code for unique monoclinic axis one of 'a','b','c' (only if 'SGLaue' is '2/m') otherwise an empty string
1202       
1203    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
1204    :return: HKL = list of 'h k l' strings sorted with largest d first; no duplicate zones
1205           
1206    """
1207    HKL = np.array(GenHLaue(1.0,SGData,A)).T[:3].T     #strip d-spacings
1208    N = min(nMax,len(HKL))
1209    return ['%d %d %d'%(h[0],h[1],h[2]) for h in HKL[:N]]       
1210
1211def GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A):
1212    'needs a doc string'
1213    HKLs = []
1214    vec = np.array(Vec)
1215    vstar = np.sqrt(calc_rDsq(vec,A))     #find extra needed for -n SS reflections
1216    dvec = 1./(maxH*vstar+1./dmin)
1217    HKL = GenHLaue(dvec,SGData,A)       
1218    SSdH = [vec*h for h in range(-maxH,maxH+1)]
1219    SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
1220    for h,k,l,d in HKL:
1221        ext = G2spc.GenHKLf([h,k,l],SGData)[0]  #h,k,l must be integral values here
1222        if not ext and d >= dmin:
1223            HKLs.append([h,k,l,0,d])
1224        for dH in SSdH:
1225            if dH:
1226                DH = SSdH[dH]
1227                H = [h+DH[0],k+DH[1],l+DH[2]]
1228                d = 1./np.sqrt(calc_rDsq(H,A))
1229                if d >= dmin:
1230                    HKLM = np.array([h,k,l,dH])
1231                    if G2spc.checkSSLaue([h,k,l,dH],SGData,SSGData) and G2spc.checkSSextc(HKLM,SSGData):
1232                        HKLs.append([h,k,l,dH,d])   
1233    return HKLs
1234   
1235def LaueUnique2(SGData,refList):
1236    ''' Impose Laue symmetry on hkl
1237   
1238    :param SGData: space group data from 'P '+Laue
1239    :param HKLF: np.array([[h,k,l,...]]) reflection set to be converted
1240   
1241    :return: HKLF new reflection array with imposed Laue symmetry
1242    '''
1243    for ref in refList:
1244        H = ref[:3]
1245        Uniq = G2spc.GenHKLf(H,SGData)[2]
1246        Uniq = G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Uniq,2),1),0)
1247        ref[:3] = Uniq[-1]
1248    return refList
1249   
1250def LaueUnique(Laue,HKLF):
1251    ''' Impose Laue symmetry on hkl
1252   
1253    :param str Laue: Laue symbol, as below
1254   
1255      centrosymmetric Laue groups::
1256       
1257            ['-1','2/m','112/m','2/m11','mmm','-42m','-4m2','4/mmm','-3',
1258            '-31m','-3m1','6/m','6/mmm','m3','m3m']
1259     
1260      noncentrosymmetric Laue groups::
1261     
1262           ['1','2','211','112','m','m11','11m','222','mm2','m2m','2mm',
1263           '4','-4','422','4mm','3','312','321','31m','3m1','6','-6',
1264           '622','6mm','-62m','-6m2','23','432','-43m']
1265     
1266    :param HKLF: np.array([[h,k,l,...]]) reflection set to be converted
1267   
1268    :returns: HKLF new reflection array with imposed Laue symmetry
1269    '''
1270   
1271    HKLFT = HKLF.T
1272    mat41 = np.array([[0,1,0],[-1,0,0],[0,0,1]])    #hkl -> k,-h,l
1273    mat43 = np.array([[0,-1,0],[1,0,0],[0,0,1]])    #hkl -> -k,h,l
1274    mat4bar = np.array([[0,-1,0],[1,0,0],[0,0,-1]]) #hkl -> k,-h,-l
1275    mat31 = np.array([[-1,-1,0],[1,0,0],[0,0,1]])   #hkl -> ihl = -h-k,h,l
1276    mat32 = np.array([[0,1,0],[-1,-1,0],[0,0,1]])   #hkl -> kil = k,-h-k,l
1277    matd3 = np.array([[0,1,0],[0,0,1],[1,0,0]])     #hkl -> k,l,h
1278    matd3q = np.array([[0,0,-1],[-1,0,0],[0,1,0]])  #hkl -> -l,-h,k
1279    matd3t = np.array([[0,0,-1],[1,0,0],[0,-1,0]])  #hkl -> -l,h,-k
1280    mat6 = np.array([[1,1,0],[-1,0,0],[0,0,1]])     #hkl -> h+k,-h,l really 65
1281    matdm = np.array([[0,1,0],[1,0,0],[0,0,1]])     #hkl -> k,h,l
1282    matdmp = np.array([[-1,-1,0],[0,1,0],[0,0,1]])  #hkl -> -h-k,k,l
1283    matkm = np.array([[-1,0,0],[1,1,0],[0,0,1]])    #hkl -> -h,h+k,l
1284    matd2 = np.array([[0,1,0],[1,0,0],[0,0,-1]])    #hkl -> k,h,-l
1285    matdm3 = np.array([[1,0,0],[0,0,1],[0,1,0]])    #hkl -> h,l,k
1286    mat2d43 = np.array([[0,1,0],[1,0,0],[0,0,1]])   #hkl -> k,-h,l
1287    matk2 = np.array([[-1,0,0],[1,1,0],[0,0,-1]])   #hkl -> -h,-i,-l
1288    #triclinic
1289    if Laue == '1': #ok
1290        pass
1291    elif Laue == '-1':  #ok
1292        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1293        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1294        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1295    #monoclinic
1296    #noncentrosymmetric - all ok
1297    elif Laue == '2': 
1298        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1299        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1300    elif Laue == '1 1 2':
1301        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1302        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1303    elif Laue == '2 1 1':   
1304        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1305        HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1306    elif Laue == 'm':
1307        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1308    elif Laue == 'm 1 1':
1309        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1310    elif Laue == '1 1 m':
1311        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1312    #centrosymmetric - all ok
1313    elif Laue == '2/m 1 1':       
1314        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1315        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1316        HKLFT[:3] = np.where((HKLFT[2]*HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1317    elif Laue == '2/m':
1318        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])
1319        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1320        HKLFT[:3] = np.where((HKLFT[0]*HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1321    elif Laue == '1 1 2/m':
1322        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1323        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1324        HKLFT[:3] = np.where((HKLFT[1]*HKLFT[2]==0)&(HKLFT[0]<0),HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1325    #orthorhombic
1326    #noncentrosymmetric - all OK
1327    elif Laue == '2 2 2':
1328        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1329        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1330        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1331        HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1332    elif Laue == 'm m 2':
1333        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1334        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1335    elif Laue == '2 m m': 
1336        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1337        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1338    elif Laue == 'm 2 m':
1339        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1340        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1341    #centrosymmetric - all ok
1342    elif Laue == 'm m m':
1343        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1344        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1345        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1346    #tetragonal
1347    #noncentrosymmetric - all ok
1348    elif Laue == '4':
1349        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1350        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1351        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1352    elif Laue == '-4': 
1353        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1354        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1355        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1356        HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1357        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1358    elif Laue == '4 2 2':
1359        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1360        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1361        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1362        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1363        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])   #in lieu od 2-fold
1364    elif Laue == '4 m m':
1365        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1366        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1367        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1368        HKLFT[:3] = np.where(HKLFT[0]<HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1369    elif Laue == '-4 2 m':
1370        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1371        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1372        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1373        HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1374        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1375        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1376        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1377    elif Laue == '-4 m 2':
1378        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1379        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1380        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<=0),np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1381        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]<0),HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1382        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]==0),np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1383        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3]) 
1384        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]>HKLFT[1]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1385    #centrosymmetric - all ok
1386    elif Laue == '4/m':
1387        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1388        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1389        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1390        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]>0),np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1391    elif Laue == '4/m m m':
1392        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1393        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1394        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])       
1395        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat41[nxs,:,:])).T,HKLFT[:3])
1396        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1397    #trigonal - all hex cell
1398    #noncentrosymmetric - all ok
1399    elif Laue == '3':
1400        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1401        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1402        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1403    elif Laue == '3 1 2':
1404        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1405        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1406        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1407        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1408        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1409    elif Laue == '3 2 1':
1410        HKLFT[:3] = np.where(HKLFT[0]<=-2*HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1411        HKLFT[:3] = np.where(HKLFT[1]<-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1412        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1413        HKLFT[:3] = np.where((HKLFT[2]>0)&(HKLFT[1]==HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1414        HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T
1415        HKLFT[:3] = np.where((HKLFT[0]!=0)&(HKLFT[2]>0)&(HKLFT[0]==-2*HKLFT[1]),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1416    elif Laue == '3 1 m':
1417        HKLFT[:3] = np.where(HKLFT[0]>=HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1418        HKLFT[:3] = np.where(2*HKLFT[1]<-HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1419        HKLFT[:3] = np.where(HKLFT[1]>-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdmp[nxs,:,:])).T,HKLFT[:3])
1420        HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T
1421    elif Laue == '3 m 1':
1422        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1423        HKLFT[:3] = np.where((HKLFT[1]+HKLFT[0])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1424        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matkm[nxs,:,:])).T,HKLFT[:3])
1425    #centrosymmetric
1426    elif Laue == '-3':  #ok
1427        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1428        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1429        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1430        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1431        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]<0),-np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1432        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],-mat31[nxs,:,:])).T,HKLFT[:3])   
1433    elif Laue == '-3 m 1':  #ok
1434        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1435        HKLFT[:3] = np.where((HKLFT[1]+HKLFT[0])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1436        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matkm[nxs,:,:])).T,HKLFT[:3])
1437        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1438        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1439    elif Laue == '-3 1 m':  #ok
1440        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,-1,-1])[:,nxs],HKLFT[:3])
1441        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1442        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1443        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1444        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],-mat31[nxs,:,:])).T,HKLFT[:3])   
1445        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1446    #hexagonal
1447    #noncentrosymmetric
1448    elif Laue == '6':   #ok
1449        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1450        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1451        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1452        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1453    elif Laue == '-6':  #ok
1454        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1455        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1456        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1457        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1458    elif Laue == '6 2 2':   #ok
1459        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1460        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1461        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1462        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1463        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1464        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[0]>HKLFT[1]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1465    elif Laue == '6 m m':   #ok
1466        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1467        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1468        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1469        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1470        HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1471    elif Laue == '-6 m 2':  #ok
1472        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1473        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1474        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1475        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat31[nxs,:,:])).T,HKLFT[:3])
1476        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],matk2[nxs,:,:])).T,HKLFT[:3])
1477        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1478    elif Laue == '-6 2 m':  #ok
1479        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1480        HKLFT[:3] = np.where(HKLFT[0]<=-2*HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1481        HKLFT[:3] = np.where(HKLFT[1]<-2*HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1482        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1483        HKLFT[:3] = np.where((HKLFT[2]>0)&(HKLFT[1]==HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1484        HKLFT[:3] = np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T
1485        HKLFT[:3] = np.where(HKLFT[2]<0,np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1486        HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1487    #centrosymmetric
1488    elif Laue == '6/m': #ok
1489        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1490        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1491        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1492        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1493        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1494    elif Laue == '6/m m m': #ok
1495        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1496        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1497        HKLFT[:3] = np.where((HKLFT[0]+HKLFT[1])<0,np.squeeze(np.inner(HKLF[:,:3],mat32[nxs,:,:])).T,HKLFT[:3])
1498        HKLFT[:3] = np.where(HKLFT[0]<0,np.squeeze(np.inner(HKLF[:,:3],mat6[nxs,:,:])).T,HKLFT[:3])
1499        HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm.T[nxs,:,:])).T,HKLFT[:3])
1500    #cubic - all ok
1501    #noncentrosymmetric -
1502    elif Laue == '2 3': 
1503        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1504        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1505        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1506        HKLFT[:3] = np.where((HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1507        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1508        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1509        HKLFT[:3] = np.where((HKLFT[2]<0)&((HKLFT[0]>-HKLFT[2])|(HKLFT[1]>-HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3t[nxs,:,:])).T,HKLFT[:3])
1510        HKLFT[:3] = np.where((HKLFT[2]<0)&((HKLFT[0]>-HKLFT[2])|(HKLFT[1]>=-HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3t[nxs,:,:])).T,HKLFT[:3])
1511        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3])       
1512    elif Laue == '4 3 2':   
1513        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,-1,-1])[:,nxs],HKLFT[:3])
1514        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])
1515        HKLFT[:3] = np.where(HKLFT[1]<0,np.squeeze(np.inner(HKLF[:,:3],mat43[nxs,:,:])).T,HKLFT[:3])
1516        HKLFT[:3] = np.where((HKLFT[2]==0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matd2[nxs,:,:])).T,HKLFT[:3])
1517        HKLFT[:3] = np.where(HKLFT[0]==0,np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])   #in lieu od 2-fold
1518        HKLFT[:3] = np.where((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1519        HKLFT[:3] = np.where((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1520        HKLFT[:3] = np.where(HKLFT[1]==0,np.squeeze(np.inner(HKLF[:,:3],mat2d43[nxs,:,:])).T,HKLFT[:3])
1521    elif Laue == '-4 3 m': 
1522        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1523        HKLFT[:3] = np.where(HKLFT[0]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1524        HKLFT[:3] = np.where(HKLFT[0]<=0,HKLFT[:3]*np.array([-1,-1,1])[:,nxs],HKLFT[:3])     
1525        HKLFT[:3] = np.where(HKLFT[1]<=0,np.squeeze(np.inner(HKLF[:,:3],mat4bar[nxs,:,:])).T,HKLFT[:3])
1526        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[1]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1527        HKLFT[:3] = np.where(HKLFT[1]<HKLFT[0],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1528        HKLFT[:3] = np.where((HKLFT[0]==0)&(HKLFT[2]<0),HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])
1529        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1530        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1531        HKLFT[:3] = np.where((HKLFT[2]>=0)&(HKLFT[1]<HKLFT[0]),np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1532        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([-1,1,-1])[:,nxs],HKLFT[:3]) 
1533        HKLFT[:3] = np.where((HKLFT[0]<0)&(HKLFT[2]<-HKLFT[0])&(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matd3q[nxs,:,:])).T,HKLFT[:3])
1534        HKLFT[:3] = np.where((HKLFT[0]<0)&(HKLFT[2]>=-HKLFT[0])&(HKLFT[1]>HKLFT[2]),np.squeeze(np.inner(HKLF[:,:3],matdm3[nxs,:,:])).T,HKLFT[:3])
1535    #centrosymmetric
1536    elif Laue == 'm 3':
1537        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1538        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1539        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])           
1540        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1541        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1542    elif Laue == 'm 3 m':
1543        HKLFT[:3] = np.where(HKLFT[0]<0,HKLFT[:3]*np.array([-1,1,1])[:,nxs],HKLFT[:3])
1544        HKLFT[:3] = np.where(HKLFT[1]<0,HKLFT[:3]*np.array([1,-1,1])[:,nxs],HKLFT[:3])
1545        HKLFT[:3] = np.where(HKLFT[2]<0,HKLFT[:3]*np.array([1,1,-1])[:,nxs],HKLFT[:3])           
1546        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1547        HKLFT[:3] = np.where((HKLFT[2]>=0)&((HKLFT[0]>=HKLFT[2])|(HKLFT[1]>HKLFT[2])),np.squeeze(np.inner(HKLF[:,:3],matd3[nxs,:,:])).T,HKLFT[:3])
1548        HKLFT[:3] = np.where(HKLFT[0]>HKLFT[1],np.squeeze(np.inner(HKLF[:,:3],matdm[nxs,:,:])).T,HKLFT[:3])
1549    return HKLFT.T
1550       
1551
1552#Spherical harmonics routines
1553def OdfChk(SGLaue,L,M):
1554    'needs doc string'
1555    if not L%2 and abs(M) <= L:
1556        if SGLaue == '0':                      #cylindrical symmetry
1557            if M == 0: return True
1558        elif SGLaue == '-1':
1559            return True
1560        elif SGLaue == '2/m':
1561            if not abs(M)%2: return True
1562        elif SGLaue == 'mmm':
1563            if not abs(M)%2 and M >= 0: return True
1564        elif SGLaue == '4/m':
1565            if not abs(M)%4: return True
1566        elif SGLaue == '4/mmm':
1567            if not abs(M)%4 and M >= 0: return True
1568        elif SGLaue in ['3R','3']:
1569            if not abs(M)%3: return True
1570        elif SGLaue in ['3mR','3m1','31m']:
1571            if not abs(M)%3 and M >= 0: return True
1572        elif SGLaue == '6/m':
1573            if not abs(M)%6: return True
1574        elif SGLaue == '6/mmm':
1575            if not abs(M)%6 and M >= 0: return True
1576        elif SGLaue == 'm3':
1577            if M > 0:
1578                if L%12 == 2:
1579                    if M <= L//12: return True
1580                else:
1581                    if M <= L//12+1: return True
1582        elif SGLaue == 'm3m':
1583            if M > 0:
1584                if L%12 == 2:
1585                    if M <= L//12: return True
1586                else:
1587                    if M <= L//12+1: return True
1588    return False
1589       
1590def GenSHCoeff(SGLaue,SamSym,L,IfLMN=True):
1591    'needs doc string'
1592    coeffNames = []
1593    for iord in [2*i+2 for i in range(L//2)]:
1594        for m in [i-iord for i in range(2*iord+1)]:
1595            if OdfChk(SamSym,iord,m):
1596                for n in [i-iord for i in range(2*iord+1)]:
1597                    if OdfChk(SGLaue,iord,n):
1598                        if IfLMN:
1599                            coeffNames.append('C(%d,%d,%d)'%(iord,m,n))
1600                        else:
1601                            coeffNames.append('C(%d,%d)'%(iord,n))
1602    return coeffNames
1603   
1604def CrsAng(H,cell,SGData):
1605    'needs doc string'
1606    a,b,c,al,be,ga = cell
1607    SQ3 = 1.732050807569
1608    H1 = np.array([1,0,0])
1609    H2 = np.array([0,1,0])
1610    H3 = np.array([0,0,1])
1611    H4 = np.array([1,1,1])
1612    G,g = cell2Gmat(cell)
1613    Laue = SGData['SGLaue']
1614    Naxis = SGData['SGUniq']
1615    if len(H.shape) == 1:
1616        DH = np.inner(H,np.inner(G,H))
1617    else:
1618        DH = np.array([np.inner(h,np.inner(G,h)) for h in H])
1619    if Laue == '2/m':
1620        if Naxis == 'a':
1621            DR = np.inner(H1,np.inner(G,H1))
1622            DHR = np.inner(H,np.inner(G,H1))
1623        elif Naxis == 'b':
1624            DR = np.inner(H2,np.inner(G,H2))
1625            DHR = np.inner(H,np.inner(G,H2))
1626        else:
1627            DR = np.inner(H3,np.inner(G,H3))
1628            DHR = np.inner(H,np.inner(G,H3))
1629    elif Laue in ['R3','R3m']:
1630        DR = np.inner(H4,np.inner(G,H4))
1631        DHR = np.inner(H,np.inner(G,H4))
1632    else:
1633        DR = np.inner(H3,np.inner(G,H3))
1634        DHR = np.inner(H,np.inner(G,H3))
1635    DHR /= np.sqrt(DR*DH)
1636    phi = np.where(DHR <= 1.0,acosd(DHR),0.0)
1637    if Laue == '-1':
1638        BA = H.T[1]*a/(b-H.T[0]*cosd(ga))
1639        BB = H.T[0]*sind(ga)**2
1640    elif Laue == '2/m':
1641        if Naxis == 'a':
1642            BA = H.T[2]*b/(c-H.T[1]*cosd(al))
1643            BB = H.T[1]*sind(al)**2
1644        elif Naxis == 'b':
1645            BA = H.T[0]*c/(a-H.T[2]*cosd(be))
1646            BB = H.T[2]*sind(be)**2
1647        else:
1648            BA = H.T[1]*a/(b-H.T[0]*cosd(ga))
1649            BB = H.T[0]*sind(ga)**2
1650    elif Laue in ['mmm','4/m','4/mmm']:
1651        BA = H.T[1]*a
1652        BB = H.T[0]*b
1653    elif Laue in ['3R','3mR']:
1654        BA = H.T[0]+H.T[1]-2.0*H.T[2]
1655        BB = SQ3*(H.T[0]-H.T[1])
1656    elif Laue in ['m3','m3m']:
1657        BA = H.T[1]
1658        BB = H.T[0]
1659    else:
1660        BA = H.T[0]+2.0*H.T[1]
1661        BB = SQ3*H.T[0]
1662    beta = atan2d(BA,BB)
1663    return phi,beta
1664   
1665def SamAng(Tth,Gangls,Sangl,IFCoup):
1666    """Compute sample orientation angles vs laboratory coord. system
1667
1668    :param Tth:        Signed theta                                   
1669    :param Gangls:     Sample goniometer angles phi,chi,omega,azmuth 
1670    :param Sangl:      Sample angle zeros om-0, chi-0, phi-0         
1671    :param IFCoup:     True if omega & 2-theta coupled in CW scan
1672    :returns: 
1673        psi,gam:    Sample odf angles                             
1674        dPSdA,dGMdA:    Angle zero derivatives
1675    """                         
1676   
1677    if IFCoup:
1678        GSomeg = sind(Gangls[2]+Tth)
1679        GComeg = cosd(Gangls[2]+Tth)
1680    else:
1681        GSomeg = sind(Gangls[2])
1682        GComeg = cosd(Gangls[2])
1683    GSTth = sind(Tth)
1684    GCTth = cosd(Tth)     
1685    GSazm = sind(Gangls[3])
1686    GCazm = cosd(Gangls[3])
1687    GSchi = sind(Gangls[1])
1688    GCchi = cosd(Gangls[1])
1689    GSphi = sind(Gangls[0]+Sangl[2])
1690    GCphi = cosd(Gangls[0]+Sangl[2])
1691    SSomeg = sind(Sangl[0])
1692    SComeg = cosd(Sangl[0])
1693    SSchi = sind(Sangl[1])
1694    SCchi = cosd(Sangl[1])
1695    AT = -GSTth*GComeg+GCTth*GCazm*GSomeg
1696    BT = GSTth*GSomeg+GCTth*GCazm*GComeg
1697    CT = -GCTth*GSazm*GSchi
1698    DT = -GCTth*GSazm*GCchi
1699   
1700    BC1 = -AT*GSphi+(CT+BT*GCchi)*GCphi
1701    BC2 = DT-BT*GSchi
1702    BC3 = AT*GCphi+(CT+BT*GCchi)*GSphi
1703     
1704    BC = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg     
1705    psi = acosd(BC)
1706   
1707    BD = 1.0-BC**2
1708    C = np.where(BD>1.e-6,rpd/np.sqrt(BD),0.)
1709    dPSdA = [-C*(-BC1*SSomeg*SCchi-BC2*SSomeg*SSchi-BC3*SComeg),
1710        -C*(-BC1*SComeg*SSchi+BC2*SComeg*SCchi),
1711        -C*(-BC1*SSomeg-BC3*SComeg*SCchi)]
1712     
1713    BA = -BC1*SSchi+BC2*SCchi
1714    BB = BC1*SSomeg*SCchi+BC2*SSomeg*SSchi+BC3*SComeg
1715    gam = atan2d(BB,BA)
1716
1717    BD = (BA**2+BB**2)/rpd
1718
1719    dBAdO = 0
1720    dBAdC = -BC1*SCchi-BC2*SSchi
1721    dBAdF = BC3*SSchi
1722   
1723    dBBdO = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg
1724    dBBdC = -BC1*SSomeg*SSchi+BC2*SSomeg*SCchi
1725    dBBdF = BC1*SComeg-BC3*SSomeg*SCchi
1726   
1727    dGMdA = np.where(BD > 1.e-6,[(BA*dBBdO-BB*dBAdO)/BD,(BA*dBBdC-BB*dBAdC)/BD, \
1728        (BA*dBBdF-BB*dBAdF)/BD],[np.zeros_like(BD),np.zeros_like(BD),np.zeros_like(BD)])
1729       
1730    return psi,gam,dPSdA,dGMdA
1731
1732BOH = {
1733'L=2':[[],[],[]],
1734'L=4':[[0.30469720,0.36418281],[],[]],
1735'L=6':[[-0.14104740,0.52775103],[],[]],
1736'L=8':[[0.28646862,0.21545346,0.32826995],[],[]],
1737'L=10':[[-0.16413497,0.33078546,0.39371345],[],[]],
1738'L=12':[[0.26141975,0.27266871,0.03277460,0.32589402],
1739    [0.09298802,-0.23773812,0.49446631,0.0],[]],
1740'L=14':[[-0.17557309,0.25821932,0.27709173,0.33645360],[],[]],
1741'L=16':[[0.24370673,0.29873515,0.06447688,0.00377,0.32574495],
1742    [0.12039646,-0.25330128,0.23950998,0.40962508,0.0],[]],
1743'L=18':[[-0.16914245,0.17017340,0.34598142,0.07433932,0.32696037],
1744    [-0.06901768,0.16006562,-0.24743528,0.47110273,0.0],[]],
1745'L=20':[[0.23067026,0.31151832,0.09287682,0.01089683,0.00037564,0.32573563],
1746    [0.13615420,-0.25048007,0.12882081,0.28642879,0.34620433,0.0],[]],
1747'L=22':[[-0.16109560,0.10244188,0.36285175,0.13377513,0.01314399,0.32585583],
1748    [-0.09620055,0.20244115,-0.22389483,0.17928946,0.42017231,0.0],[]],
1749'L=24':[[0.22050742,0.31770654,0.11661736,0.02049853,0.00150861,0.00003426,0.32573505],
1750    [0.13651722,-0.21386648,0.00522051,0.33939435,0.10837396,0.32914497,0.0],
1751    [0.05378596,-0.11945819,0.16272298,-0.26449730,0.44923956,0.0,0.0]],
1752'L=26':[[-0.15435003,0.05261630,0.35524646,0.18578869,0.03259103,0.00186197,0.32574594],
1753    [-0.11306511,0.22072681,-0.18706142,0.05439948,0.28122966,0.35634355,0.0],[]],
1754'L=28':[[0.21225019,0.32031716,0.13604702,0.03132468,0.00362703,0.00018294,0.00000294,0.32573501],
1755    [0.13219496,-0.17206256,-0.08742608,0.32671661,0.17973107,0.02567515,0.32619598,0.0],
1756    [0.07989184,-0.16735346,0.18839770,-0.20705337,0.12926808,0.42715602,0.0,0.0]],
1757'L=30':[[-0.14878368,0.01524973,0.33628434,0.22632587,0.05790047,0.00609812,0.00022898,0.32573594],
1758    [-0.11721726,0.20915005,-0.11723436,-0.07815329,0.31318947,0.13655742,0.33241385,0.0],
1759    [-0.04297703,0.09317876,-0.11831248,0.17355132,-0.28164031,0.42719361,0.0,0.0]],
1760'L=32':[[0.20533892,0.32087437,0.15187897,0.04249238,0.00670516,0.00054977,0.00002018,0.00000024,0.32573501],
1761    [0.12775091,-0.13523423,-0.14935701,0.28227378,0.23670434,0.05661270,0.00469819,0.32578978,0.0],
1762    [0.09703829,-0.19373733,0.18610682,-0.14407046,0.00220535,0.26897090,0.36633402,0.0,0.0]],
1763'L=34':[[-0.14409234,-0.01343681,0.31248977,0.25557722,0.08571889,0.01351208,0.00095792,0.00002550,0.32573508],
1764    [-0.11527834,0.18472133,-0.04403280,-0.16908618,0.27227021,0.21086614,0.04041752,0.32688152,0.0],
1765    [-0.06773139,0.14120811,-0.15835721,0.18357456,-0.19364673,0.08377174,0.43116318,0.0,0.0]]
1766}
1767
1768Lnorm = lambda L: 4.*np.pi/(2.0*L+1.)
1769
1770def GetKcl(L,N,SGLaue,phi,beta):
1771    'needs doc string'
1772    import pytexture as ptx
1773    if SGLaue in ['m3','m3m']:
1774        if 'array' in str(type(phi)) and np.any(phi.shape):
1775            Kcl = np.zeros_like(phi)
1776        else:
1777            Kcl = 0.
1778        for j in range(0,L+1,4):
1779            im = j//4
1780            if 'array' in str(type(phi)) and np.any(phi.shape):
1781                pcrs = ptx.pyplmpsi(L,j,len(phi),phi)[0]
1782            else:
1783                pcrs = ptx.pyplmpsi(L,j,1,phi)[0]
1784            Kcl += BOH['L=%d'%(L)][N-1][im]*pcrs*cosd(j*beta)       
1785    else:
1786        if 'array' in str(type(phi)) and np.any(phi.shape):
1787            pcrs = ptx.pyplmpsi(L,N,len(phi),phi)[0]
1788        else:
1789            pcrs = ptx.pyplmpsi(L,N,1,phi)[0]
1790        pcrs *= RSQ2PI
1791        if N:
1792            pcrs *= SQ2
1793        if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1794            if SGLaue in ['3mR','3m1','31m']: 
1795                if N%6 == 3:
1796                    Kcl = pcrs*sind(N*beta)
1797                else:
1798                    Kcl = pcrs*cosd(N*beta)
1799            else:
1800                Kcl = pcrs*cosd(N*beta)
1801        else:
1802            Kcl = pcrs*(cosd(N*beta)+sind(N*beta))
1803    return Kcl
1804   
1805def GetKsl(L,M,SamSym,psi,gam):
1806    'needs doc string'
1807    import pytexture as ptx
1808    if 'array' in str(type(psi)) and np.any(psi.shape):
1809        psrs,dpdps = ptx.pyplmpsi(L,M,len(psi),psi)
1810    else:
1811        psrs,dpdps = ptx.pyplmpsi(L,M,1,psi)
1812    psrs *= RSQ2PI
1813    dpdps *= RSQ2PI
1814    if M:
1815        psrs *= SQ2
1816        dpdps *= SQ2
1817    if SamSym in ['mmm',]:
1818        dum = cosd(M*gam)
1819        Ksl = psrs*dum
1820        dKsdp = dpdps*dum
1821        dKsdg = -psrs*M*sind(M*gam)
1822    else:
1823        dum = cosd(M*gam)+sind(M*gam)
1824        Ksl = psrs*dum
1825        dKsdp = dpdps*dum
1826        dKsdg = psrs*M*(-sind(M*gam)+cosd(M*gam))
1827    return Ksl,dKsdp,dKsdg
1828   
1829def GetKclKsl(L,N,SGLaue,psi,phi,beta):
1830    """
1831    This is used for spherical harmonics description of preferred orientation;
1832        cylindrical symmetry only (M=0) and no sample angle derivatives returned
1833    """
1834    import pytexture as ptx
1835    Ksl,x = ptx.pyplmpsi(L,0,1,psi)
1836    Ksl *= RSQ2PI
1837    if SGLaue in ['m3','m3m']:
1838        Kcl = 0.0
1839        for j in range(0,L+1,4):
1840            im = j//4
1841            pcrs,dum = ptx.pyplmpsi(L,j,1,phi)
1842            Kcl += BOH['L=%d'%(L)][N-1][im]*pcrs*cosd(j*beta)       
1843    else:
1844        pcrs,dum = ptx.pyplmpsi(L,N,1,phi)
1845        pcrs *= RSQ2PI
1846        if N:
1847            pcrs *= SQ2
1848        if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1849            if SGLaue in ['3mR','3m1','31m']: 
1850                if N%6 == 3:
1851                    Kcl = pcrs*sind(N*beta)
1852                else:
1853                    Kcl = pcrs*cosd(N*beta)
1854            else:
1855                Kcl = pcrs*cosd(N*beta)
1856        else:
1857            Kcl = pcrs*(cosd(N*beta)+sind(N*beta))
1858    return Kcl*Ksl,Lnorm(L)
1859   
1860def Glnh(Start,SHCoef,psi,gam,SamSym):
1861    'needs doc string'
1862    import pytexture as ptx
1863
1864    if Start:
1865        ptx.pyqlmninit()
1866        Start = False
1867    Fln = np.zeros(len(SHCoef))
1868    for i,term in enumerate(SHCoef):
1869        l,m,n = eval(term.strip('C'))
1870        pcrs,dum = ptx.pyplmpsi(l,m,1,psi)
1871        pcrs *= RSQPI
1872        if m == 0:
1873            pcrs /= SQ2
1874        if SamSym in ['mmm',]:
1875            Ksl = pcrs*cosd(m*gam)
1876        else:
1877            Ksl = pcrs*(cosd(m*gam)+sind(m*gam))
1878        Fln[i] = SHCoef[term]*Ksl*Lnorm(l)
1879    ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
1880    return ODFln
1881
1882def Flnh(Start,SHCoef,phi,beta,SGData):
1883    'needs doc string'
1884    import pytexture as ptx
1885   
1886    if Start:
1887        ptx.pyqlmninit()
1888        Start = False
1889    Fln = np.zeros(len(SHCoef))
1890    for i,term in enumerate(SHCoef):
1891        l,m,n = eval(term.strip('C'))
1892        if SGData['SGLaue'] in ['m3','m3m']:
1893            Kcl = 0.0
1894            for j in range(0,l+1,4):
1895                im = j//4
1896                pcrs,dum = ptx.pyplmpsi(l,j,1,phi)
1897                Kcl += BOH['L='+str(l)][n-1][im]*pcrs*cosd(j*beta)       
1898        else:                #all but cubic
1899            pcrs,dum = ptx.pyplmpsi(l,n,1,phi)
1900            pcrs *= RSQPI
1901            if n == 0:
1902                pcrs /= SQ2
1903            if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1904               if SGData['SGLaue'] in ['3mR','3m1','31m']: 
1905                   if n%6 == 3:
1906                       Kcl = pcrs*sind(n*beta)
1907                   else:
1908                       Kcl = pcrs*cosd(n*beta)
1909               else:
1910                   Kcl = pcrs*cosd(n*beta)
1911            else:
1912                Kcl = pcrs*(cosd(n*beta)+sind(n*beta))
1913        Fln[i] = SHCoef[term]*Kcl*Lnorm(l)
1914    ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
1915    return ODFln
1916   
1917def polfcal(ODFln,SamSym,psi,gam):
1918    '''Perform a pole figure computation.
1919    Note that the the number of gam values must either be 1 or must
1920    match psi. Updated for numpy 1.8.0
1921    '''
1922    import pytexture as ptx
1923    PolVal = np.ones_like(psi)
1924    for term in ODFln:
1925        if abs(ODFln[term][1]) > 1.e-3:
1926            l,m,n = eval(term.strip('C'))
1927            psrs,dum = ptx.pyplmpsi(l,m,len(psi),psi)
1928            if SamSym in ['-1','2/m']:
1929                if m:
1930                    Ksl = RSQPI*psrs*(cosd(m*gam)+sind(m*gam))
1931                else:
1932                    Ksl = RSQPI*psrs/SQ2
1933            else:
1934                if m:
1935                    Ksl = RSQPI*psrs*cosd(m*gam)
1936                else:
1937                    Ksl = RSQPI*psrs/SQ2
1938            PolVal += ODFln[term][1]*Ksl
1939    return PolVal
1940   
1941def invpolfcal(ODFln,SGData,phi,beta):
1942    'needs doc string'
1943    import pytexture as ptx
1944   
1945    invPolVal = np.ones_like(beta)
1946    for term in ODFln:
1947        if abs(ODFln[term][1]) > 1.e-3:
1948            l,m,n = eval(term.strip('C'))
1949            if SGData['SGLaue'] in ['m3','m3m']:
1950                Kcl = 0.0
1951                for j in range(0,l+1,4):
1952                    im = j//4
1953                    pcrs,dum = ptx.pyplmpsi(l,j,len(beta),phi)
1954                    Kcl += BOH['L=%d'%(l)][n-1][im]*pcrs*cosd(j*beta)       
1955            else:                #all but cubic
1956                pcrs,dum = ptx.pyplmpsi(l,n,len(beta),phi)
1957                pcrs *= RSQPI
1958                if n == 0:
1959                    pcrs /= SQ2
1960                if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1961                   if SGData['SGLaue'] in ['3mR','3m1','31m']: 
1962                       if n%6 == 3:
1963                           Kcl = pcrs*sind(n*beta)
1964                       else:
1965                           Kcl = pcrs*cosd(n*beta)
1966                   else:
1967                       Kcl = pcrs*cosd(n*beta)
1968                else:
1969                    Kcl = pcrs*(cosd(n*beta)+sind(n*beta))
1970            invPolVal += ODFln[term][1]*Kcl
1971    return invPolVal
1972   
1973   
1974def textureIndex(SHCoef):
1975    'needs doc string'
1976    Tindx = 1.0
1977    for term in SHCoef:
1978        l = eval(term.strip('C'))[0]
1979        Tindx += SHCoef[term]**2/(2.0*l+1.)
1980    return Tindx
1981   
1982# self-test materials follow.
1983selftestlist = []
1984'''Defines a list of self-tests'''
1985selftestquiet = True
1986def _ReportTest():
1987    'Report name and doc string of current routine when ``selftestquiet`` is False'
1988    if not selftestquiet:
1989        import inspect
1990        caller = inspect.stack()[1][3]
1991        doc = eval(caller).__doc__
1992        if doc is not None:
1993            print('testing '+__file__+' with '+caller+' ('+doc+')')
1994        else:
1995            print('testing '+__file__()+" with "+caller)
1996NeedTestData = True
1997def TestData():
1998    array = np.array
1999    global NeedTestData
2000    NeedTestData = False
2001    global CellTestData
2002    # output from uctbx computed on platform darwin on 2010-05-28
2003    CellTestData = [
2004# cell, g, G, cell*, V, V*
2005  [(4, 4, 4, 90, 90, 90), 
2006   array([[  1.60000000e+01,   9.79717439e-16,   9.79717439e-16],
2007       [  9.79717439e-16,   1.60000000e+01,   9.79717439e-16],
2008       [  9.79717439e-16,   9.79717439e-16,   1.60000000e+01]]), array([[  6.25000000e-02,   3.82702125e-18,   3.82702125e-18],
2009       [  3.82702125e-18,   6.25000000e-02,   3.82702125e-18],
2010       [  3.82702125e-18,   3.82702125e-18,   6.25000000e-02]]), (0.25, 0.25, 0.25, 90.0, 90.0, 90.0), 64.0, 0.015625],
2011# cell, g, G, cell*, V, V*
2012  [(4.0999999999999996, 5.2000000000000002, 6.2999999999999998, 100, 80, 130), 
2013   array([[ 16.81      , -13.70423184,   4.48533243],
2014       [-13.70423184,  27.04      ,  -5.6887143 ],
2015       [  4.48533243,  -5.6887143 ,  39.69      ]]), array([[ 0.10206349,  0.05083339, -0.00424823],
2016       [ 0.05083339,  0.06344997,  0.00334956],
2017       [-0.00424823,  0.00334956,  0.02615544]]), (0.31947376387537696, 0.25189277536327803, 0.16172643497798223, 85.283666420376008, 94.716333579624006, 50.825714168082683), 100.98576357983838, 0.0099023858863968445],
2018# cell, g, G, cell*, V, V*
2019  [(3.5, 3.5, 6, 90, 90, 120), 
2020   array([[  1.22500000e+01,  -6.12500000e+00,   1.28587914e-15],
2021       [ -6.12500000e+00,   1.22500000e+01,   1.28587914e-15],
2022       [  1.28587914e-15,   1.28587914e-15,   3.60000000e+01]]), array([[  1.08843537e-01,   5.44217687e-02,   3.36690552e-18],
2023       [  5.44217687e-02,   1.08843537e-01,   3.36690552e-18],
2024       [  3.36690552e-18,   3.36690552e-18,   2.77777778e-02]]), (0.32991443953692895, 0.32991443953692895, 0.16666666666666669, 90.0, 90.0, 60.000000000000021), 63.652867178156257, 0.015710211406520427],
2025  ]
2026    global CoordTestData
2027    CoordTestData = [
2028# cell, ((frac, ortho),...)
2029  ((4,4,4,90,90,90,), [
2030 ((0.10000000000000001, 0.0, 0.0),(0.40000000000000002, 0.0, 0.0)),
2031 ((0.0, 0.10000000000000001, 0.0),(2.4492935982947065e-17, 0.40000000000000002, 0.0)),
2032 ((0.0, 0.0, 0.10000000000000001),(2.4492935982947065e-17, -2.4492935982947065e-17, 0.40000000000000002)),
2033 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.40000000000000013, 0.79999999999999993, 1.2)),
2034 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.80000000000000016, 1.2, 0.40000000000000002)),
2035 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(1.2, 0.80000000000000004, 0.40000000000000002)),
2036 ((0.5, 0.5, 0.5),(2.0, 1.9999999999999998, 2.0)),
2037]),
2038# cell, ((frac, ortho),...)
2039  ((4.1,5.2,6.3,100,80,130,), [
2040 ((0.10000000000000001, 0.0, 0.0),(0.40999999999999998, 0.0, 0.0)),
2041 ((0.0, 0.10000000000000001, 0.0),(-0.33424955703700043, 0.39834311042186865, 0.0)),
2042 ((0.0, 0.0, 0.10000000000000001),(0.10939835193016617, -0.051013289294572106, 0.6183281045774256)),
2043 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.069695941716497567, 0.64364635296002093, 1.8549843137322766)),
2044 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(-0.073350319180835066, 1.1440160419710339, 0.6183281045774256)),
2045 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.67089923785616512, 0.74567293154916525, 0.6183281045774256)),
2046 ((0.5, 0.5, 0.5),(0.92574397446582857, 1.7366491056364828, 3.0916405228871278)),
2047]),
2048# cell, ((frac, ortho),...)
2049  ((3.5,3.5,6,90,90,120,), [
2050 ((0.10000000000000001, 0.0, 0.0),(0.35000000000000003, 0.0, 0.0)),
2051 ((0.0, 0.10000000000000001, 0.0),(-0.17499999999999993, 0.3031088913245536, 0.0)),
2052 ((0.0, 0.0, 0.10000000000000001),(3.6739403974420595e-17, -3.6739403974420595e-17, 0.60000000000000009)),
2053 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(2.7675166561703527e-16, 0.60621778264910708, 1.7999999999999998)),
2054 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.17500000000000041, 0.90932667397366063, 0.60000000000000009)),
2055 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.70000000000000018, 0.6062177826491072, 0.60000000000000009)),
2056 ((0.5, 0.5, 0.5),(0.87500000000000067, 1.5155444566227676, 3.0)),
2057]),
2058]
2059    global LaueTestData             #generated by GSAS
2060    LaueTestData = {
2061    'R 3 m':[(4.,4.,6.,90.,90.,120.),((1,0,1,6),(1,0,-2,6),(0,0,3,2),(1,1,0,6),(2,0,-1,6),(2,0,2,6),
2062        (1,1,3,12),(1,0,4,6),(2,1,1,12),(2,1,-2,12),(3,0,0,6),(1,0,-5,6),(2,0,-4,6),(3,0,-3,6),(3,0,3,6),
2063        (0,0,6,2),(2,2,0,6),(2,1,4,12),(2,0,5,6),(3,1,-1,12),(3,1,2,12),(1,1,6,12),(2,2,3,12),(2,1,-5,12))],
2064    'R 3':[(4.,4.,6.,90.,90.,120.),((1,0,1,6),(1,0,-2,6),(0,0,3,2),(1,1,0,6),(2,0,-1,6),(2,0,2,6),(1,1,3,6),
2065        (1,1,-3,6),(1,0,4,6),(3,-1,1,6),(2,1,1,6),(3,-1,-2,6),(2,1,-2,6),(3,0,0,6),(1,0,-5,6),(2,0,-4,6),
2066        (2,2,0,6),(3,0,3,6),(3,0,-3,6),(0,0,6,2),(3,-1,4,6),(2,0,5,6),(2,1,4,6),(4,-1,-1,6),(3,1,-1,6),
2067        (3,1,2,6),(4,-1,2,6),(2,2,-3,6),(1,1,-6,6),(1,1,6,6),(2,2,3,6),(2,1,-5,6),(3,-1,-5,6))],
2068    'P 3':[(4.,4.,6.,90.,90.,120.),((0,0,1,2),(1,0,0,6),(1,0,1,6),(0,0,2,2),(1,0,-1,6),(1,0,2,6),(1,0,-2,6),
2069        (1,1,0,6),(0,0,3,2),(1,1,1,6),(1,1,-1,6),(1,0,3,6),(1,0,-3,6),(2,0,0,6),(2,0,-1,6),(1,1,-2,6),
2070        (1,1,2,6),(2,0,1,6),(2,0,-2,6),(2,0,2,6),(0,0,4,2),(1,1,-3,6),(1,1,3,6),(1,0,-4,6),(1,0,4,6),
2071        (2,0,-3,6),(2,1,0,6),(2,0,3,6),(3,-1,0,6),(2,1,1,6),(3,-1,-1,6),(2,1,-1,6),(3,-1,1,6),(1,1,4,6),
2072        (3,-1,2,6),(3,-1,-2,6),(1,1,-4,6),(0,0,5,2),(2,1,2,6),(2,1,-2,6),(3,0,0,6),(3,0,1,6),(2,0,4,6),
2073        (2,0,-4,6),(3,0,-1,6),(1,0,-5,6),(1,0,5,6),(3,-1,-3,6),(2,1,-3,6),(2,1,3,6),(3,-1,3,6),(3,0,-2,6),
2074        (3,0,2,6),(1,1,5,6),(1,1,-5,6),(2,2,0,6),(3,0,3,6),(3,0,-3,6),(0,0,6,2),(2,0,-5,6),(2,1,-4,6),
2075        (2,2,-1,6),(3,-1,-4,6),(2,2,1,6),(3,-1,4,6),(2,1,4,6),(2,0,5,6),(1,0,-6,6),(1,0,6,6),(4,-1,0,6),
2076        (3,1,0,6),(3,1,-1,6),(3,1,1,6),(4,-1,-1,6),(2,2,2,6),(4,-1,1,6),(2,2,-2,6),(3,1,2,6),(3,1,-2,6),
2077        (3,0,4,6),(3,0,-4,6),(4,-1,-2,6),(4,-1,2,6),(2,2,-3,6),(1,1,6,6),(1,1,-6,6),(2,2,3,6),(3,-1,5,6),
2078        (2,1,5,6),(2,1,-5,6),(3,-1,-5,6))],
2079    'P 3 m 1':[(4.,4.,6.,90.,90.,120.),((0,0,1,2),(1,0,0,6),(1,0,-1,6),(1,0,1,6),(0,0,2,2),(1,0,-2,6),
2080        (1,0,2,6),(1,1,0,6),(0,0,3,2),(1,1,1,12),(1,0,-3,6),(1,0,3,6),(2,0,0,6),(1,1,2,12),(2,0,1,6),
2081        (2,0,-1,6),(0,0,4,2),(2,0,-2,6),(2,0,2,6),(1,1,3,12),(1,0,-4,6),(1,0,4,6),(2,0,3,6),(2,1,0,12),
2082        (2,0,-3,6),(2,1,1,12),(2,1,-1,12),(1,1,4,12),(2,1,2,12),(0,0,5,2),(2,1,-2,12),(3,0,0,6),(1,0,-5,6),
2083        (3,0,1,6),(3,0,-1,6),(1,0,5,6),(2,0,4,6),(2,0,-4,6),(2,1,3,12),(2,1,-3,12),(3,0,-2,6),(3,0,2,6),
2084        (1,1,5,12),(3,0,-3,6),(0,0,6,2),(2,2,0,6),(3,0,3,6),(2,1,4,12),(2,2,1,12),(2,0,5,6),(2,1,-4,12),
2085        (2,0,-5,6),(1,0,-6,6),(1,0,6,6),(3,1,0,12),(3,1,-1,12),(3,1,1,12),(2,2,2,12),(3,1,2,12),
2086        (3,0,4,6),(3,1,-2,12),(3,0,-4,6),(1,1,6,12),(2,2,3,12))],
2087    'P 3 1 m':[(4.,4.,6.,90.,90.,120.),((0,0,1,2),(1,0,0,6),(0,0,2,2),(1,0,1,12),(1,0,2,12),(1,1,0,6),
2088        (0,0,3,2),(1,1,-1,6),(1,1,1,6),(1,0,3,12),(2,0,0,6),(2,0,1,12),(1,1,2,6),(1,1,-2,6),(2,0,2,12),
2089        (0,0,4,2),(1,1,-3,6),(1,1,3,6),(1,0,4,12),(2,1,0,12),(2,0,3,12),(2,1,1,12),(2,1,-1,12),(1,1,-4,6),
2090        (1,1,4,6),(0,0,5,2),(2,1,-2,12),(2,1,2,12),(3,0,0,6),(1,0,5,12),(2,0,4,12),(3,0,1,12),(2,1,-3,12),
2091        (2,1,3,12),(3,0,2,12),(1,1,5,6),(1,1,-5,6),(3,0,3,12),(0,0,6,2),(2,2,0,6),(2,1,-4,12),(2,0,5,12),
2092        (2,2,-1,6),(2,2,1,6),(2,1,4,12),(3,1,0,12),(1,0,6,12),(2,2,2,6),(3,1,-1,12),(2,2,-2,6),(3,1,1,12),
2093        (3,1,-2,12),(3,0,4,12),(3,1,2,12),(1,1,-6,6),(2,2,3,6),(2,2,-3,6),(1,1,6,6))],
2094    }
2095   
2096    global FLnhTestData
2097    FLnhTestData = [{
2098    'C(4,0,0)': (0.965, 0.42760447),
2099    'C(2,0,0)': (1.0122, -0.80233610),
2100    'C(2,0,2)': (0.0061, 8.37491546E-03),
2101    'C(6,0,4)': (-0.0898, 4.37985696E-02),
2102    'C(6,0,6)': (-0.1369, -9.04081762E-02),
2103    'C(6,0,0)': (0.5935, -0.18234928),
2104    'C(4,0,4)': (0.1872, 0.16358127),
2105    'C(6,0,2)': (0.6193, 0.27573633),
2106    'C(4,0,2)': (-0.1897, 0.12530720)},[1,0,0]]
2107def test0():
2108    if NeedTestData: TestData()
2109    msg = 'test cell2Gmat, fillgmat, Gmat2cell'
2110    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2111        G, g = cell2Gmat(cell)
2112        assert np.allclose(G,tG),msg
2113        assert np.allclose(g,tg),msg
2114        tcell = Gmat2cell(g)
2115        assert np.allclose(cell,tcell),msg
2116        tcell = Gmat2cell(G)
2117        assert np.allclose(tcell,trcell),msg
2118if __name__ == '__main__': selftestlist.append(test0)
2119
2120def test1():
2121    'test cell2A and A2Gmat'
2122    _ReportTest()
2123    if NeedTestData: TestData()
2124    msg = 'test cell2A and A2Gmat'
2125    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2126        G, g = A2Gmat(cell2A(cell))
2127        assert np.allclose(G,tG),msg
2128        assert np.allclose(g,tg),msg
2129if __name__ == '__main__': selftestlist.append(test1)
2130
2131def test2():
2132    'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
2133    _ReportTest()
2134    if NeedTestData: TestData()
2135    msg = 'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
2136    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2137        G, g = cell2Gmat(cell)
2138        tcell = A2cell(Gmat2A(G))
2139        assert np.allclose(cell,tcell),msg
2140if __name__ == '__main__': selftestlist.append(test2)
2141
2142def test3():
2143    'test invcell2Gmat'
2144    _ReportTest()
2145    if NeedTestData: TestData()
2146    msg = 'test invcell2Gmat'
2147    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2148        G, g = invcell2Gmat(trcell)
2149        assert np.allclose(G,tG),msg
2150        assert np.allclose(g,tg),msg
2151if __name__ == '__main__': selftestlist.append(test3)
2152
2153def test4():
2154    'test calc_rVsq, calc_rV, calc_V'
2155    _ReportTest()
2156    if NeedTestData: TestData()
2157    msg = 'test calc_rVsq, calc_rV, calc_V'
2158    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2159        assert np.allclose(calc_rV(cell2A(cell)),trV), msg
2160        assert np.allclose(calc_V(cell2A(cell)),tV), msg
2161if __name__ == '__main__': selftestlist.append(test4)
2162
2163def test5():
2164    'test A2invcell'
2165    _ReportTest()
2166    if NeedTestData: TestData()
2167    msg = 'test A2invcell'
2168    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
2169        rcell = A2invcell(cell2A(cell))
2170        assert np.allclose(rcell,trcell),msg
2171if __name__ == '__main__': selftestlist.append(test5)
2172
2173def test6():
2174    'test cell2AB'
2175    _ReportTest()
2176    if NeedTestData: TestData()
2177    msg = 'test cell2AB'
2178    for (cell,coordlist) in CoordTestData:
2179        A,B = cell2AB(cell)
2180        for (frac,ortho) in coordlist:
2181            to = np.inner(A,frac)
2182            tf = np.inner(B,to)
2183            assert np.allclose(ortho,to), msg
2184            assert np.allclose(frac,tf), msg
2185            to = np.sum(A*frac,axis=1)
2186            tf = np.sum(B*to,axis=1)
2187            assert np.allclose(ortho,to), msg
2188            assert np.allclose(frac,tf), msg
2189if __name__ == '__main__': selftestlist.append(test6)
2190
2191def test7():
2192    'test GetBraviasNum(...) and GenHBravais(...)'
2193    _ReportTest()
2194    import os.path
2195    import sys
2196    import GSASIIspc as spc
2197    testdir = os.path.join(os.path.split(os.path.abspath( __file__ ))[0],'testinp')
2198    if os.path.exists(testdir):
2199        if testdir not in sys.path: sys.path.insert(0,testdir)
2200    import sgtbxlattinp
2201    derror = 1e-4
2202    def indexmatch(hklin, hkllist, system):
2203        for hklref in hkllist:
2204            hklref = list(hklref)
2205            # these permutations are far from complete, but are sufficient to
2206            # allow the test to complete
2207            if system == 'cubic':
2208                permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
2209            elif system == 'monoclinic':
2210                permlist = [(1,2,3),(-1,2,-3)]
2211            else:
2212                permlist = [(1,2,3)]
2213
2214            for perm in permlist:
2215                hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
2216                if hkl == hklref: return True
2217                if [-i for i in hkl] == hklref: return True
2218        else:
2219            return False
2220
2221    for key in sgtbxlattinp.sgtbx7:
2222        spdict = spc.SpcGroup(key)
2223        cell = sgtbxlattinp.sgtbx7[key][0]
2224        system = spdict[1]['SGSys']
2225        center = spdict[1]['SGLatt']
2226
2227        bravcode = GetBraviasNum(center, system)
2228
2229        g2list = GenHBravais(sgtbxlattinp.dmin, bravcode, cell2A(cell))
2230
2231        assert len(sgtbxlattinp.sgtbx7[key][1]) == len(g2list), 'Reflection lists differ for %s' % key
2232        for h,k,l,d,num in g2list:
2233            for hkllist,dref in sgtbxlattinp.sgtbx7[key][1]: 
2234                if abs(d-dref) < derror:
2235                    if indexmatch((h,k,l,), hkllist, system):
2236                        break
2237            else:
2238                assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
2239if __name__ == '__main__': selftestlist.append(test7)
2240
2241def test8():
2242    'test GenHLaue'
2243    _ReportTest()
2244    import GSASIIspc as spc
2245    import sgtbxlattinp
2246    derror = 1e-4
2247    dmin = sgtbxlattinp.dmin
2248
2249    def indexmatch(hklin, hklref, system, axis):
2250        # these permutations are far from complete, but are sufficient to
2251        # allow the test to complete
2252        if system == 'cubic':
2253            permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
2254        elif system == 'monoclinic' and axis=='b':
2255            permlist = [(1,2,3),(-1,2,-3)]
2256        elif system == 'monoclinic' and axis=='a':
2257            permlist = [(1,2,3),(1,-2,-3)]
2258        elif system == 'monoclinic' and axis=='c':
2259            permlist = [(1,2,3),(-1,-2,3)]
2260        elif system == 'trigonal':
2261            permlist = [(1,2,3),(2,1,3),(-1,-2,3),(-2,-1,3)]
2262        elif system == 'rhombohedral':
2263            permlist = [(1,2,3),(2,3,1),(3,1,2)]
2264        else:
2265            permlist = [(1,2,3)]
2266
2267        hklref = list(hklref)
2268        for perm in permlist:
2269            hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
2270            if hkl == hklref: return True
2271            if [-i for i in hkl] == hklref: return True
2272        return False
2273
2274    for key in sgtbxlattinp.sgtbx8:
2275        spdict = spc.SpcGroup(key)[1]
2276        cell = sgtbxlattinp.sgtbx8[key][0]
2277        Axis = spdict['SGUniq']
2278        system = spdict['SGSys']
2279
2280        g2list = GenHLaue(dmin,spdict,cell2A(cell))
2281        #if len(g2list) != len(sgtbxlattinp.sgtbx8[key][1]):
2282        #    print 'failed',key,':' ,len(g2list),'vs',len(sgtbxlattinp.sgtbx8[key][1])
2283        #    print 'GSAS-II:'
2284        #    for h,k,l,d in g2list: print '  ',(h,k,l),d
2285        #    print 'SGTBX:'
2286        #    for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print '  ',hkllist,dref
2287        assert len(g2list) == len(sgtbxlattinp.sgtbx8[key][1]), (
2288            'Reflection lists differ for %s' % key
2289            )
2290        #match = True
2291        for h,k,l,d in g2list:
2292            for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: 
2293                if abs(d-dref) < derror:
2294                    if indexmatch((h,k,l,), hkllist, system, Axis): break
2295            else:
2296                assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
2297                #match = False
2298        #if not match:
2299            #for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print '  ',hkllist,dref
2300            #print center, Laue, Axis, system
2301if __name__ == '__main__': selftestlist.append(test8)
2302           
2303def test9():
2304    'test GenHLaue'
2305    _ReportTest()
2306    import GSASIIspc as G2spc
2307    if NeedTestData: TestData()
2308    for spc in LaueTestData:
2309        data = LaueTestData[spc]
2310        cell = data[0]
2311        hklm = np.array(data[1])
2312        H = hklm[-1][:3]
2313        hklO = hklm.T[:3].T
2314        A = cell2A(cell)
2315        dmin = 1./np.sqrt(calc_rDsq(H,A))
2316        SGData = G2spc.SpcGroup(spc)[1]
2317        hkls = np.array(GenHLaue(dmin,SGData,A))
2318        hklN = hkls.T[:3].T
2319        #print spc,hklO.shape,hklN.shape
2320        err = True
2321        for H in hklO:
2322            if H not in hklN:
2323                print ('%d %s'%(H,' missing from hkl from GSASII'))
2324                err = False
2325        assert(err)
2326if __name__ == '__main__': selftestlist.append(test9)
2327       
2328       
2329   
2330
2331if __name__ == '__main__':
2332    # run self-tests
2333    selftestquiet = False
2334    for test in selftestlist:
2335        test()
2336    print ("OK")
Note: See TracBrowser for help on using the repository browser.