source: trunk/GSASIIlattice.py @ 989

Last change on this file since 989 was 989, checked in by toby, 10 years ago

add constraint derivs for single xtal; allow cellFill to work without esds; update docs; fix controls copy bug
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 56.7 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: 2013-07-14 03:21:12 +0000 (Sun, 14 Jul 2013) $
27# $Author: toby $
28# $Revision: 989 $
29# $URL: trunk/GSASIIlattice.py $
30# $Id: GSASIIlattice.py 989 2013-07-14 03:21:12Z toby $
31########### SVN repository information ###################
32import math
33import numpy as np
34import numpy.linalg as nl
35
36# trig functions in degrees
37sind = lambda x: np.sin(x*np.pi/180.)
38asind = lambda x: 180.*np.arcsin(x)/np.pi
39tand = lambda x: np.tan(x*np.pi/180.)
40atand = lambda x: 180.*np.arctan(x)/np.pi
41atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
42cosd = lambda x: np.cos(x*np.pi/180.)
43acosd = lambda x: 180.*np.arccos(x)/np.pi
44rdsq2d = lambda x,p: round(1.0/np.sqrt(x),p)
45
46def sec2HMS(sec):
47    """Convert time in sec to H:M:S string
48   
49    :param sec: time in seconds
50    :return: H:M:S string (to nearest 100th second)
51   
52    """
53    H = int(sec/3600)
54    M = int(sec/60-H*60)
55    S = sec-3600*H-60*M
56    return '%d:%2d:%.2f'%(H,M,S)
57   
58def rotdMat(angle,axis=0):
59    """Prepare rotation matrix for angle in degrees about axis(=0,1,2)
60
61    :param angle: angle in degrees
62    :param axis:  axis (0,1,2 = x,y,z) about which for the rotation
63    :return: rotation matrix - 3x3 numpy array
64
65    """
66    if axis == 2:
67        return np.array([[cosd(angle),-sind(angle),0],[sind(angle),cosd(angle),0],[0,0,1]])
68    elif axis == 1:
69        return np.array([[cosd(angle),0,-sind(angle)],[0,1,0],[sind(angle),0,cosd(angle)]])
70    else:
71        return np.array([[1,0,0],[0,cosd(angle),-sind(angle)],[0,sind(angle),cosd(angle)]])
72       
73def rotdMat4(angle,axis=0):
74    """Prepare rotation matrix for angle in degrees about axis(=0,1,2) with scaling for OpenGL
75
76    :param angle: angle in degrees
77    :param axis:  axis (0,1,2 = x,y,z) about which for the rotation
78    :return: rotation matrix - 4x4 numpy array (last row/column for openGL scaling)
79
80    """
81    Mat = rotdMat(angle,axis)
82    return np.concatenate((np.concatenate((Mat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
83   
84def fillgmat(cell):
85    """Compute lattice metric tensor from unit cell constants
86
87    :param cell: tuple with a,b,c,alpha, beta, gamma (degrees)
88    :return: 3x3 numpy array
89
90    """
91    a,b,c,alp,bet,gam = cell
92    g = np.array([
93        [a*a,  a*b*cosd(gam),  a*c*cosd(bet)],
94        [a*b*cosd(gam),  b*b,  b*c*cosd(alp)],
95        [a*c*cosd(bet) ,b*c*cosd(alp),   c*c]])
96    return g
97           
98def cell2Gmat(cell):
99    """Compute real and reciprocal lattice metric tensor from unit cell constants
100
101    :param cell: tuple with a,b,c,alpha, beta, gamma (degrees)
102    :return: reciprocal (G) & real (g) metric tensors (list of two numpy 3x3 arrays)
103
104    """
105    g = fillgmat(cell)
106    G = nl.inv(g)       
107    return G,g
108
109def A2Gmat(A,inverse=True):
110    """Fill real & reciprocal metric tensor (G) from A.
111
112    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
113    :param bool inverse: if True return both G and g; else just G
114    :return: reciprocal (G) & real (g) metric tensors (list of two numpy 3x3 arrays)
115
116    """
117    G = np.zeros(shape=(3,3))
118    G = [
119        [A[0],  A[3]/2.,  A[4]/2.], 
120        [A[3]/2.,A[1],    A[5]/2.], 
121        [A[4]/2.,A[5]/2.,    A[2]]]
122    if inverse:
123        g = nl.inv(G)
124        return G,g
125    else:
126        return G
127
128def Gmat2A(G):
129    """Extract A from reciprocal metric tensor (G)
130
131    :param G: reciprocal maetric tensor (3x3 numpy array
132    :return: A = [G11,G22,G33,2*G12,2*G13,2*G23]
133
134    """
135    return [G[0][0],G[1][1],G[2][2],2.*G[0][1],2.*G[0][2],2.*G[1][2]]
136   
137def cell2A(cell):
138    """Obtain A = [G11,G22,G33,2*G12,2*G13,2*G23] from lattice parameters
139
140    :param cell: [a,b,c,alpha,beta,gamma] (degrees)
141    :return: G reciprocal metric tensor as 3x3 numpy array
142
143    """
144    G,g = cell2Gmat(cell)
145    return Gmat2A(G)
146
147def A2cell(A):
148    """Compute unit cell constants from A
149
150    :param A: [G11,G22,G33,2*G12,2*G13,2*G23] G - reciprocal metric tensor
151    :return: a,b,c,alpha, beta, gamma (degrees) - lattice parameters
152
153    """
154    G,g = A2Gmat(A)
155    return Gmat2cell(g)
156
157def Gmat2cell(g):
158    """Compute real/reciprocal lattice parameters from real/reciprocal metric tensor (g/G)
159    The math works the same either way.
160
161    :param g (or G): real (or reciprocal) metric tensor 3x3 array
162    :return: a,b,c,alpha, beta, gamma (degrees) (or a*,b*,c*,alpha*,beta*,gamma* degrees)
163
164    """
165    oldset = np.seterr('raise')
166    a = np.sqrt(max(0,g[0][0]))
167    b = np.sqrt(max(0,g[1][1]))
168    c = np.sqrt(max(0,g[2][2]))
169    alp = acosd(g[2][1]/(b*c))
170    bet = acosd(g[2][0]/(a*c))
171    gam = acosd(g[0][1]/(a*b))
172    np.seterr(**oldset)
173    return a,b,c,alp,bet,gam
174
175def invcell2Gmat(invcell):
176    """Compute real and reciprocal lattice metric tensor from reciprocal
177       unit cell constants
178       
179    :param invcell: [a*,b*,c*,alpha*, beta*, gamma*] (degrees)
180    :return: reciprocal (G) & real (g) metric tensors (list of two 3x3 arrays)
181
182    """
183    G = fillgmat(invcell)
184    g = nl.inv(G)
185    return G,g
186       
187def calc_rVsq(A):
188    """Compute the square of the reciprocal lattice volume (1/V**2) from A'
189
190    """
191    G,g = A2Gmat(A)
192    rVsq = nl.det(G)
193    if rVsq < 0:
194        return 1
195    return rVsq
196   
197def calc_rV(A):
198    """Compute the reciprocal lattice volume (V*) from A
199    """
200    return np.sqrt(calc_rVsq(A))
201   
202def calc_V(A):
203    """Compute the real lattice volume (V) from A
204    """
205    return 1./calc_rV(A)
206
207def A2invcell(A):
208    """Compute reciprocal unit cell constants from A
209    returns tuple with a*,b*,c*,alpha*, beta*, gamma* (degrees)
210    """
211    G,g = A2Gmat(A)
212    return Gmat2cell(G)
213   
214def Gmat2AB(G):
215    """Computes orthogonalization matrix from reciprocal metric tensor G
216
217    :returns: tuple of two 3x3 numpy arrays (A,B)
218
219       * A for crystal to Cartesian transformations A*x = np.inner(A,x) = X
220       * B (= inverse of A) for Cartesian to crystal transformation B*X = np.inner(B,X) = x
221
222    """
223    cellstar = Gmat2cell(G)
224    g = nl.inv(G)
225    cell = Gmat2cell(g)
226    A = np.zeros(shape=(3,3))
227    # from Giacovazzo (Fundamentals 2nd Ed.) p.75
228    A[0][0] = cell[0]                # a
229    A[0][1] = cell[1]*cosd(cell[5])  # b cos(gamma)
230    A[0][2] = cell[2]*cosd(cell[4])  # c cos(beta)
231    A[1][1] = cell[1]*sind(cell[5])  # b sin(gamma)
232    A[1][2] = -cell[2]*cosd(cellstar[3])*sind(cell[4]) # - c cos(alpha*) sin(beta)
233    A[2][2] = 1/cellstar[2]         # 1/c*
234    B = nl.inv(A)
235    return A,B
236   
237
238def cell2AB(cell):
239    """Computes orthogonalization matrix from unit cell constants
240
241    :param tuple cell: a,b,c, alpha, beta, gamma (degrees)
242    :returns: tuple of two 3x3 numpy arrays (A,B)
243       A for crystal to Cartesian transformations A*x = np.inner(A,x) = X
244       B (= inverse of A) for Cartesian to crystal transformation B*X = np.inner(B,X) = x
245    """
246    G,g = cell2Gmat(cell) 
247    cellstar = Gmat2cell(G)
248    A = np.zeros(shape=(3,3))
249    # from Giacovazzo (Fundamentals 2nd Ed.) p.75
250    A[0][0] = cell[0]                # a
251    A[0][1] = cell[1]*cosd(cell[5])  # b cos(gamma)
252    A[0][2] = cell[2]*cosd(cell[4])  # c cos(beta)
253    A[1][1] = cell[1]*sind(cell[5])  # b sin(gamma)
254    A[1][2] = -cell[2]*cosd(cellstar[3])*sind(cell[4]) # - c cos(alpha*) sin(beta)
255    A[2][2] = 1/cellstar[2]         # 1/c*
256    B = nl.inv(A)
257    return A,B
258   
259def U6toUij(U6):
260    """Fill matrix (Uij) from U6 = [U11,U22,U33,U12,U13,U23]
261    NB: there is a non numpy version in GSASIIspc: U2Uij
262
263    :param list U6: 6 terms of u11,u22,...
264    :returns:
265        Uij - numpy [3][3] array of uij
266    """
267    U = np.array([
268        [U6[0],  U6[3],  U6[4]], 
269        [U6[3],  U6[1],  U6[5]], 
270        [U6[4],  U6[5],  U6[2]]])
271    return U
272
273def UijtoU6(U):
274    """Fill vector [U11,U22,U33,U12,U13,U23] from Uij
275    NB: there is a non numpy version in GSASIIspc: Uij2U
276    """
277    U6 = np.array([U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]])
278    return U6
279
280def Uij2betaij(Uij,G):
281    """
282    Convert Uij to beta-ij tensors -- stub for eventual completion
283   
284    :param Uij: numpy array [Uij]
285    :param G: reciprocal metric tensor
286    :returns: beta-ij - numpy array [beta-ij]
287    """
288    pass
289   
290def cell2GS(cell):
291    ''' returns Uij to betaij conversion matrix'''
292    G,g = cell2Gmat(cell)
293    GS = G
294    GS[0][1] = GS[1][0] = math.sqrt(GS[0][0]*GS[1][1])
295    GS[0][2] = GS[2][0] = math.sqrt(GS[0][0]*GS[2][2])
296    GS[1][2] = GS[2][1] = math.sqrt(GS[1][1]*GS[2][2])
297    return GS   
298   
299def Uij2Ueqv(Uij,GS,Amat):
300    ''' returns 1/3 trace of diagonalized U matrix'''
301    U = np.multiply(U6toUij(Uij),GS)
302    U = np.inner(Amat,np.inner(U,Amat).T)
303    E,R = nl.eigh(U)
304    return np.sum(E)/3.
305       
306def CosSinAngle(U,V,G):
307    """ calculate sin & cos of angle between U & V in generalized coordinates
308    defined by metric tensor G
309
310    :param U: 3-vectors assume numpy arrays
311    :param V: 3-vectors assume numpy arrays
312    :param G: metric tensor for U & V defined space assume numpy array
313    :returns:
314        cos(phi) & sin(phi)
315    """
316    u = U/np.sqrt(np.inner(U,np.inner(G,U)))
317    v = V/np.sqrt(np.inner(V,np.inner(G,V)))
318    cosP = np.inner(u,np.inner(G,v))
319    sinP = np.sqrt(max(0.0,1.0-cosP**2))
320    return cosP,sinP
321   
322def criticalEllipse(prob):
323    """
324    Calculate critical values for probability ellipsoids from probability
325    """
326    if not ( 0.01 <= prob < 1.0):
327        return 1.54 
328    coeff = np.array([6.44988E-09,4.16479E-07,1.11172E-05,1.58767E-04,0.00130554,
329        0.00604091,0.0114921,-0.040301,-0.6337203,1.311582])
330    llpr = math.log(-math.log(prob))
331    return np.polyval(coeff,llpr)
332   
333def CellBlock(nCells):
334    """
335    Generate block of unit cells n*n*n on a side; [0,0,0] centered, n = 2*nCells+1
336    currently only works for nCells = 0 or 1 (not >1)
337    """
338    if nCells:
339        N = 2*nCells+1
340        N2 = N*N
341        N3 = N*N*N
342        cellArray = []
343        A = np.array(range(N3))
344        cellGen = np.array([A/N2-1,A/N%N-1,A%N-1]).T
345        for cell in cellGen:
346            cellArray.append(cell)
347        return cellArray
348    else:
349        return [0,0,0]
350       
351def CellAbsorption(ElList,Volume):
352    '''Compute unit cell absorption
353
354    :param dict ElList: dictionary of element contents including mu and
355      number of atoms be cell
356    :param float Volume: unit cell volume
357    :returns: mu-total/Volume
358    '''
359    muT = 0
360    for El in ElList:
361        muT += ElList[El]['mu']*ElList[El]['FormulaNo']
362    return muT/Volume
363   
364#Permutations and Combinations
365# Four routines: combinations,uniqueCombinations, selections & permutations
366#These taken from Python Cookbook, 2nd Edition. 19.15 p724-726
367#   
368def _combinators(_handle, items, n):
369    """ factored-out common structure of all following combinators """
370    if n==0:
371        yield [ ]
372        return
373    for i, item in enumerate(items):
374        this_one = [ item ]
375        for cc in _combinators(_handle, _handle(items, i), n-1):
376            yield this_one + cc
377def combinations(items, n):
378    """ take n distinct items, order matters """
379    def skipIthItem(items, i):
380        return items[:i] + items[i+1:]
381    return _combinators(skipIthItem, items, n)
382def uniqueCombinations(items, n):
383    """ take n distinct items, order is irrelevant """
384    def afterIthItem(items, i):
385        return items[i+1:]
386    return _combinators(afterIthItem, items, n)
387def selections(items, n):
388    """ take n (not necessarily distinct) items, order matters """
389    def keepAllItems(items, i):
390        return items
391    return _combinators(keepAllItems, items, n)
392def permutations(items):
393    """ take all items, order matters """
394    return combinations(items, len(items))
395
396#reflection generation routines
397#for these: H = [h,k,l]; A is as used in calc_rDsq; G - inv metric tensor, g - metric tensor;
398#           cell - a,b,c,alp,bet,gam in A & deg
399   
400def calc_rDsq(H,A):
401    'needs doc string'
402    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]
403    return rdsq
404   
405def calc_rDsq2(H,G):
406    'needs doc string'
407    return np.inner(H,np.inner(G,H))
408   
409def calc_rDsqZ(H,A,Z,tth,lam):
410    'needs doc string'
411    rpd = np.pi/180.
412    rdsq = calc_rDsq(H,A)+Z*sind(tth)*2.0*rpd/lam**2
413    return rdsq
414       
415def MaxIndex(dmin,A):
416    'needs doc string'
417    Hmax = [0,0,0]
418    try:
419        cell = A2cell(A)
420    except:
421        cell = [1,1,1,90,90,90]
422    for i in range(3):
423        Hmax[i] = int(round(cell[i]/dmin))
424    return Hmax
425   
426def sortHKLd(HKLd,ifreverse,ifdup):
427    '''needs doc string
428
429    :param HKLd: a list of [h,k,l,d,...];
430    :param ifreverse: True for largest d first
431    :param ifdup: True if duplicate d-spacings allowed
432    '''
433    T = []
434    for i,H in enumerate(HKLd):
435        if ifdup:
436            T.append((H[3],i))
437        else:
438            T.append(H[3])           
439    D = dict(zip(T,HKLd))
440    T.sort()
441    if ifreverse:
442        T.reverse()
443    X = []
444    okey = ''
445    for key in T: 
446        if key != okey: X.append(D[key])    #remove duplicate d-spacings
447        okey = key
448    return X
449   
450def SwapIndx(Axis,H):
451    'needs doc string'
452    if Axis in [1,-1]:
453        return H
454    elif Axis in [2,-3]:
455        return [H[1],H[2],H[0]]
456    else:
457        return [H[2],H[0],H[1]]
458       
459def Rh2Hx(Rh):
460    'needs doc string'
461    Hx = [0,0,0]
462    Hx[0] = Rh[0]-Rh[1]
463    Hx[1] = Rh[1]-Rh[2]
464    Hx[2] = np.sum(Rh)
465    return Hx
466   
467def Hx2Rh(Hx):
468    'needs doc string'
469    Rh = [0,0,0]
470    itk = -Hx[0]+Hx[1]+Hx[2]
471    if itk%3 != 0:
472        return 0        #error - not rhombohedral reflection
473    else:
474        Rh[1] = itk/3
475        Rh[0] = Rh[1]+Hx[0]
476        Rh[2] = Rh[1]-Hx[1]
477        if Rh[0] < 0:
478            for i in range(3):
479                Rh[i] = -Rh[i]
480        return Rh
481       
482def CentCheck(Cent,H):
483    'needs doc string'
484    h,k,l = H
485    if Cent == 'A' and (k+l)%2:
486        return False
487    elif Cent == 'B' and (h+l)%2:
488        return False
489    elif Cent == 'C' and (h+k)%2:
490        return False
491    elif Cent == 'I' and (h+k+l)%2:
492        return False
493    elif Cent == 'F' and ((h+k)%2 or (h+l)%2 or (k+l)%2):
494        return False
495    elif Cent == 'R' and (-h+k+l)%3:
496        return False
497    else:
498        return True
499                                   
500def GetBraviasNum(center,system):
501    """Determine the Bravais lattice number, as used in GenHBravais
502   
503    :param center: one of: 'P', 'C', 'I', 'F', 'R' (see SGLatt from GSASIIspc.SpcGroup)
504    :param system: one of 'cubic', 'hexagonal', 'tetragonal', 'orthorhombic', 'trigonal' (for R)
505      'monoclinic', 'triclinic' (see SGSys from GSASIIspc.SpcGroup)
506    :return: a number between 0 and 13
507      or throws a ValueError exception if the combination of center, system is not found (i.e. non-standard)
508
509    """
510    if center.upper() == 'F' and system.lower() == 'cubic':
511        return 0
512    elif center.upper() == 'I' and system.lower() == 'cubic':
513        return 1
514    elif center.upper() == 'P' and system.lower() == 'cubic':
515        return 2
516    elif center.upper() == 'R' and system.lower() == 'trigonal':
517        return 3
518    elif center.upper() == 'P' and system.lower() == 'hexagonal':
519        return 4
520    elif center.upper() == 'I' and system.lower() == 'tetragonal':
521        return 5
522    elif center.upper() == 'P' and system.lower() == 'tetragonal':
523        return 6
524    elif center.upper() == 'F' and system.lower() == 'orthorhombic':
525        return 7
526    elif center.upper() == 'I' and system.lower() == 'orthorhombic':
527        return 8
528    elif center.upper() == 'C' and system.lower() == 'orthorhombic':
529        return 9
530    elif center.upper() == 'P' and system.lower() == 'orthorhombic':
531        return 10
532    elif center.upper() == 'C' and system.lower() == 'monoclinic':
533        return 11
534    elif center.upper() == 'P' and system.lower() == 'monoclinic':
535        return 12
536    elif center.upper() == 'P' and system.lower() == 'triclinic':
537        return 13
538    raise ValueError,'non-standard Bravais lattice center=%s, cell=%s' % (center,system)
539
540def GenHBravais(dmin,Bravais,A):
541    """Generate the positionally unique powder diffraction reflections
542     
543    :param dmin: minimum d-spacing in A
544    :param Bravais: lattice type (see GetBraviasNum). Bravais is one of::
545             0 F cubic
546             1 I cubic
547             2 P cubic
548             3 R hexagonal (trigonal not rhombohedral)
549             4 P hexagonal
550             5 I tetragonal
551             6 P tetragonal
552             7 F orthorhombic
553             8 I orthorhombic
554             9 C orthorhombic
555             10 P orthorhombic
556             11 C monoclinic
557             12 P monoclinic
558             13 P triclinic
559           
560    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
561    :return: HKL unique d list of [h,k,l,d,-1] sorted with largest d first
562           
563    """
564    import math
565    if Bravais in [9,11]:
566        Cent = 'C'
567    elif Bravais in [1,5,8]:
568        Cent = 'I'
569    elif Bravais in [0,7]:
570        Cent = 'F'
571    elif Bravais in [3]:
572        Cent = 'R'
573    else:
574        Cent = 'P'
575    Hmax = MaxIndex(dmin,A)
576    dminsq = 1./(dmin**2)
577    HKL = []
578    if Bravais == 13:                       #triclinic
579        for l in range(-Hmax[2],Hmax[2]+1):
580            for k in range(-Hmax[1],Hmax[1]+1):
581                hmin = 0
582                if (k < 0): hmin = 1
583                if (k ==0 and l < 0): hmin = 1
584                for h in range(hmin,Hmax[0]+1):
585                    H=[h,k,l]
586                    rdsq = calc_rDsq(H,A)
587                    if 0 < rdsq <= dminsq:
588                        HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
589    elif Bravais in [11,12]:                #monoclinic - b unique
590        Hmax = SwapIndx(2,Hmax)
591        for h in range(Hmax[0]+1):
592            for k in range(-Hmax[1],Hmax[1]+1):
593                lmin = 0
594                if k < 0:lmin = 1
595                for l in range(lmin,Hmax[2]+1):
596                    [h,k,l] = SwapIndx(-2,[h,k,l])
597                    H = []
598                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
599                    if H:
600                        rdsq = calc_rDsq(H,A)
601                        if 0 < rdsq <= dminsq:
602                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
603                    [h,k,l] = SwapIndx(2,[h,k,l])
604    elif Bravais in [7,8,9,10]:            #orthorhombic
605        for h in range(Hmax[0]+1):
606            for k in range(Hmax[1]+1):
607                for l in range(Hmax[2]+1):
608                    H = []
609                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
610                    if H:
611                        rdsq = calc_rDsq(H,A)
612                        if 0 < rdsq <= dminsq:
613                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
614    elif Bravais in [5,6]:                  #tetragonal
615        for l in range(Hmax[2]+1):
616            for k in range(Hmax[1]+1):
617                for h in range(k,Hmax[0]+1):
618                    H = []
619                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
620                    if H:
621                        rdsq = calc_rDsq(H,A)
622                        if 0 < rdsq <= dminsq:
623                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
624    elif Bravais in [3,4]:
625        lmin = 0
626        if Bravais == 3: lmin = -Hmax[2]                  #hexagonal/trigonal
627        for l in range(lmin,Hmax[2]+1):
628            for k in range(Hmax[1]+1):
629                hmin = k
630                if l < 0: hmin += 1
631                for h in range(hmin,Hmax[0]+1):
632                    H = []
633                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
634                    if H:
635                        rdsq = calc_rDsq(H,A)
636                        if 0 < rdsq <= dminsq:
637                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
638
639    else:                                   #cubic
640        for l in range(Hmax[2]+1):
641            for k in range(l,Hmax[1]+1):
642                for h in range(k,Hmax[0]+1):
643                    H = []
644                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
645                    if H:
646                        rdsq = calc_rDsq(H,A)
647                        if 0 < rdsq <= dminsq:
648                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
649    return sortHKLd(HKL,True,False)
650   
651def getHKLmax(dmin,SGData,A):
652    'finds maximum allowed hkl for given A within dmin'
653    SGLaue = SGData['SGLaue']
654    if SGLaue in ['3R','3mR']:        #Rhombohedral axes
655        Hmax = [0,0,0]
656        cell = A2cell(A)
657        aHx = cell[0]*math.sqrt(2.0*(1.0-cosd(cell[3])))
658        cHx = cell[0]*math.sqrt(3.0*(1.0+2.0*cosd(cell[3])))
659        Hmax[0] = Hmax[1] = int(round(aHx/dmin))
660        Hmax[2] = int(round(cHx/dmin))
661        #print Hmax,aHx,cHx
662    else:                           # all others
663        Hmax = MaxIndex(dmin,A)
664    return Hmax
665   
666def GenHLaue(dmin,SGData,A):
667    """Generate the crystallographically unique powder diffraction reflections
668    for a lattice and Bravais type
669   
670    :param dmin: minimum d-spacing
671    :param SGData: space group dictionary with at least
672   
673        * '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'
674        * 'SGLatt': lattice centering: one of 'P','A','B','C','I','F'
675        * 'SGUniq': code for unique monoclinic axis one of 'a','b','c' (only if 'SGLaue' is '2/m') otherwise an empty string
676       
677    :param A: reciprocal metric tensor elements as [G11,G22,G33,2*G12,2*G13,2*G23]
678    :return: HKL = list of [h,k,l,d] sorted with largest d first and is unique
679            part of reciprocal space ignoring anomalous dispersion
680           
681    """
682    import math
683    SGLaue = SGData['SGLaue']
684    SGLatt = SGData['SGLatt']
685    SGUniq = SGData['SGUniq']
686    #finds maximum allowed hkl for given A within dmin
687    Hmax = getHKLmax(dmin,SGData,A)
688       
689    dminsq = 1./(dmin**2)
690    HKL = []
691    if SGLaue == '-1':                       #triclinic
692        for l in range(-Hmax[2],Hmax[2]+1):
693            for k in range(-Hmax[1],Hmax[1]+1):
694                hmin = 0
695                if (k < 0) or (k ==0 and l < 0): hmin = 1
696                for h in range(hmin,Hmax[0]+1):
697                    H = []
698                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
699                    if H:
700                        rdsq = calc_rDsq(H,A)
701                        if 0 < rdsq <= dminsq:
702                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
703    elif SGLaue == '2/m':                #monoclinic
704        axisnum = 1 + ['a','b','c'].index(SGUniq)
705        Hmax = SwapIndx(axisnum,Hmax)
706        for h in range(Hmax[0]+1):
707            for k in range(-Hmax[1],Hmax[1]+1):
708                lmin = 0
709                if k < 0:lmin = 1
710                for l in range(lmin,Hmax[2]+1):
711                    [h,k,l] = SwapIndx(-axisnum,[h,k,l])
712                    H = []
713                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
714                    if H:
715                        rdsq = calc_rDsq(H,A)
716                        if 0 < rdsq <= dminsq:
717                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
718                    [h,k,l] = SwapIndx(axisnum,[h,k,l])
719    elif SGLaue in ['mmm','4/m','6/m']:            #orthorhombic
720        for l in range(Hmax[2]+1):
721            for h in range(Hmax[0]+1):
722                kmin = 1
723                if SGLaue == 'mmm' or h ==0: kmin = 0
724                for k in range(kmin,Hmax[1]+1):
725                    H = []
726                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
727                    if H:
728                        rdsq = calc_rDsq(H,A)
729                        if 0 < rdsq <= dminsq:
730                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
731    elif SGLaue in ['4/mmm','6/mmm']:                  #tetragonal & hexagonal
732        for l in range(Hmax[2]+1):
733            for h in range(Hmax[0]+1):
734                for k in range(h+1):
735                    H = []
736                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
737                    if H:
738                        rdsq = calc_rDsq(H,A)
739                        if 0 < rdsq <= dminsq:
740                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
741    elif SGLaue in ['3m1','31m','3','3R','3mR']:                  #trigonals
742        for l in range(-Hmax[2],Hmax[2]+1):
743            hmin = 0
744            if l < 0: hmin = 1
745            for h in range(hmin,Hmax[0]+1):
746                if SGLaue in ['3R','3']:
747                    kmax = h
748                    kmin = -int((h-1.)/2.)
749                else:
750                    kmin = 0
751                    kmax = h
752                    if SGLaue in ['3m1','3mR'] and l < 0: kmax = h-1
753                    if SGLaue == '31m' and l < 0: kmin = 1
754                for k in range(kmin,kmax+1):
755                    H = []
756                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
757                    if SGLaue in ['3R','3mR']:
758                        H = Hx2Rh(H)
759                    if H:
760                        rdsq = calc_rDsq(H,A)
761                        if 0 < rdsq <= dminsq:
762                            HKL.append([H[0],H[1],H[2],1/math.sqrt(rdsq)])
763    else:                                   #cubic
764        for h in range(Hmax[0]+1):
765            for k in range(h+1):
766                lmin = 0
767                lmax = k
768                if SGLaue =='m3':
769                    lmax = h-1
770                    if h == k: lmax += 1
771                for l in range(lmin,lmax+1):
772                    H = []
773                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
774                    if H:
775                        rdsq = calc_rDsq(H,A)
776                        if 0 < rdsq <= dminsq:
777                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
778    return sortHKLd(HKL,True,True)
779
780#Spherical harmonics routines
781def OdfChk(SGLaue,L,M):
782    'needs doc string'
783    if not L%2 and abs(M) <= L:
784        if SGLaue == '0':                      #cylindrical symmetry
785            if M == 0: return True
786        elif SGLaue == '-1':
787            return True
788        elif SGLaue == '2/m':
789            if not abs(M)%2: return True
790        elif SGLaue == 'mmm':
791            if not abs(M)%2 and M >= 0: return True
792        elif SGLaue == '4/m':
793            if not abs(M)%4: return True
794        elif SGLaue == '4/mmm':
795            if not abs(M)%4 and M >= 0: return True
796        elif SGLaue in ['3R','3']:
797            if not abs(M)%3: return True
798        elif SGLaue in ['3mR','3m1','31m']:
799            if not abs(M)%3 and M >= 0: return True
800        elif SGLaue == '6/m':
801            if not abs(M)%6: return True
802        elif SGLaue == '6/mmm':
803            if not abs(M)%6 and M >= 0: return True
804        elif SGLaue == 'm3':
805            if M > 0:
806                if L%12 == 2:
807                    if M <= L/12: return True
808                else:
809                    if M <= L/12+1: return True
810        elif SGLaue == 'm3m':
811            if M > 0:
812                if L%12 == 2:
813                    if M <= L/12: return True
814                else:
815                    if M <= L/12+1: return True
816    return False
817       
818def GenSHCoeff(SGLaue,SamSym,L,IfLMN=True):
819    'needs doc string'
820    coeffNames = []
821    for iord in [2*i+2 for i in range(L/2)]:
822        for m in [i-iord for i in range(2*iord+1)]:
823            if OdfChk(SamSym,iord,m):
824                for n in [i-iord for i in range(2*iord+1)]:
825                    if OdfChk(SGLaue,iord,n):
826                        if IfLMN:
827                            coeffNames.append('C(%d,%d,%d)'%(iord,m,n))
828                        else:
829                            coeffNames.append('C(%d,%d)'%(iord,n))
830    return coeffNames
831   
832def CrsAng(H,cell,SGData):
833    'needs doc string'
834    a,b,c,al,be,ga = cell
835    SQ3 = 1.732050807569
836    H1 = np.array([1,0,0])
837    H2 = np.array([0,1,0])
838    H3 = np.array([0,0,1])
839    H4 = np.array([1,1,1])
840    G,g = cell2Gmat(cell)
841    Laue = SGData['SGLaue']
842    Naxis = SGData['SGUniq']
843    DH = np.inner(H,np.inner(G,H))
844    if Laue == '2/m':
845        if Naxis == 'a':
846            DR = np.inner(H1,np.inner(G,H1))
847            DHR = np.inner(H,np.inner(G,H1))
848        elif Naxis == 'b':
849            DR = np.inner(H2,np.inner(G,H2))
850            DHR = np.inner(H,np.inner(G,H2))
851        else:
852            DR = np.inner(H3,np.inner(G,H3))
853            DHR = np.inner(H,np.inner(G,H3))
854    elif Laue in ['R3','R3m']:
855        DR = np.inner(H4,np.inner(G,H4))
856        DHR = np.inner(H,np.inner(G,H4))
857       
858    else:
859        DR = np.inner(H3,np.inner(G,H3))
860        DHR = np.inner(H,np.inner(G,H3))
861    DHR /= np.sqrt(DR*DH)
862    phi = np.where(DHR <= 1.0,acosd(DHR),0.0)
863    if Laue == '-1':
864        BA = H[1]*a/(b-H[0]*cosd(ga))
865        BB = H[0]*sind(ga)**2
866    elif Laue == '2/m':
867        if Naxis == 'a':
868            BA = H[2]*b/(c-H[1]*cosd(al))
869            BB = H[1]*sind(al)**2
870        elif Naxis == 'b':
871            BA = H[0]*c/(a-H[2]*cosd(be))
872            BB = H[2]*sind(be)**2
873        else:
874            BA = H[1]*a/(b-H[0]*cosd(ga))
875            BB = H[0]*sind(ga)**2
876    elif Laue in ['mmm','4/m','4/mmm']:
877        BA = H[1]*a
878        BB = H[0]*b
879   
880    elif Laue in ['3R','3mR']:
881        BA = H[0]+H[1]-2.0*H[2]
882        BB = SQ3*(H[0]-H[1])
883    elif Laue in ['m3','m3m']:
884        BA = H[1]
885        BB = H[0]
886    else:
887        BA = H[0]+2.0*H[1]
888        BB = SQ3*H[0]
889    beta = atan2d(BA,BB)
890    return phi,beta
891   
892def SamAng(Tth,Gangls,Sangl,IFCoup):
893    """Compute sample orientation angles vs laboratory coord. system
894
895    :param Tth:        Signed theta                                   
896    :param Gangls:     Sample goniometer angles phi,chi,omega,azmuth 
897    :param Sangl:      Sample angle zeros om-0, chi-0, phi-0         
898    :param IFCoup:     True if omega & 2-theta coupled in CW scan
899    :returns: 
900        psi,gam:    Sample odf angles                             
901        dPSdA,dGMdA:    Angle zero derivatives
902    """                         
903   
904    rpd = math.pi/180.
905    if IFCoup:
906        GSomeg = sind(Gangls[2]+Tth)
907        GComeg = cosd(Gangls[2]+Tth)
908    else:
909        GSomeg = sind(Gangls[2])
910        GComeg = cosd(Gangls[2])
911    GSTth = sind(Tth)
912    GCTth = cosd(Tth)     
913    GSazm = sind(Gangls[3])
914    GCazm = cosd(Gangls[3])
915    GSchi = sind(Gangls[1])
916    GCchi = cosd(Gangls[1])
917    GSphi = sind(Gangls[0]+Sangl[2])
918    GCphi = cosd(Gangls[0]+Sangl[2])
919    SSomeg = sind(Sangl[0])
920    SComeg = cosd(Sangl[0])
921    SSchi = sind(Sangl[1])
922    SCchi = cosd(Sangl[1])
923    AT = -GSTth*GComeg+GCTth*GCazm*GSomeg
924    BT = GSTth*GSomeg+GCTth*GCazm*GComeg
925    CT = -GCTth*GSazm*GSchi
926    DT = -GCTth*GSazm*GCchi
927   
928    BC1 = -AT*GSphi+(CT+BT*GCchi)*GCphi
929    BC2 = DT-BT*GSchi
930    BC3 = AT*GCphi+(CT+BT*GCchi)*GSphi
931     
932    BC = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg     
933    psi = acosd(BC)
934   
935    BD = 1.0-BC**2
936    if BD > 0.:
937        C = rpd/math.sqrt(BD)
938    else:
939        C = 0.
940    dPSdA = [-C*(-BC1*SSomeg*SCchi-BC2*SSomeg*SSchi-BC3*SComeg),
941        -C*(-BC1*SComeg*SSchi+BC2*SComeg*SCchi),
942        -C*(-BC1*SSomeg-BC3*SComeg*SCchi)]
943     
944    BA = -BC1*SSchi+BC2*SCchi
945    BB = BC1*SSomeg*SCchi+BC2*SSomeg*SSchi+BC3*SComeg
946    gam = atan2d(BB,BA)
947
948    BD = (BA**2+BB**2)/rpd
949
950    dBAdO = 0
951    dBAdC = -BC1*SCchi-BC2*SSchi
952    dBAdF = BC3*SSchi
953   
954    dBBdO = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg
955    dBBdC = -BC1*SSomeg*SSchi+BC2*SSomeg*SCchi
956    dBBdF = BC1*SComeg-BC3*SSomeg*SCchi
957   
958    if BD > 0.:
959        dGMdA = [(BA*dBBdO-BB*dBAdO)/BD,(BA*dBBdC-BB*dBAdC)/BD,(BA*dBBdF-BB*dBAdF)/BD]
960    else:
961        dGMdA = [0.0,0.0,0.0]
962
963       
964    return psi,gam,dPSdA,dGMdA
965
966BOH = {
967'L=2':[[],[],[]],
968'L=4':[[0.30469720,0.36418281],[],[]],
969'L=6':[[-0.14104740,0.52775103],[],[]],
970'L=8':[[0.28646862,0.21545346,0.32826995],[],[]],
971'L=10':[[-0.16413497,0.33078546,0.39371345],[],[]],
972'L=12':[[0.26141975,0.27266871,0.03277460,0.32589402],
973    [0.09298802,-0.23773812,0.49446631,0.0],[]],
974'L=14':[[-0.17557309,0.25821932,0.27709173,0.33645360],[],[]],
975'L=16':[[0.24370673,0.29873515,0.06447688,0.00377,0.32574495],
976    [0.12039646,-0.25330128,0.23950998,0.40962508,0.0],[]],
977'L=18':[[-0.16914245,0.17017340,0.34598142,0.07433932,0.32696037],
978    [-0.06901768,0.16006562,-0.24743528,0.47110273,0.0],[]],
979'L=20':[[0.23067026,0.31151832,0.09287682,0.01089683,0.00037564,0.32573563],
980    [0.13615420,-0.25048007,0.12882081,0.28642879,0.34620433,0.0],[]],
981'L=22':[[-0.16109560,0.10244188,0.36285175,0.13377513,0.01314399,0.32585583],
982    [-0.09620055,0.20244115,-0.22389483,0.17928946,0.42017231,0.0],[]],
983'L=24':[[0.22050742,0.31770654,0.11661736,0.02049853,0.00150861,0.00003426,0.32573505],
984    [0.13651722,-0.21386648,0.00522051,0.33939435,0.10837396,0.32914497,0.0],
985    [0.05378596,-0.11945819,0.16272298,-0.26449730,0.44923956,0.0,0.0]],
986'L=26':[[-0.15435003,0.05261630,0.35524646,0.18578869,0.03259103,0.00186197,0.32574594],
987    [-0.11306511,0.22072681,-0.18706142,0.05439948,0.28122966,0.35634355,0.0],[]],
988'L=28':[[0.21225019,0.32031716,0.13604702,0.03132468,0.00362703,0.00018294,0.00000294,0.32573501],
989    [0.13219496,-0.17206256,-0.08742608,0.32671661,0.17973107,0.02567515,0.32619598,0.0],
990    [0.07989184,-0.16735346,0.18839770,-0.20705337,0.12926808,0.42715602,0.0,0.0]],
991'L=30':[[-0.14878368,0.01524973,0.33628434,0.22632587,0.05790047,0.00609812,0.00022898,0.32573594],
992    [-0.11721726,0.20915005,-0.11723436,-0.07815329,0.31318947,0.13655742,0.33241385,0.0],
993    [-0.04297703,0.09317876,-0.11831248,0.17355132,-0.28164031,0.42719361,0.0,0.0]],
994'L=32':[[0.20533892,0.32087437,0.15187897,0.04249238,0.00670516,0.00054977,0.00002018,0.00000024,0.32573501],
995    [0.12775091,-0.13523423,-0.14935701,0.28227378,0.23670434,0.05661270,0.00469819,0.32578978,0.0],
996    [0.09703829,-0.19373733,0.18610682,-0.14407046,0.00220535,0.26897090,0.36633402,0.0,0.0]],
997'L=34':[[-0.14409234,-0.01343681,0.31248977,0.25557722,0.08571889,0.01351208,0.00095792,0.00002550,0.32573508],
998    [-0.11527834,0.18472133,-0.04403280,-0.16908618,0.27227021,0.21086614,0.04041752,0.32688152,0.0],
999    [-0.06773139,0.14120811,-0.15835721,0.18357456,-0.19364673,0.08377174,0.43116318,0.0,0.0]]
1000}
1001
1002Lnorm = lambda L: 4.*np.pi/(2.0*L+1.)
1003
1004def GetKcl(L,N,SGLaue,phi,beta):
1005    'needs doc string'
1006    import pytexture as ptx
1007    RSQ2PI = 0.3989422804014
1008    SQ2 = 1.414213562373
1009    if SGLaue in ['m3','m3m']:
1010        Kcl = 0.0
1011        for j in range(0,L+1,4):
1012            im = j/4+1
1013            pcrs,dum = ptx.pyplmpsi(L,j,1,phi)
1014            Kcl += BOH['L='+str(L)][N-1][im-1]*pcrs*cosd(j*beta)       
1015    else:
1016        pcrs,dum = ptx.pyplmpsi(L,N,1,phi)
1017        pcrs *= RSQ2PI
1018        if N:
1019            pcrs *= SQ2
1020        if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1021            if SGLaue in ['3mR','3m1','31m']: 
1022                if N%6 == 3:
1023                    Kcl = pcrs*sind(N*beta)
1024                else:
1025                    Kcl = pcrs*cosd(N*beta)
1026            else:
1027                Kcl = pcrs*cosd(N*beta)
1028        else:
1029            Kcl = pcrs*(cosd(N*beta)+sind(N*beta))
1030    return Kcl
1031   
1032def GetKsl(L,M,SamSym,psi,gam):
1033    'needs doc string'
1034    import pytexture as ptx
1035    RSQPI = 0.5641895835478
1036    SQ2 = 1.414213562373
1037    psrs,dpdps = ptx.pyplmpsi(L,M,1,psi)
1038    psrs *= RSQPI
1039    dpdps *= RSQPI
1040    if M == 0:
1041        psrs /= SQ2
1042        dpdps /= SQ2
1043    if SamSym in ['mmm',]:
1044        dum = cosd(M*gam)
1045        Ksl = psrs*dum
1046        dKsdp = dpdps*dum
1047        dKsdg = -psrs*M*sind(M*gam)
1048    else:
1049        dum = cosd(M*gam)+sind(M*gam)
1050        Ksl = psrs*dum
1051        dKsdp = dpdps*dum
1052        dKsdg = psrs*M*(-sind(M*gam)+cosd(M*gam))
1053    return Ksl,dKsdp,dKsdg
1054   
1055def GetKclKsl(L,N,SGLaue,psi,phi,beta):
1056    """
1057    This is used for spherical harmonics description of preferred orientation;
1058        cylindrical symmetry only (M=0) and no sample angle derivatives returned
1059    """
1060    import pytexture as ptx
1061    RSQ2PI = 0.3989422804014
1062    SQ2 = 1.414213562373
1063    Ksl,x = ptx.pyplmpsi(L,0,1,psi)
1064    Ksl *= RSQ2PI
1065    if SGLaue in ['m3','m3m']:
1066        Kcl = 0.0
1067        for j in range(0,L+1,4):
1068            im = j/4+1
1069            pcrs,dum = ptx.pyplmpsi(L,j,1,phi)
1070            Kcl += BOH['L='+str(L)][N-1][im-1]*pcrs*cosd(j*beta)       
1071    else:
1072        pcrs,dum = ptx.pyplmpsi(L,N,1,phi)
1073        pcrs *= RSQ2PI
1074        if N:
1075            pcrs *= SQ2
1076        if SGLaue in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1077            if SGLaue in ['3mR','3m1','31m']: 
1078                if N%6 == 3:
1079                    Kcl = pcrs*sind(N*beta)
1080                else:
1081                    Kcl = pcrs*cosd(N*beta)
1082            else:
1083                Kcl = pcrs*cosd(N*beta)
1084        else:
1085            Kcl = pcrs*(cosd(N*beta)+sind(N*beta))
1086    return Kcl*Ksl,Lnorm(L)
1087   
1088def Glnh(Start,SHCoef,psi,gam,SamSym):
1089    'needs doc string'
1090    import pytexture as ptx
1091    RSQPI = 0.5641895835478
1092    SQ2 = 1.414213562373
1093
1094    if Start:
1095        ptx.pyqlmninit()
1096        Start = False
1097    Fln = np.zeros(len(SHCoef))
1098    for i,term in enumerate(SHCoef):
1099        l,m,n = eval(term.strip('C'))
1100        pcrs,dum = ptx.pyplmpsi(l,m,1,psi)
1101        pcrs *= RSQPI
1102        if m == 0:
1103            pcrs /= SQ2
1104        if SamSym in ['mmm',]:
1105            Ksl = pcrs*cosd(m*gam)
1106        else:
1107            Ksl = pcrs*(cosd(m*gam)+sind(m*gam))
1108        Fln[i] = SHCoef[term]*Ksl*Lnorm(l)
1109    ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
1110    return ODFln
1111
1112def Flnh(Start,SHCoef,phi,beta,SGData):
1113    'needs doc string'
1114    import pytexture as ptx
1115   
1116    FORPI = 12.5663706143592
1117    RSQPI = 0.5641895835478
1118    SQ2 = 1.414213562373
1119
1120    if Start:
1121        ptx.pyqlmninit()
1122        Start = False
1123    Fln = np.zeros(len(SHCoef))
1124    for i,term in enumerate(SHCoef):
1125        l,m,n = eval(term.strip('C'))
1126        if SGData['SGLaue'] in ['m3','m3m']:
1127            Kcl = 0.0
1128            for j in range(0,l+1,4):
1129                im = j/4+1
1130                pcrs,dum = ptx.pyplmpsi(l,j,1,phi)
1131                Kcl += BOH['L='+str(l)][n-1][im-1]*pcrs*cosd(j*beta)       
1132        else:                #all but cubic
1133            pcrs,dum = ptx.pyplmpsi(l,n,1,phi)
1134            pcrs *= RSQPI
1135            if n == 0:
1136                pcrs /= SQ2
1137            if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1138               if SGData['SGLaue'] in ['3mR','3m1','31m']: 
1139                   if n%6 == 3:
1140                       Kcl = pcrs*sind(n*beta)
1141                   else:
1142                       Kcl = pcrs*cosd(n*beta)
1143               else:
1144                   Kcl = pcrs*cosd(n*beta)
1145            else:
1146                Kcl = pcrs*(cosd(n*beta)+sind(n*beta))
1147        Fln[i] = SHCoef[term]*Kcl*Lnorm(l)
1148    ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
1149    return ODFln
1150   
1151def polfcal(ODFln,SamSym,psi,gam):
1152    'needs doc string'
1153    import pytexture as ptx
1154    RSQPI = 0.5641895835478
1155    SQ2 = 1.414213562373
1156    PolVal = np.ones_like(gam)
1157    for term in ODFln:
1158        if abs(ODFln[term][1]) > 1.e-3:
1159            l,m,n = eval(term.strip('C'))
1160            psrs,dum = ptx.pyplmpsi(l,m,len(psi),psi)
1161            if SamSym in ['-1','2/m']:
1162                if m != 0:
1163                    Ksl = RSQPI*psrs*(cosd(m*gam)+sind(m*gam))
1164                else:
1165                    Ksl = RSQPI*psrs/SQ2
1166            else:
1167                if m != 0:
1168                    Ksl = RSQPI*psrs*cosd(m*gam)
1169                else:
1170                    Ksl = RSQPI*psrs/SQ2
1171            PolVal += ODFln[term][1]*Ksl
1172    return PolVal
1173   
1174def invpolfcal(ODFln,SGData,phi,beta):
1175    'needs doc string'
1176    import pytexture as ptx
1177   
1178    FORPI = 12.5663706143592
1179    RSQPI = 0.5641895835478
1180    SQ2 = 1.414213562373
1181
1182    invPolVal = np.ones_like(beta)
1183    for term in ODFln:
1184        if abs(ODFln[term][1]) > 1.e-3:
1185            l,m,n = eval(term.strip('C'))
1186            if SGData['SGLaue'] in ['m3','m3m']:
1187                Kcl = 0.0
1188                for j in range(0,l+1,4):
1189                    im = j/4+1
1190                    pcrs,dum = ptx.pyplmpsi(l,j,len(beta),phi)
1191                    Kcl += BOH['L='+str(l)][n-1][im-1]*pcrs*cosd(j*beta)       
1192            else:                #all but cubic
1193                pcrs,dum = ptx.pyplmpsi(l,n,len(beta),phi)
1194                pcrs *= RSQPI
1195                if n == 0:
1196                    pcrs /= SQ2
1197                if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm','R3mR','3m1','31m']:
1198                   if SGData['SGLaue'] in ['3mR','3m1','31m']: 
1199                       if n%6 == 3:
1200                           Kcl = pcrs*sind(n*beta)
1201                       else:
1202                           Kcl = pcrs*cosd(n*beta)
1203                   else:
1204                       Kcl = pcrs*cosd(n*beta)
1205                else:
1206                    Kcl = pcrs*(cosd(n*beta)+sind(n*beta))
1207            invPolVal += ODFln[term][1]*Kcl
1208    return invPolVal
1209   
1210   
1211def textureIndex(SHCoef):
1212    'needs doc string'
1213    Tindx = 1.0
1214    for term in SHCoef:
1215        l = eval(term.strip('C'))[0]
1216        Tindx += SHCoef[term]**2/(2.0*l+1.)
1217    return Tindx
1218   
1219# self-test materials follow.
1220selftestlist = []
1221'''Defines a list of self-tests'''
1222selftestquiet = True
1223def _ReportTest():
1224    'Report name and doc string of current routine when ``selftestquiet`` is False'
1225    if not selftestquiet:
1226        import inspect
1227        caller = inspect.stack()[1][3]
1228        doc = eval(caller).__doc__
1229        if doc is not None:
1230            print('testing '+__file__+' with '+caller+' ('+doc+')')
1231        else:
1232            print('testing '+__file__()+" with "+caller)
1233NeedTestData = True
1234def TestData():
1235    array = np.array
1236    global NeedTestData
1237    NeedTestData = False
1238    global CellTestData
1239    # output from uctbx computed on platform darwin on 2010-05-28
1240    CellTestData = [
1241# cell, g, G, cell*, V, V*
1242  [(4, 4, 4, 90, 90, 90), 
1243   array([[  1.60000000e+01,   9.79717439e-16,   9.79717439e-16],
1244       [  9.79717439e-16,   1.60000000e+01,   9.79717439e-16],
1245       [  9.79717439e-16,   9.79717439e-16,   1.60000000e+01]]), array([[  6.25000000e-02,   3.82702125e-18,   3.82702125e-18],
1246       [  3.82702125e-18,   6.25000000e-02,   3.82702125e-18],
1247       [  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],
1248# cell, g, G, cell*, V, V*
1249  [(4.0999999999999996, 5.2000000000000002, 6.2999999999999998, 100, 80, 130), 
1250   array([[ 16.81      , -13.70423184,   4.48533243],
1251       [-13.70423184,  27.04      ,  -5.6887143 ],
1252       [  4.48533243,  -5.6887143 ,  39.69      ]]), array([[ 0.10206349,  0.05083339, -0.00424823],
1253       [ 0.05083339,  0.06344997,  0.00334956],
1254       [-0.00424823,  0.00334956,  0.02615544]]), (0.31947376387537696, 0.25189277536327803, 0.16172643497798223, 85.283666420376008, 94.716333579624006, 50.825714168082683), 100.98576357983838, 0.0099023858863968445],
1255# cell, g, G, cell*, V, V*
1256  [(3.5, 3.5, 6, 90, 90, 120), 
1257   array([[  1.22500000e+01,  -6.12500000e+00,   1.28587914e-15],
1258       [ -6.12500000e+00,   1.22500000e+01,   1.28587914e-15],
1259       [  1.28587914e-15,   1.28587914e-15,   3.60000000e+01]]), array([[  1.08843537e-01,   5.44217687e-02,   3.36690552e-18],
1260       [  5.44217687e-02,   1.08843537e-01,   3.36690552e-18],
1261       [  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],
1262  ]
1263    global CoordTestData
1264    CoordTestData = [
1265# cell, ((frac, ortho),...)
1266  ((4,4,4,90,90,90,), [
1267 ((0.10000000000000001, 0.0, 0.0),(0.40000000000000002, 0.0, 0.0)),
1268 ((0.0, 0.10000000000000001, 0.0),(2.4492935982947065e-17, 0.40000000000000002, 0.0)),
1269 ((0.0, 0.0, 0.10000000000000001),(2.4492935982947065e-17, -2.4492935982947065e-17, 0.40000000000000002)),
1270 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.40000000000000013, 0.79999999999999993, 1.2)),
1271 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.80000000000000016, 1.2, 0.40000000000000002)),
1272 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(1.2, 0.80000000000000004, 0.40000000000000002)),
1273 ((0.5, 0.5, 0.5),(2.0, 1.9999999999999998, 2.0)),
1274]),
1275# cell, ((frac, ortho),...)
1276  ((4.1,5.2,6.3,100,80,130,), [
1277 ((0.10000000000000001, 0.0, 0.0),(0.40999999999999998, 0.0, 0.0)),
1278 ((0.0, 0.10000000000000001, 0.0),(-0.33424955703700043, 0.39834311042186865, 0.0)),
1279 ((0.0, 0.0, 0.10000000000000001),(0.10939835193016617, -0.051013289294572106, 0.6183281045774256)),
1280 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.069695941716497567, 0.64364635296002093, 1.8549843137322766)),
1281 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(-0.073350319180835066, 1.1440160419710339, 0.6183281045774256)),
1282 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.67089923785616512, 0.74567293154916525, 0.6183281045774256)),
1283 ((0.5, 0.5, 0.5),(0.92574397446582857, 1.7366491056364828, 3.0916405228871278)),
1284]),
1285# cell, ((frac, ortho),...)
1286  ((3.5,3.5,6,90,90,120,), [
1287 ((0.10000000000000001, 0.0, 0.0),(0.35000000000000003, 0.0, 0.0)),
1288 ((0.0, 0.10000000000000001, 0.0),(-0.17499999999999993, 0.3031088913245536, 0.0)),
1289 ((0.0, 0.0, 0.10000000000000001),(3.6739403974420595e-17, -3.6739403974420595e-17, 0.60000000000000009)),
1290 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(2.7675166561703527e-16, 0.60621778264910708, 1.7999999999999998)),
1291 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.17500000000000041, 0.90932667397366063, 0.60000000000000009)),
1292 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.70000000000000018, 0.6062177826491072, 0.60000000000000009)),
1293 ((0.5, 0.5, 0.5),(0.87500000000000067, 1.5155444566227676, 3.0)),
1294]),
1295]
1296    global LaueTestData             #generated by GSAS
1297    LaueTestData = {
1298    '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),
1299        (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),
1300        (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))],
1301    '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),
1302        (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),
1303        (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),
1304        (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))],
1305    '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),
1306        (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),
1307        (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),
1308        (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),
1309        (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),
1310        (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),
1311        (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),
1312        (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),
1313        (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),
1314        (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),
1315        (2,1,5,6),(2,1,-5,6),(3,-1,-5,6))],
1316    '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),
1317        (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),
1318        (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),
1319        (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),
1320        (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),
1321        (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),
1322        (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),
1323        (3,0,4,6),(3,1,-2,12),(3,0,-4,6),(1,1,6,12),(2,2,3,12))],
1324    '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),
1325        (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),
1326        (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),
1327        (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),
1328        (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),
1329        (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),
1330        (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))],
1331    }
1332   
1333    global FLnhTestData
1334    FLnhTestData = [{
1335    'C(4,0,0)': (0.965, 0.42760447),
1336    'C(2,0,0)': (1.0122, -0.80233610),
1337    'C(2,0,2)': (0.0061, 8.37491546E-03),
1338    'C(6,0,4)': (-0.0898, 4.37985696E-02),
1339    'C(6,0,6)': (-0.1369, -9.04081762E-02),
1340    'C(6,0,0)': (0.5935, -0.18234928),
1341    'C(4,0,4)': (0.1872, 0.16358127),
1342    'C(6,0,2)': (0.6193, 0.27573633),
1343    'C(4,0,2)': (-0.1897, 0.12530720)},[1,0,0]]
1344def test0():
1345    if NeedTestData: TestData()
1346    msg = 'test cell2Gmat, fillgmat, Gmat2cell'
1347    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
1348        G, g = cell2Gmat(cell)
1349        assert np.allclose(G,tG),msg
1350        assert np.allclose(g,tg),msg
1351        tcell = Gmat2cell(g)
1352        assert np.allclose(cell,tcell),msg
1353        tcell = Gmat2cell(G)
1354        assert np.allclose(tcell,trcell),msg
1355selftestlist.append(test0)
1356
1357def test1():
1358    'test cell2A and A2Gmat'
1359    _ReportTest()
1360    if NeedTestData: TestData()
1361    msg = 'test cell2A and A2Gmat'
1362    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
1363        G, g = A2Gmat(cell2A(cell))
1364        assert np.allclose(G,tG),msg
1365        assert np.allclose(g,tg),msg
1366selftestlist.append(test1)
1367
1368def test2():
1369    'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
1370    _ReportTest()
1371    if NeedTestData: TestData()
1372    msg = 'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
1373    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
1374        G, g = cell2Gmat(cell)
1375        tcell = A2cell(Gmat2A(G))
1376        assert np.allclose(cell,tcell),msg
1377selftestlist.append(test2)
1378
1379def test3():
1380    'test invcell2Gmat'
1381    _ReportTest()
1382    if NeedTestData: TestData()
1383    msg = 'test invcell2Gmat'
1384    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
1385        G, g = invcell2Gmat(trcell)
1386        assert np.allclose(G,tG),msg
1387        assert np.allclose(g,tg),msg
1388selftestlist.append(test3)
1389
1390def test4():
1391    'test calc_rVsq, calc_rV, calc_V'
1392    _ReportTest()
1393    if NeedTestData: TestData()
1394    msg = 'test calc_rVsq, calc_rV, calc_V'
1395    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
1396        assert np.allclose(calc_rV(cell2A(cell)),trV), msg
1397        assert np.allclose(calc_V(cell2A(cell)),tV), msg
1398selftestlist.append(test4)
1399
1400def test5():
1401    'test A2invcell'
1402    _ReportTest()
1403    if NeedTestData: TestData()
1404    msg = 'test A2invcell'
1405    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
1406        rcell = A2invcell(cell2A(cell))
1407        assert np.allclose(rcell,trcell),msg
1408selftestlist.append(test5)
1409
1410def test6():
1411    'test cell2AB'
1412    _ReportTest()
1413    if NeedTestData: TestData()
1414    msg = 'test cell2AB'
1415    for (cell,coordlist) in CoordTestData:
1416        A,B = cell2AB(cell)
1417        for (frac,ortho) in coordlist:
1418            to = np.inner(A,frac)
1419            tf = np.inner(B,to)
1420            assert np.allclose(ortho,to), msg
1421            assert np.allclose(frac,tf), msg
1422            to = np.sum(A*frac,axis=1)
1423            tf = np.sum(B*to,axis=1)
1424            assert np.allclose(ortho,to), msg
1425            assert np.allclose(frac,tf), msg
1426selftestlist.append(test6)
1427
1428def test7():
1429    'test GetBraviasNum(...) and GenHBravais(...)'
1430    _ReportTest()
1431    import os.path
1432    import sys
1433    import GSASIIspc as spc
1434    testdir = os.path.join(os.path.split(os.path.abspath( __file__ ))[0],'testinp')
1435    if os.path.exists(testdir):
1436        if testdir not in sys.path: sys.path.insert(0,testdir)
1437    import sgtbxlattinp
1438    derror = 1e-4
1439    def indexmatch(hklin, hkllist, system):
1440        for hklref in hkllist:
1441            hklref = list(hklref)
1442            # these permutations are far from complete, but are sufficient to
1443            # allow the test to complete
1444            if system == 'cubic':
1445                permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
1446            elif system == 'monoclinic':
1447                permlist = [(1,2,3),(-1,2,-3)]
1448            else:
1449                permlist = [(1,2,3)]
1450
1451            for perm in permlist:
1452                hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
1453                if hkl == hklref: return True
1454                if [-i for i in hkl] == hklref: return True
1455        else:
1456            return False
1457
1458    for key in sgtbxlattinp.sgtbx7:
1459        spdict = spc.SpcGroup(key)
1460        cell = sgtbxlattinp.sgtbx7[key][0]
1461        system = spdict[1]['SGSys']
1462        center = spdict[1]['SGLatt']
1463
1464        bravcode = GetBraviasNum(center, system)
1465
1466        g2list = GenHBravais(sgtbxlattinp.dmin, bravcode, cell2A(cell))
1467
1468        assert len(sgtbxlattinp.sgtbx7[key][1]) == len(g2list), 'Reflection lists differ for %s' % key
1469        for h,k,l,d,num in g2list:
1470            for hkllist,dref in sgtbxlattinp.sgtbx7[key][1]: 
1471                if abs(d-dref) < derror:
1472                    if indexmatch((h,k,l,), hkllist, system):
1473                        break
1474            else:
1475                assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
1476selftestlist.append(test7)
1477
1478def test8():
1479    'test GenHLaue'
1480    _ReportTest()
1481    import GSASIIspc as spc
1482    import sgtbxlattinp
1483    derror = 1e-4
1484    dmin = sgtbxlattinp.dmin
1485
1486    def indexmatch(hklin, hklref, system, axis):
1487        # these permutations are far from complete, but are sufficient to
1488        # allow the test to complete
1489        if system == 'cubic':
1490            permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
1491        elif system == 'monoclinic' and axis=='b':
1492            permlist = [(1,2,3),(-1,2,-3)]
1493        elif system == 'monoclinic' and axis=='a':
1494            permlist = [(1,2,3),(1,-2,-3)]
1495        elif system == 'monoclinic' and axis=='c':
1496            permlist = [(1,2,3),(-1,-2,3)]
1497        elif system == 'trigonal':
1498            permlist = [(1,2,3),(2,1,3),(-1,-2,3),(-2,-1,3)]
1499        elif system == 'rhombohedral':
1500            permlist = [(1,2,3),(2,3,1),(3,1,2)]
1501        else:
1502            permlist = [(1,2,3)]
1503
1504        hklref = list(hklref)
1505        for perm in permlist:
1506            hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
1507            if hkl == hklref: return True
1508            if [-i for i in hkl] == hklref: return True
1509        return False
1510
1511    for key in sgtbxlattinp.sgtbx8:
1512        spdict = spc.SpcGroup(key)[1]
1513        cell = sgtbxlattinp.sgtbx8[key][0]
1514        center = spdict['SGLatt']
1515        Laue = spdict['SGLaue']
1516        Axis = spdict['SGUniq']
1517        system = spdict['SGSys']
1518
1519        g2list = GenHLaue(dmin,spdict,cell2A(cell))
1520        #if len(g2list) != len(sgtbxlattinp.sgtbx8[key][1]):
1521        #    print 'failed',key,':' ,len(g2list),'vs',len(sgtbxlattinp.sgtbx8[key][1])
1522        #    print 'GSAS-II:'
1523        #    for h,k,l,d in g2list: print '  ',(h,k,l),d
1524        #    print 'SGTBX:'
1525        #    for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print '  ',hkllist,dref
1526        assert len(g2list) == len(sgtbxlattinp.sgtbx8[key][1]), (
1527            'Reflection lists differ for %s' % key
1528            )
1529        #match = True
1530        for h,k,l,d in g2list:
1531            for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: 
1532                if abs(d-dref) < derror:
1533                    if indexmatch((h,k,l,), hkllist, system, Axis): break
1534            else:
1535                assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
1536                #match = False
1537        #if not match:
1538            #for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print '  ',hkllist,dref
1539            #print center, Laue, Axis, system
1540selftestlist.append(test8)
1541           
1542def test9():
1543    'test GenHLaue'
1544    _ReportTest()
1545    import GSASIIspc as G2spc
1546    if NeedTestData: TestData()
1547    for spc in LaueTestData:
1548        data = LaueTestData[spc]
1549        cell = data[0]
1550        hklm = np.array(data[1])
1551        H = hklm[-1][:3]
1552        hklO = hklm.T[:3].T
1553        A = cell2A(cell)
1554        dmin = 1./np.sqrt(calc_rDsq(H,A))
1555        SGData = G2spc.SpcGroup(spc)[1]
1556        hkls = np.array(GenHLaue(dmin,SGData,A))
1557        hklN = hkls.T[:3].T
1558        #print spc,hklO.shape,hklN.shape
1559        err = True
1560        for H in hklO:
1561            if H not in hklN:
1562                print H,' missing from hkl from GSASII'
1563                err = False
1564        assert(err)
1565selftestlist.append(test9)
1566       
1567       
1568   
1569
1570if __name__ == '__main__':
1571    # run self-tests
1572    selftestquiet = False
1573    for test in selftestlist:
1574        test()
1575    print "OK"
Note: See TracBrowser for help on using the repository browser.