source: trunk/GSASIIlattice.py @ 959

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

added cell2GS & Uij2Ueqv

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