source: trunk/GSASIIlattice.py @ 342

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

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