source: trunk/GSASIIlattice.py @ 2154

Last change on this file since 2154 was 2154, checked in by vondreele, 7 years ago

remove range restriction on LGmix - inhibited doing refinements
do wx.CallAfter? for size/strain plots in G2ddataGUI
Add TransformDialog?, prodMGMT, TransformCell?, TransformU6, TransformXYZ to G2lattice
A start on TransformAtoms? in G2math
fill in OnTransform?

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