source: trunk/GSASIIlattice.py @ 303

Last change on this file since 303 was 303, checked in by vondreele, 12 years ago

Allow reading of xye (Topas style) files
Begin implementation of spherical harmonics texture
Refactor indexing; replace cell refinement from peak positions
  • Property svn:keywords set to Date Author Revision URL Id
File size: 40.9 KB
Line 
1'''Perform lattice-related computations'''
2
3########### SVN repository information ###################
4# $Date: 2011-06-20 19:37:07 +0000 (Mon, 20 Jun 2011) $
5# $Author: vondreele $
6# $Revision: 303 $
7# $URL: trunk/GSASIIlattice.py $
8# $Id: GSASIIlattice.py 303 2011-06-20 19:37:07Z vondreele $
9########### SVN repository information ###################
10import math
11import numpy as np
12import numpy.linalg as nl
13
14# trig functions in degrees
15sind = lambda x: np.sin(x*np.pi/180.)
16asind = lambda x: 180.*np.arcsin(x)/np.pi
17tand = lambda x: np.tan(x*np.pi/180.)
18atand = lambda x: 180.*np.arctan(x)/np.pi
19atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
20cosd = lambda x: np.cos(x*np.pi/180.)
21acosd = lambda x: 180.*np.arccos(x)/np.pi
22rdsq2d = lambda x,p: round(1.0/np.sqrt(x),p)
23
24def sec2HMS(sec):
25    H = int(sec/3600)
26    M = int(sec/60-H*60)
27    S = sec-3600*H-60*M
28    return '%d:%2d:%.2f'%(H,M,S)
29   
30def rotdMat(angle,axis=0):
31    '''Prepare rotation matrix for angle in degrees about axis(=0,1,2)
32    Returns numpy 3,3 array
33    '''
34    if axis == 2:
35        return np.array([[cosd(angle),-sind(angle),0],[sind(angle),cosd(angle),0],[0,0,1]])
36    elif axis == 1:
37        return np.array([[cosd(angle),0,-sind(angle)],[0,1,0],[sind(angle),0,cosd(angle)]])
38    else:
39        return np.array([[1,0,0],[0,cosd(angle),-sind(angle)],[0,sind(angle),cosd(angle)]])
40       
41def rotdMat4(angle,axis=0):
42    '''Prepare rotation matrix for angle in degrees about axis(=0,1,2) with scaling for OpenGL
43    Returns numpy 4,4 array
44    '''
45    Mat = rotdMat(angle,axis)
46    return np.concatenate((np.concatenate((Mat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
47   
48def fillgmat(cell):
49    '''Compute lattice metric tensor from unit cell constants
50    cell is tuple with a,b,c,alpha, beta, gamma (degrees)
51    returns 3x3 numpy array
52    '''
53    a,b,c,alp,bet,gam = cell
54    g = np.array([
55        [a*a,  a*b*cosd(gam),  a*c*cosd(bet)],
56        [a*b*cosd(gam),  b*b,  b*c*cosd(alp)],
57        [a*c*cosd(bet) ,b*c*cosd(alp),   c*c]])
58    return g
59           
60def cell2Gmat(cell):
61    '''Compute real and reciprocal lattice metric tensor from unit cell constants
62    cell is tuple with a,b,c,alpha, beta, gamma (degrees)
63    returns reciprocal (G) & real (g) metric tensors (list of two 3x3 arrays)
64    '''
65    g = fillgmat(cell)
66    G = nl.inv(g)       
67    return G,g
68
69def A2Gmat(A):
70    '''Fill reciprocal metric tensor (G) from A
71    returns reciprocal (G) & real (g) metric tensors (list of two 3x3 arrays)
72    '''
73    G = np.zeros(shape=(3,3))
74    G = [
75        [A[0],  A[3]/2.,  A[4]/2.], 
76        [A[3]/2.,A[1],    A[5]/2.], 
77        [A[4]/2.,A[5]/2.,    A[2]]]
78    g = nl.inv(G)
79    return G,g
80
81def Gmat2A(G):
82    'Extract A from reciprocal metric tensor (G)'
83    return [G[0][0],G[1][1],G[2][2],2.*G[0][1],2.*G[0][2],2.*G[1][2]]
84   
85def cell2A(cell):
86    G,g = cell2Gmat(cell)
87    return Gmat2A(G)
88
89def A2cell(A):
90    '''Compute unit cell constants from A tensor
91    returns tuple with a,b,c,alpha, beta, gamma (degrees)
92    '''
93    G,g = A2Gmat(A)
94    return Gmat2cell(g)
95
96def Gmat2cell(g):
97    '''Compute lattice parameters from real metric tensor (g)
98    returns tuple with a,b,c,alpha, beta, gamma (degrees)
99    Alternatively,compute reciprocal lattice parameters from inverse metric tensor (G)
100    returns tuple with a*,b*,c*,alpha*, beta*, gamma* (degrees)
101    '''
102    a = np.sqrt(max(0,g[0][0]))
103    b = np.sqrt(max(0,g[1][1]))
104    c = np.sqrt(max(0,g[2][2]))
105    alp = acosd(g[2][1]/(b*c))
106    bet = acosd(g[2][0]/(a*c))
107    gam = acosd(g[0][1]/(a*b))
108    return a,b,c,alp,bet,gam
109
110def invcell2Gmat(invcell):
111    '''Compute real and reciprocal lattice metric tensor from reciprocal
112       unit cell constants
113    invcell is tuple with a*,b*,c*,alpha*, beta*, gamma* (degrees)
114    returns reciprocal (G) & real (g) metric tensors (list of two 3x3 arrays)
115    '''
116    G = fillgmat(invcell)
117    g = nl.inv(G)
118    return G,g
119       
120def calc_rVsq(A):
121    'Compute the square of the reciprocal lattice volume (V* **2) from A'
122    G,g = A2Gmat(A)
123    rVsq = nl.det(G)
124    if rVsq < 0:
125        return 1
126    return rVsq
127   
128def calc_rV(A):
129    'Compute the reciprocal lattice volume (V*) from A'
130    return np.sqrt(calc_rVsq(A))
131   
132def calc_V(A):
133    'Compute the real lattice volume (V) from A'
134    return 1./calc_rV(A)
135
136def A2invcell(A):
137    '''Compute reciprocal unit cell constants from A
138    returns tuple with a*,b*,c*,alpha*, beta*, gamma* (degrees)
139    '''
140    G,g = A2Gmat(A)
141    return Gmat2cell(G)
142
143def cell2AB(cell):
144    '''Computes orthogonalization matrix from unit cell constants
145    cell is tuple with a,b,c,alpha, beta, gamma (degrees)
146    returns tuple of two 3x3 numpy arrays (A,B)
147       A for crystal to Cartesian transformations A*x = np.inner(A,x) = X
148       B (= inverse of A) for Cartesian to crystal transformation B*X = np.inner(B*x) = x
149    '''
150    G,g = cell2Gmat(cell) 
151    cellstar = Gmat2cell(G)
152    A = np.zeros(shape=(3,3))
153    # from Giacovazzo (Fundamentals 2nd Ed.) p.75
154    A[0][0] = cell[0]                # a
155    A[0][1] = cell[1]*cosd(cell[5])  # b cos(gamma)
156    A[0][2] = cell[2]*cosd(cell[4])  # c cos(beta)
157    A[1][1] = cell[1]*sind(cell[5])  # b sin(gamma)
158    A[1][2] = -cell[2]*cosd(cellstar[3])*sind(cell[4]) # - c cos(alpha*) sin(beta)
159    A[2][2] = 1/cellstar[2]         # 1/c*
160    B = nl.inv(A)
161    return A,B
162       
163def Uij2betaij(Uij,G):
164    '''
165    Convert Uij to beta-ij tensors
166    input:
167    Uij - numpy array [Uij]
168    G - reciprocal metric tensor
169    returns:
170    beta-ij - numpy array [beta-ij]
171    '''
172    pass
173   
174def criticalEllipse(prob):
175    '''
176    Calculate critical values for probability ellipsoids from probability
177    '''
178    if not ( 0.01 <= prob < 1.0):
179        return 1.54 
180    coeff = np.array([6.44988E-09,4.16479E-07,1.11172E-05,1.58767E-04,0.00130554,
181        0.00604091,0.0114921,-0.040301,-0.6337203,1.311582])
182    llpr = math.log(-math.log(prob))
183    return np.polyval(coeff,llpr)
184   
185def CellBlock(nCells):
186    '''
187    Generate block of unit cells n*n*n on a side; [0,0,0] centered, n = 2*nCells+1
188    currently only works for nCells = 0 or 1 (not >1)
189    '''
190    if nCells:
191        N = 2*nCells+1
192        N2 = N*N
193        N3 = N*N*N
194        cellArray = []
195        A = np.array(range(N3))
196        cellGen = np.array([A/N2-1,A/N%N-1,A%N-1]).T
197        for cell in cellGen:
198            cellArray.append(cell)
199        return cellArray
200    else:
201        return [0,0,0]
202       
203def CellAbsorption(ElList,Volume):
204# ElList = dictionary of element contents including mu
205    muT = 0
206    for El in ElList:
207        muT += ElList[El]['mu']*ElList[El]['FormulaNo']
208    return muT/Volume
209   
210#Permutations and Combinations
211# Four routines: combinations,uniqueCombinations, selections & permutations
212#These taken from Python Cookbook, 2nd Edition. 19.15 p724-726
213#   
214def _combinators(_handle, items, n):
215    ''' factored-out common structure of all following combinators '''
216    if n==0:
217        yield [ ]
218        return
219    for i, item in enumerate(items):
220        this_one = [ item ]
221        for cc in _combinators(_handle, _handle(items, i), n-1):
222            yield this_one + cc
223def combinations(items, n):
224    ''' take n distinct items, order matters '''
225    def skipIthItem(items, i):
226        return items[:i] + items[i+1:]
227    return _combinators(skipIthItem, items, n)
228def uniqueCombinations(items, n):
229    ''' take n distinct items, order is irrelevant '''
230    def afterIthItem(items, i):
231        return items[i+1:]
232    return _combinators(afterIthItem, items, n)
233def selections(items, n):
234    ''' take n (not necessarily distinct) items, order matters '''
235    def keepAllItems(items, i):
236        return items
237    return _combinators(keepAllItems, items, n)
238def permutations(items):
239    ''' take all items, order matters '''
240    return combinations(items, len(items))
241
242#reflection generation routines
243#for these: H = [h,k,l]; A is as used in calc_rDsq; G - inv metric tensor, g - metric tensor;
244#           cell - a,b,c,alp,bet,gam in A & deg
245   
246def calc_rDsq(H,A):
247    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]
248    return rdsq
249   
250def calc_rDsq2(H,G):
251    return np.inner(H,np.inner(G,H))
252   
253def calc_rDsqZ(H,A,Z,tth,lam):
254    rpd = math.pi/180.
255    rdsq = calc_rDsq(H,A)+Z*math.sin(tth*rpd)*2.0*rpd/(lam*lam)
256    return rdsq
257       
258def MaxIndex(dmin,A):
259    Hmax = [0,0,0]
260    try:
261        cell = A2cell(A)
262    except:
263        cell = [1,1,1,90,90,90]
264    for i in range(3):
265        Hmax[i] = int(round(cell[i]/dmin))
266    return Hmax
267   
268def sortHKLd(HKLd,ifreverse,ifdup):
269    #HKLd is a list of [h,k,l,d,...]; ifreverse=True for largest d first
270    #ifdup = True if duplicate d-spacings allowed
271    T = []
272    for i,H in enumerate(HKLd):
273        if ifdup:
274            T.append((H[3],i))
275        else:
276            T.append(H[3])           
277    D = dict(zip(T,HKLd))
278    T.sort()
279    if ifreverse:
280        T.reverse()
281    X = []
282    okey = ''
283    for key in T: 
284        if key != okey: X.append(D[key])    #remove duplicate d-spacings
285        okey = key
286    return X
287   
288def SwapIndx(Axis,H):
289    if Axis in [1,-1]:
290        return H
291    elif Axis in [2,-3]:
292        return [H[1],H[2],H[0]]
293    else:
294        return [H[2],H[0],H[1]]
295       
296def Rh2Hx(Rh):
297    Hx = [0,0,0]
298    Hx[0] = Rh[0]-Rh[1]
299    Hx[1] = Rh[1]-Rh[2]
300    Hx[2] = np.sum(Rh)
301    return Hx
302   
303def Hx2Rh(Hx):
304        Rh = [0,0,0]
305        itk = -Hx[0]+Hx[1]+Hx[2]
306        if itk%3 != 0:
307            return 0        #error - not rhombohedral reflection
308        else:
309            Rh[1] = itk/3
310            Rh[0] = Rh[1]+Hx[0]
311            Rh[2] = Rh[1]-Hx[1]
312            if Rh[0] < 0:
313                for i in range(3):
314                    Rh[i] = -Rh[i]
315            return Rh
316       
317def CentCheck(Cent,H):
318    h,k,l = H
319    if Cent == 'A' and (k+l)%2:
320        return False
321    elif Cent == 'B' and (h+l)%2:
322        return False
323    elif Cent == 'C' and (h+k)%2:
324        return False
325    elif Cent == 'I' and (h+k+l)%2:
326        return False
327    elif Cent == 'F' and ((h+k)%2 or (h+l)%2 or (k+l)%2):
328        return False
329    elif Cent == 'R' and (-h+k+l)%3:
330        return False
331    else:
332        return True
333                                   
334def GetBraviasNum(center,system):
335    '''Determine the Bravais lattice number, as used in GenHBravais
336         center = one of: P, C, I, F, R (see SGLatt from GSASIIspc.SpcGroup)
337         lattice = is cubic, hexagonal, tetragonal, orthorhombic, trigonal (R)
338             monoclinic, triclinic (see SGSys from GSASIIspc.SpcGroup)
339       Returns a number between 0 and 13
340          or throws an exception if the setting is non-standard
341       '''
342    if center.upper() == 'F' and system.lower() == 'cubic':
343        return 0
344    elif center.upper() == 'I' and system.lower() == 'cubic':
345        return 1
346    elif center.upper() == 'P' and system.lower() == 'cubic':
347        return 2
348    elif center.upper() == 'R' and system.lower() == 'trigonal':
349        return 3
350    elif center.upper() == 'P' and system.lower() == 'hexagonal':
351        return 4
352    elif center.upper() == 'I' and system.lower() == 'tetragonal':
353        return 5
354    elif center.upper() == 'P' and system.lower() == 'tetragonal':
355        return 6
356    elif center.upper() == 'F' and system.lower() == 'orthorhombic':
357        return 7
358    elif center.upper() == 'I' and system.lower() == 'orthorhombic':
359        return 8
360    elif center.upper() == 'C' and system.lower() == 'orthorhombic':
361        return 9
362    elif center.upper() == 'P' and system.lower() == 'orthorhombic':
363        return 10
364    elif center.upper() == 'C' and system.lower() == 'monoclinic':
365        return 11
366    elif center.upper() == 'P' and system.lower() == 'monoclinic':
367        return 12
368    elif center.upper() == 'P' and system.lower() == 'triclinic':
369        return 13
370    raise ValueError,'non-standard Bravais lattice center=%s, cell=%s' % (center,system)
371
372def GenHBravais(dmin,Bravais,A):
373    '''Generate the positionally unique powder diffraction reflections
374    input:
375       dmin is minimum d-space
376       Bravais is 0-13 to indicate lattice type (see GetBraviasNum)
377       A is reciprocal cell tensor (see Gmat2A or cell2A)
378    returns:
379       a list of tuples containing: h,k,l,d-space,-1   
380    '''
381# Bravais in range(14) to indicate Bravais lattice:
382#   0 F cubic
383#   1 I cubic
384#   2 P cubic
385#   3 R hexagonal (trigonal not rhombohedral)
386#   4 P hexagonal
387#   5 I tetragonal
388#   6 P tetragonal
389#   7 F orthorhombic
390#   8 I orthorhombic
391#   9 C orthorhombic
392#  10 P orthorhombic
393#  11 C monoclinic
394#  12 P monoclinic
395#  13 P triclinic
396# A - as defined in calc_rDsq
397# returns HKL = [h,k,l,d,0] sorted so d largest first
398    import math
399    if Bravais in [9,11]:
400        Cent = 'C'
401    elif Bravais in [1,5,8]:
402        Cent = 'I'
403    elif Bravais in [0,7]:
404        Cent = 'F'
405    elif Bravais in [3]:
406        Cent = 'R'
407    else:
408        Cent = 'P'
409    Hmax = MaxIndex(dmin,A)
410    dminsq = 1./(dmin**2)
411    HKL = []
412    if Bravais == 13:                       #triclinic
413        for l in range(-Hmax[2],Hmax[2]+1):
414            for k in range(-Hmax[1],Hmax[1]+1):
415                hmin = 0
416                if (k < 0): hmin = 1
417                if (k ==0 and l < 0): hmin = 1
418                for h in range(hmin,Hmax[0]+1):
419                    H=[h,k,l]
420                    rdsq = calc_rDsq(H,A)
421                    if 0 < rdsq <= dminsq:
422                        HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
423    elif Bravais in [11,12]:                #monoclinic - b unique
424        Hmax = SwapIndx(2,Hmax)
425        for h in range(Hmax[0]+1):
426            for k in range(-Hmax[1],Hmax[1]+1):
427                lmin = 0
428                if k < 0:lmin = 1
429                for l in range(lmin,Hmax[2]+1):
430                    [h,k,l] = SwapIndx(-2,[h,k,l])
431                    H = []
432                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
433                    if H:
434                        rdsq = calc_rDsq(H,A)
435                        if 0 < rdsq <= dminsq:
436                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
437                    [h,k,l] = SwapIndx(2,[h,k,l])
438    elif Bravais in [7,8,9,10]:            #orthorhombic
439        for h in range(Hmax[0]+1):
440            for k in range(Hmax[1]+1):
441                for l in range(Hmax[2]+1):
442                    H = []
443                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
444                    if H:
445                        rdsq = calc_rDsq(H,A)
446                        if 0 < rdsq <= dminsq:
447                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
448    elif Bravais in [5,6]:                  #tetragonal
449        for l in range(Hmax[2]+1):
450            for k in range(Hmax[1]+1):
451                for h in range(k,Hmax[0]+1):
452                    H = []
453                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
454                    if H:
455                        rdsq = calc_rDsq(H,A)
456                        if 0 < rdsq <= dminsq:
457                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
458    elif Bravais in [3,4]:
459        lmin = 0
460        if Bravais == 3: lmin = -Hmax[2]                  #hexagonal/trigonal
461        for l in range(lmin,Hmax[2]+1):
462            for k in range(Hmax[1]+1):
463                hmin = k
464                if l < 0: hmin += 1
465                for h in range(hmin,Hmax[0]+1):
466                    H = []
467                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
468                    if H:
469                        rdsq = calc_rDsq(H,A)
470                        if 0 < rdsq <= dminsq:
471                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
472
473    else:                                   #cubic
474        for l in range(Hmax[2]+1):
475            for k in range(l,Hmax[1]+1):
476                for h in range(k,Hmax[0]+1):
477                    H = []
478                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
479                    if H:
480                        rdsq = calc_rDsq(H,A)
481                        if 0 < rdsq <= dminsq:
482                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
483    return sortHKLd(HKL,True,False)
484   
485def GenHLaue(dmin,SGLaue,SGLatt,SGUniq,A):
486    '''Generate the crystallographically unique powder diffraction reflections
487    for a lattice and Bravais type
488    '''
489# dmin - minimum d-spacing
490# SGLaue - Laue group symbol = '-1','2/m','mmm','4/m','6/m','4/mmm','6/mmm',
491#                            '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
492# SGLatt - lattice centering = 'P','A','B','C','I','F'
493# SGUniq - code for unique monoclinic axis = 'a','b','c'
494# A - 6 terms as defined in calc_rDsq
495# returns - HKL = list of [h,k,l,d] sorted with largest d first and is unique
496# part of reciprocal space ignoring anomalous dispersion
497    import math
498    #finds maximum allowed hkl for given A within dmin
499    if SGLaue in ['3R','3mR']:        #Rhombohedral axes
500        Hmax = [0,0,0]
501        cell = A2cell(A)
502        aHx = cell[0]*math.sqrt(2.0*(1.0-cosd(cell[3])))
503        cHx = cell[0]*math.sqrt(3.0*(1.0+2.0*cosd(cell[3])))
504        Hmax[0] = Hmax[1] = int(round(aHx/dmin))
505        Hmax[2] = int(round(cHx/dmin))
506        #print Hmax,aHx,cHx
507    else:                           # all others
508        Hmax = MaxIndex(dmin,A)
509       
510    dminsq = 1./(dmin**2)
511    HKL = []
512    if SGLaue == '-1':                       #triclinic
513        for l in range(-Hmax[2],Hmax[2]+1):
514            for k in range(-Hmax[1],Hmax[1]+1):
515                hmin = 0
516                if (k < 0) or (k ==0 and l < 0): hmin = 1
517                for h in range(hmin,Hmax[0]+1):
518                    H = []
519                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
520                    rdsq = calc_rDsq(H,A)
521                    if 0 < rdsq <= dminsq:
522                        HKL.append([h,k,l,1/math.sqrt(rdsq)])
523    elif SGLaue == '2/m':                #monoclinic
524        axisnum = 1 + ['a','b','c'].index(SGUniq)
525        Hmax = SwapIndx(axisnum,Hmax)
526        for h in range(Hmax[0]+1):
527            for k in range(-Hmax[1],Hmax[1]+1):
528                lmin = 0
529                if k < 0:lmin = 1
530                for l in range(lmin,Hmax[2]+1):
531                    [h,k,l] = SwapIndx(-axisnum,[h,k,l])
532                    H = []
533                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
534                    if H:
535                        rdsq = calc_rDsq(H,A)
536                        if 0 < rdsq <= dminsq:
537                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
538                    [h,k,l] = SwapIndx(axisnum,[h,k,l])
539    elif SGLaue in ['mmm','4/m','6/m']:            #orthorhombic
540        for l in range(Hmax[2]+1):
541            for h in range(Hmax[0]+1):
542                kmin = 1
543                if SGLaue == 'mmm' or h ==0: kmin = 0
544                for k in range(kmin,Hmax[1]+1):
545                    H = []
546                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
547                    if H:
548                        rdsq = calc_rDsq(H,A)
549                        if 0 < rdsq <= dminsq:
550                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
551    elif SGLaue in ['4/mmm','6/mmm']:                  #tetragonal & hexagonal
552        for l in range(Hmax[2]+1):
553            for h in range(Hmax[0]+1):
554                for k in range(h+1):
555                    H = []
556                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
557                    if H:
558                        rdsq = calc_rDsq(H,A)
559                        if 0 < rdsq <= dminsq:
560                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
561    elif SGLaue in ['3m1','31m','3','3R','3mR']:                  #trigonals
562        for l in range(-Hmax[2],Hmax[2]+1):
563            hmin = 0
564            if l < 0: hmin = 1
565            for h in range(hmin,Hmax[0]+1):
566                if SGLaue in ['3R','3']:
567                    kmax = h
568                    kmin = -int((h-1.)/2.)
569                else:
570                    kmin = 0
571                    kmax = h
572                    if SGLaue in ['3m1','3mR'] and l < 0: kmax = h-1
573                    if SGLaue == '31m' and l < 0: kmin = 1
574                for k in range(kmin,kmax+1):
575                    H = []
576                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
577                    if SGLaue in ['3R','3mR']:
578                        H = Hx2Rh(H)
579                    if H:
580                        rdsq = calc_rDsq(H,A)
581                        if 0 < rdsq <= dminsq:
582                            HKL.append([H[0],H[1],H[2],1/math.sqrt(rdsq)])
583    else:                                   #cubic
584        for h in range(Hmax[0]+1):
585            for k in range(h+1):
586                lmin = 0
587                lmax = k
588                if SGLaue =='m3':
589                    lmax = h-1
590                    if h == k: lmax += 1
591                for l in range(lmin,lmax+1):
592                    H = []
593                    if CentCheck(SGLatt,[h,k,l]): H=[h,k,l]
594                    if H:
595                        rdsq = calc_rDsq(H,A)
596                        if 0 < rdsq <= dminsq:
597                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
598    return sortHKLd(HKL,True,True)
599
600#Spherical harmonics routines
601def OdfChk(SGLaue,L,M):
602    if not L%2 and abs(M) <= L:
603        if SGLaue == '0':                      #cylindrical symmetry
604            if M == 0: return True
605        elif SGLaue == '-1':
606            return True
607        elif SGLaue == '2/m':
608            if not abs(M)%2: return True
609        elif SGLaue == 'mmm':
610            if not abs(M)%2 and M >= 0: return True
611        elif SGLaue == '4/m':
612            if not abs(M)%4: return True
613        elif SGLaue == '4/mmm':
614            if not abs(M)%4 and M >= 0: return True
615        elif SGLaue in ['3R','3']:
616            if not abs(M)%3: return True
617        elif SGLaue in ['3mR','3m1','31m']:
618            if not abs(M)%3 and M >= 0: return True
619        elif SGLaue == '6/m':
620            if not abs(M)%6: return True
621        elif SGLaue == '6/mmm':
622            if not abs(M)%6 and M >= 0: return True
623        elif SGLaue == 'm3':
624            if M > 0:
625                if L%12 == 2:
626                    if M <= L/12: return True
627                else:
628                    if M <= L/12+1: return True
629        elif SGLaue == 'm3m':
630            if M > 0:
631                if L%12 == 2:
632                    if M <= L/12: return True
633                else:
634                    if M <= L/12+1: return True
635    return False
636       
637def GenSHCoeff(SGLaue,SamSym,L):
638    coeffNames = []
639    for iord in [2*i+2 for i in range(L/2)]:
640        for m in [i-iord for i in range(2*iord+1)]:
641            if OdfChk(SamSym,iord,m):
642                for n in [i-iord for i in range(2*iord+1)]:
643                    if OdfChk(SGLaue,iord,n):
644                        coeffNames.append('C(%d,%d,%d)'%(iord,m,n))
645    return coeffNames
646
647def CrsAng(H,cell,SGData):
648    a,b,c,al,be,ga = cell
649    SQ3 = 1.732050807569
650    H1 = np.array([1,0,0])
651    H2 = np.array([0,1,0])
652    H3 = np.array([0,0,1])
653    H4 = np.array([1,1,1])
654    G,g = cell2Gmat(cell)
655    Laue = SGData['SGLaue']
656    Naxis = SGData['SGUniq']
657    DH = np.inner(H,np.inner(G,H))
658    if Laue == '2/m':
659        if Naxis == 'a':
660            DR = np.inner(H1,np.inner(G,H1))
661            DHR = np.inner(H,np.inner(G,H1))
662        elif Naxis == 'b':
663            DR = np.inner(H2,np.inner(G,H2))
664            DHR = np.inner(H,np.inner(G,H2))
665        else:
666            DR = np.inner(H3,np.inner(G,H3))
667            DHR = np.inner(H,np.inner(G,H3))
668    elif Laue in ['R3','R3m']:
669        DR = np.inner(H4,np.inner(G,H4))
670        DRH = np.inner(H,np.inner(G,H4))
671       
672    else:
673        DR = np.inner(H3,np.inner(G,H3))
674        DHR = np.inner(H,np.inner(G,H3))
675    DHR /= np.sqrt(DR*DH)
676    phi = acosd(DHR)
677    if Laue == '-1':
678        BA = H[1]*a/(b-H[0]*cosd(ga))
679        BB = H[0]*sind(ga)**2
680    elif Laue == '2/m':
681        if Naxis == 'a':
682            BA = H[2]*b/(c-H[1]*cods(al))
683            BB = H[1]*sind(al)**2
684        elif Naxis == 'b':
685            BA = H[0]*c/(a-H[2]*cosd(be))
686            BB = H[2]*sind(be)**2
687        else:
688            BA = H[1]*a/(b-H[0]*cosd(ga))
689            BB = H[0]*sind(ga)**2
690    elif Laue in ['mmm','4/m','4/mmm']:
691        BA = H[1]*a
692        BB = H[0]*b
693   
694    elif Laue in ['3R','3mR']:
695        BA = H[0]+H[1]-2.0*H[2]
696        BB = SQ3*(H[0]-H[1])
697    elif Laue in ['m3','m3m']:
698        BA = H[1]
699        BB = H[0]
700    else:
701        BA = H[0]+2.0*H[1]
702        BB = SQ3*H[0]
703    beta = atan2d(BA,BB)
704    return phi,beta
705   
706def SamAng(Tth,Gangls,Sangl,IFCoup):
707    if IFCoup:
708        GSomeg = sind(Gangls[2]+Tth)
709        GComeg = cosd(Gangls[2]+Tth)
710    else:
711        GSomeg = sind(Gangls[2])
712        GComeg = cosd(Gangls[2])
713    GSTth = sind(Tth)
714    GCTth = cosd(Tth)     
715    GSazm = sind(Gangls[3])
716    GCazm = cosd(Gangls[3])
717    GSchi = sind(Gangls[1])
718    GCchi = cosd(Gangls[1])
719    GSphi = sind(Gangls[0]+Sangl[2])
720    GCphi = cosd(Gangls[0]+Sangl[2])
721    SSomeg = sind(Sangl[0])
722    SComeg = cosd(Sangl[0])
723    SSchi = sind(Sangl[1])
724    SCchi = cosd(Sangl[1])
725    AT = -GSTth*GComeg+GCTth*GCazm*GSomeg
726    BT = GSTth*GSomeg+GCTth*GCazm*GComeg
727    CT = -GCTth*GSazm*GSchi
728    DT = -GCTth*GSazm*GCchi
729   
730    BC1 = -AT*GSphi+(CT+BT*GCchi)*GCphi
731    BC2 = DT-BT*GSchi
732    BC3 = AT*GCphi+(CT+BT*GCchi)*GSphi
733     
734    BC = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg     
735    psi = acosd(BC)
736     
737    BA = -BC1*SSchi+BC2*SCchi
738    BB = BC1*SSomeg*SCchi+BC2*SSomeg*SSchi+BC3*SComeg
739    gam = atand2(BB,BA)
740       
741    return psi,gam
742   
743def Flnh(Start,SHCoef,phi,beta,SGData):
744    import pytexture as ptx
745    BOH = {
746    'L=2':[[],[],[]],
747    'L=4':[[0.30469720,0.36418281],[],[]],
748    'L=6':[[-0.14104740,0.52775103],[],[]],
749    'L=8':[[0.28646862,0.21545346,0.32826995],[],[]],
750    'L=10':[[-0.16413497,0.33078546,0.39371345],[],[]],
751    'L=12':[[0.26141975,0.27266871,0.03277460,0.32589402],
752        [0.09298802,-0.23773812,0.49446631],[]],
753    'L=14':[[-0.17557309,0.25821932,0.27709173,0.33645360],[],[]],
754    'L=16':[[0.24370673,0.29873515,0.06447688,0.00377,0.32574495],
755        [0.12039646,-0.25330128,0.23950998,0.40962508],[]],
756    'L=18':[[-0.16914245,0.17017340,0.34598142,0.07433932,0.32696037],
757        [-0.06901768,0.16006562,-0.24743528,0.47110273],[]],
758    'L=20':[[0.23067026,0.31151832,0.09287682,0.01089683,0.00037564,0.32573563],
759        [0.13615420,-0.25048007,0.12882081,0.28642879,0.34620433],[]],
760    'L=22':[[-0.16109560,0.10244188,0.36285175,0.13377513,0.01314399,0.32585583],
761        [-0.09620055,0.20244115,-0.22389483,0.17928946,0.42017231],[]],
762    'L=24':[[0.22050742,0.31770654,0.11661736,0.02049853,0.00150861,0.00003426,0.32573505],
763        [0.13651722,-0.21386648,0.00522051,0.33939435,0.10837396,0.32914497],
764        [0.05378596,-0.11945819,0.16272298,-0.26449730,0.44923956]],
765    'L=26':[[-0.15435003,0.05261630,0.35524646,0.18578869,0.03259103,0.00186197,0.32574594],
766        [-0.11306511,0.22072681,-0.18706142,0.05439948,0.28122966,0.35634355],[]],
767    'L=28':[[0.21225019,0.32031716,0.13604702,0.03132468,0.00362703,0.00018294,0.00000294,0.32573501],
768        [0.13219496,-0.17206256,-0.08742608,0.32671661,0.17973107,0.02567515,0.32619598],
769        [0.07989184,-0.16735346,0.18839770,-0.20705337,0.12926808,0.42715602]],
770    'L=30':[[-0.14878368,0.01524973,0.33628434,0.22632587,0.05790047,0.00609812,0.00022898,0.32573594],
771        [-0.11721726,0.20915005,-0.11723436,-0.07815329,0.31318947,0.13655742,0.33241385],
772        [-0.04297703,0.09317876,-0.11831248,0.17355132,-0.28164031,0.42719361]],
773    'L=32':[[0.20533892,0.32087437,0.15187897,0.04249238,0.00670516,0.00054977,0.00002018,0.00000024,0.32573501],
774        [0.12775091,-0.13523423,-0.14935701,0.28227378,0.23670434,0.05661270,0.00469819,0.32578978],
775        [0.09703829,-0.19373733,0.18610682,-0.14407046,0.00220535,0.26897090,0.36633402]],
776    'L=34':[[-0.14409234,-0.01343681,0.31248977,0.25557722,0.08571889,0.01351208,0.00095792,0.00002550,0.32573508],
777        [-0.11527834,0.18472133,-0.04403280,-0.16908618,0.27227021,0.21086614,0.04041752,0.32688152],
778        [-0.06773139,0.14120811,-0.15835721,0.18357456,-0.19364673,0.08377174,0.43116318]]
779    }
780   
781    FORPI = 12.5663706143592
782    RSQPI = 0.5641895835478
783    SQ2 = 1.414213562373
784
785    if Start:
786        ptx.pyqlmninit()
787        Start = False
788    Fln = np.zeros(len(SHCoef))
789    for i,term in enumerate(SHCoef):
790        if abs(SHCoef[term]) > 1e-6:
791            l,m,n = eval(term.strip('C'))
792            lNorm = 4.*np.pi*(2.*l+1.)
793            if SGData['SGLaue'] in ['m3','m3m']:
794                L = [i*4 for i in range(l/4)]
795                Kcl = 0.0
796                for i in L:
797                    im = i/4
798                    pcrs = ptx.pyplmpsi(l,i,phi)
799                    Kcl += BOH['L='+str(l)][n-1][l/2-1]*pcrs*cosd(i*beta)       
800            else:                #all but cubic
801                pcrs = ptx.pyplmpsi(l,n,phi)*RSQPI
802                if not n:
803                    pcrs /= SQ2
804                if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm']:
805                    if n%6 == 3:
806                        Kcl = pcrs*sind(n*beta)
807                    else:
808                        Kcl = pcrs*cosd(n*beta)
809                elif SGData['SGLaue'] in ['3mR','3m1','31m']:
810                    Kcl = pcrs*cosd(n*beta)
811                else:
812                    Kcl = pcrs*(cosd(n*beta)+sind(n*beta))
813                                   
814            Fln[i] = SHCoef[term]*Kcl*lNorm
815    ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
816    return ODFln
817   
818def polfcal(ODFln,SamSym,psi,gam):
819    RSQPI = 0.5641895835478
820    SQ2 = 1.414213562373
821    PolVal = 1.0
822    for term in ODFln:
823        if ODFln[term][1] > 1.e-3:
824            l,m,n = eval(term.strip('C'))
825            psrs = ptx.pyplmpsi(l,m,psi)
826            if SamSym in ['-1','2/m']:
827                if m:
828                    Ksl = RSQPI*psrs*(cosd(m*gam)+sind(m*gam))
829                else:
830                    ksl = RSQPI*psrs/SQ2
831            else:
832                if m:
833                    Ksl = RSQPI*psrs*cosd(m*gam)
834                else:
835                    Ksl = RSQPI*psrs/SQ2
836        PolVal += ODFln[term][1]*Ksl
837    return PolVal
838   
839# output from uctbx computed on platform darwin on 2010-05-28
840NeedTestData = True
841def TestData():
842    array = np.array
843    global NeedTestData
844    NeedTestData = False
845    global CellTestData
846    CellTestData = [
847# cell, g, G, cell*, V, V*
848  [(4, 4, 4, 90, 90, 90), 
849   array([[  1.60000000e+01,   9.79717439e-16,   9.79717439e-16],
850       [  9.79717439e-16,   1.60000000e+01,   9.79717439e-16],
851       [  9.79717439e-16,   9.79717439e-16,   1.60000000e+01]]), array([[  6.25000000e-02,   3.82702125e-18,   3.82702125e-18],
852       [  3.82702125e-18,   6.25000000e-02,   3.82702125e-18],
853       [  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],
854# cell, g, G, cell*, V, V*
855  [(4.0999999999999996, 5.2000000000000002, 6.2999999999999998, 100, 80, 130), 
856   array([[ 16.81      , -13.70423184,   4.48533243],
857       [-13.70423184,  27.04      ,  -5.6887143 ],
858       [  4.48533243,  -5.6887143 ,  39.69      ]]), array([[ 0.10206349,  0.05083339, -0.00424823],
859       [ 0.05083339,  0.06344997,  0.00334956],
860       [-0.00424823,  0.00334956,  0.02615544]]), (0.31947376387537696, 0.25189277536327803, 0.16172643497798223, 85.283666420376008, 94.716333579624006, 50.825714168082683), 100.98576357983838, 0.0099023858863968445],
861# cell, g, G, cell*, V, V*
862  [(3.5, 3.5, 6, 90, 90, 120), 
863   array([[  1.22500000e+01,  -6.12500000e+00,   1.28587914e-15],
864       [ -6.12500000e+00,   1.22500000e+01,   1.28587914e-15],
865       [  1.28587914e-15,   1.28587914e-15,   3.60000000e+01]]), array([[  1.08843537e-01,   5.44217687e-02,   3.36690552e-18],
866       [  5.44217687e-02,   1.08843537e-01,   3.36690552e-18],
867       [  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],
868  ]
869    global CoordTestData
870    CoordTestData = [
871# cell, ((frac, ortho),...)
872  ((4,4,4,90,90,90,), [
873 ((0.10000000000000001, 0.0, 0.0),(0.40000000000000002, 0.0, 0.0)),
874 ((0.0, 0.10000000000000001, 0.0),(2.4492935982947065e-17, 0.40000000000000002, 0.0)),
875 ((0.0, 0.0, 0.10000000000000001),(2.4492935982947065e-17, -2.4492935982947065e-17, 0.40000000000000002)),
876 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.40000000000000013, 0.79999999999999993, 1.2)),
877 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.80000000000000016, 1.2, 0.40000000000000002)),
878 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(1.2, 0.80000000000000004, 0.40000000000000002)),
879 ((0.5, 0.5, 0.5),(2.0, 1.9999999999999998, 2.0)),
880]),
881# cell, ((frac, ortho),...)
882  ((4.1,5.2,6.3,100,80,130,), [
883 ((0.10000000000000001, 0.0, 0.0),(0.40999999999999998, 0.0, 0.0)),
884 ((0.0, 0.10000000000000001, 0.0),(-0.33424955703700043, 0.39834311042186865, 0.0)),
885 ((0.0, 0.0, 0.10000000000000001),(0.10939835193016617, -0.051013289294572106, 0.6183281045774256)),
886 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.069695941716497567, 0.64364635296002093, 1.8549843137322766)),
887 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(-0.073350319180835066, 1.1440160419710339, 0.6183281045774256)),
888 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.67089923785616512, 0.74567293154916525, 0.6183281045774256)),
889 ((0.5, 0.5, 0.5),(0.92574397446582857, 1.7366491056364828, 3.0916405228871278)),
890]),
891# cell, ((frac, ortho),...)
892  ((3.5,3.5,6,90,90,120,), [
893 ((0.10000000000000001, 0.0, 0.0),(0.35000000000000003, 0.0, 0.0)),
894 ((0.0, 0.10000000000000001, 0.0),(-0.17499999999999993, 0.3031088913245536, 0.0)),
895 ((0.0, 0.0, 0.10000000000000001),(3.6739403974420595e-17, -3.6739403974420595e-17, 0.60000000000000009)),
896 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(2.7675166561703527e-16, 0.60621778264910708, 1.7999999999999998)),
897 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.17500000000000041, 0.90932667397366063, 0.60000000000000009)),
898 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.70000000000000018, 0.6062177826491072, 0.60000000000000009)),
899 ((0.5, 0.5, 0.5),(0.87500000000000067, 1.5155444566227676, 3.0)),
900]),
901]
902
903def test0():
904    if NeedTestData: TestData()
905    msg = 'test cell2Gmat, fillgmat, Gmat2cell'
906    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
907        G, g = cell2Gmat(cell)
908        assert np.allclose(G,tG),msg
909        assert np.allclose(g,tg),msg
910        tcell = Gmat2cell(g)
911        assert np.allclose(cell,tcell),msg
912        tcell = Gmat2cell(G)
913        assert np.allclose(tcell,trcell),msg
914
915def test1():
916    if NeedTestData: TestData()
917    msg = 'test cell2A and A2Gmat'
918    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
919        G, g = A2Gmat(cell2A(cell))
920        assert np.allclose(G,tG),msg
921        assert np.allclose(g,tg),msg
922
923def test2():
924    if NeedTestData: TestData()
925    msg = 'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
926    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
927        G, g = cell2Gmat(cell)
928        tcell = A2cell(Gmat2A(G))
929        assert np.allclose(cell,tcell),msg
930
931def test3():
932    if NeedTestData: TestData()
933    msg = 'test invcell2Gmat'
934    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
935        G, g = invcell2Gmat(trcell)
936        assert np.allclose(G,tG),msg
937        assert np.allclose(g,tg),msg
938
939def test4():
940    if NeedTestData: TestData()
941    msg = 'test calc_rVsq, calc_rV, calc_V'
942    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
943        assert np.allclose(calc_rV(cell2A(cell)),trV), msg
944        assert np.allclose(calc_V(cell2A(cell)),tV), msg
945
946def test5():
947    if NeedTestData: TestData()
948    msg = 'test A2invcell'
949    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
950        rcell = A2invcell(cell2A(cell))
951        assert np.allclose(rcell,trcell),msg
952
953def test6():
954    if NeedTestData: TestData()
955    msg = 'test cell2AB'
956    for (cell,coordlist) in CoordTestData:
957        A,B = cell2AB(cell)
958        for (frac,ortho) in coordlist:
959            to = np.inner(A,frac)
960            tf = np.inner(B,to)
961            assert np.allclose(ortho,to), msg
962            assert np.allclose(frac,tf), msg
963            to = np.sum(A*frac,axis=1)
964            tf = np.sum(B*to,axis=1)
965            assert np.allclose(ortho,to), msg
966            assert np.allclose(frac,tf), msg
967
968# test GetBraviasNum(...) and GenHBravais(...)
969def test7():
970    import os.path
971    import sys
972    import GSASIIspc as spc
973    testdir = os.path.join(os.path.split(os.path.abspath( __file__ ))[0],'testinp')
974    if os.path.exists(testdir):
975        if testdir not in sys.path: sys.path.insert(0,testdir)
976    import sgtbxlattinp
977    derror = 1e-4
978    def indexmatch(hklin, hkllist, system):
979        for hklref in hkllist:
980            hklref = list(hklref)
981            # these permutations are far from complete, but are sufficient to
982            # allow the test to complete
983            if system == 'cubic':
984                permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
985            elif system == 'monoclinic':
986                permlist = [(1,2,3),(-1,2,-3)]
987            else:
988                permlist = [(1,2,3)]
989
990            for perm in permlist:
991                hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
992                if hkl == hklref: return True
993                if [-i for i in hkl] == hklref: return True
994        else:
995            return False
996
997    for key in sgtbxlattinp.sgtbx7:
998        spdict = spc.SpcGroup(key)
999        cell = sgtbxlattinp.sgtbx7[key][0]
1000        system = spdict[1]['SGSys']
1001        center = spdict[1]['SGLatt']
1002
1003        bravcode = GetBraviasNum(center, system)
1004
1005        g2list = GenHBravais(sgtbxlattinp.dmin, bravcode, cell2A(cell))
1006
1007        assert len(sgtbxlattinp.sgtbx7[key][1]) == len(g2list), 'Reflection lists differ for %s' % key
1008        for h,k,l,d,num in g2list:
1009            for hkllist,dref in sgtbxlattinp.sgtbx7[key][1]: 
1010                if abs(d-dref) < derror:
1011                    if indexmatch((h,k,l,), hkllist, system):
1012                        break
1013            else:
1014                assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
1015
1016def test8():
1017    import GSASIIspc as spc
1018    import sgtbxlattinp
1019    derror = 1e-4
1020    dmin = sgtbxlattinp.dmin
1021
1022    def indexmatch(hklin, hklref, system, axis):
1023        # these permutations are far from complete, but are sufficient to
1024        # allow the test to complete
1025        if system == 'cubic':
1026            permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
1027        elif system == 'monoclinic' and axis=='b':
1028            permlist = [(1,2,3),(-1,2,-3)]
1029        elif system == 'monoclinic' and axis=='a':
1030            permlist = [(1,2,3),(1,-2,-3)]
1031        elif system == 'monoclinic' and axis=='c':
1032            permlist = [(1,2,3),(-1,-2,3)]
1033        elif system == 'trigonal':
1034            permlist = [(1,2,3),(2,1,3),(-1,-2,3),(-2,-1,3)]
1035        elif system == 'rhombohedral':
1036            permlist = [(1,2,3),(2,3,1),(3,1,2)]
1037        else:
1038            permlist = [(1,2,3)]
1039
1040        hklref = list(hklref)
1041        for perm in permlist:
1042            hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
1043            if hkl == hklref: return True
1044            if [-i for i in hkl] == hklref: return True
1045        return False
1046
1047    for key in sgtbxlattinp.sgtbx8:
1048        spdict = spc.SpcGroup(key)[1]
1049        cell = sgtbxlattinp.sgtbx8[key][0]
1050        center = spdict['SGLatt']
1051        Laue = spdict['SGLaue']
1052        Axis = spdict['SGUniq']
1053        system = spdict['SGSys']
1054
1055        g2list = GenHLaue(dmin,Laue,center,Axis,cell2A(cell))
1056        #if len(g2list) != len(sgtbxlattinp.sgtbx8[key][1]):
1057        #    print 'failed',key,':' ,len(g2list),'vs',len(sgtbxlattinp.sgtbx8[key][1])
1058        #    print 'GSAS-II:'
1059        #    for h,k,l,d in g2list: print '  ',(h,k,l),d
1060        #    print 'SGTBX:'
1061        #    for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print '  ',hkllist,dref
1062        assert len(g2list) == len(sgtbxlattinp.sgtbx8[key][1]), (
1063            'Reflection lists differ for %s' % key
1064            )
1065        #match = True
1066        for h,k,l,d in g2list:
1067            for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: 
1068                if abs(d-dref) < derror:
1069                    if indexmatch((h,k,l,), hkllist, system, Axis): break
1070            else:
1071                assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
1072                #match = False
1073        #if not match:
1074            #for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print '  ',hkllist,dref
1075            #print center, Laue, Axis, system
1076
1077if __name__ == '__main__':
1078    test0()
1079    test1()
1080    test2()
1081    test3()
1082    test4()
1083    test5()
1084    test6()
1085    test7()
1086    test8()
1087    print "OK"
Note: See TracBrowser for help on using the repository browser.